%try to do spectrogram of Hs
clear all; fclose all;


%try #2 - interpolate to 36000 and then computes Spectra
Spectro_freq1=ones(31*24,62)*NaN;
Spectro_Hs1=ones(31*24,62)*NaN;
Spectro_freq2=ones(31*24,62)*NaN;
Spectro_Hs2=ones(31*24,62)*NaN;
doy_all=ones(31*24,1)*NaN;
lHs_all=ones(31*24,1)*NaN;
Hs_stdispline=ones(31*24,1)*NaN;
Hs_stdiinterp=ones(31*24,1)*NaN;
Hs_std=ones(31*24,1)*NaN;
lFHs1=ones(31*24,1)*NaN;
lFHs2=ones(31*24,1)*NaN;
for ddd=311:341
    jd=sprintf('%03i',ddd);
    for hhh=00:23
        ic=(ddd-311)*24+hhh+1;
        [ddd hhh ic]
        hr=sprintf('%02i',hhh);
        if exist(['C:\Data\CRUISES\DYNAMO_2011\data\Riegl_corrected\Riegl_heavesplineinterpolatedtowave_v012012_leg3\Wavecorrected_',jd hr,'.prn'])
            Hs_raw=load(['C:\Data\CRUISES\DYNAMO_2011\data\Riegl_corrected\Riegl_heavesplineinterpolatedtowave_v012012_leg3\Wavecorrected_',jd hr,'.prn']);
            Hs=Hs_raw(:,end);
            Hs_std(ic)=nanstd(Hs);
            lHs_all(ic)=length(Hs(~isnan(Hs)));lHs_all(ic)
            dref=Hs_raw(:,1)+Hs_raw(:,2)/24+Hs_raw(:,3)/3600/24;
            
            %using spline
            Hs_spline=interp1(dref(~isnan(Hs)),Hs(~isnan(Hs)),dref,'spline');              Hs_spline(end-10:end)=Hs(end-10);            Hs_spline(1:10)=Hs(10);
            Hs_stdispline(ic)=nanstd(Hs_spline);
                [SHs,FHs]=psd2(detrend(Hs_spline(~isnan(Hs_spline))),length(Hs_spline(~isnan(Hs_spline))),10,hamming(length(Hs_spline(~isnan(Hs_spline)))));[SHs1,FHs1]=specsmoo(SHs,10); 
                Spectro_freq1(ic,1:length(SHs1))=FHs1;  
                Spectro_Hs1(ic,1:length(SHs1))=SHs1;
            lFHs1(ic)=length(FHs1);  
            %using nearest and patching NaN with nanmedian
            Hs_interp=interp1(dref(~isnan(Hs)),Hs(~isnan(Hs)),dref);              Hs_interp(isnan(Hs_interp))=nanmedian(Hs_interp);
            Hs_stdiinterp(ic)=nanstd(Hs_interp);
                [SHs,FHs]=psd2(detrend(Hs_interp(~isnan(Hs_interp))),length(Hs_interp(~isnan(Hs_interp))),10,hamming(length(Hs_interp(~isnan(Hs_interp)))));[SHs2,FHs2]=specsmoo(SHs,10); 
                Spectro_freq2(ic,1:length(SHs1))=FHs2;  
                Spectro_Hs2(ic,1:length(SHs1))=SHs2;
           lFHs2(ic)=length(FHs2); 
            
    %         doy_all(ic)=Hs_raw(:,1)+Hs_raw(:,2)/24+Hs_raw(:,3)/3600/24;
            doy_all(ic)=ddd+hhh/24;
        end    
    end
end    

Spectro_freq1(isnan(doy_all),:)=[];Spectro_Hs1(isnan(doy_all),:)=[];
Spectro_freq2(isnan(doy_all),:)=[];Spectro_Hs2(isnan(doy_all),:)=[];
lHs_all(isnan(doy_all))=[];
Hs_stdispline(isnan(doy_all))=[];
Hs_stdiinterp(isnan(doy_all))=[];
Hs_std(isnan(doy_all))=[];

doy_all(isnan(doy_all))=[];


% figure;imagesc(doy_all,Spectro_freq2(1,:), real(log10(Spectro_Hs2)'),[-6 1]);  %
% figure;pcolor(doy_all,Spectro_freq2(1,:), real(log10(Spectro_Hs2)'));  %
% figure;surf(doy_all,Spectro_freq2(1,:),zeros(size(doy_all)),real(log10(Spectro_Hs2)')); view([0 90])
figure(1);pcolor(doy_all,Spectro_freq1(1,:), real(log10(Spectro_Hs1)'));  %
set(gca,'CLim',[-5 0.5])
shading(gca,'interp')
h = colorbar;

figure(2);pcolor(doy_all,Spectro_freq2(1,:), real(log10(Spectro_Hs2)'));  %
set(gca,'CLim',[-5 0.5])
shading(gca,'interp')
h = colorbar;

% stdt=floor(doy_all(1)+datenum(2011,1,0));
% endt=ceil(doy_all(end)+datenum(2011,1,0)); 
% stdt=floor(doy_all(1)-304);
% endt=ceil(doy_all(end)-304); 
stdt=floor(doy_all(1));
endt=ceil(doy_all(end));
xs=2;
% figure(3);pcolor(doy_all+datenum(2011,1,0),Spectro_freq2(1,:), real(log10(Spectro_Hs2)'));  %
% figure(3);pcolor(doy_all-304,Spectro_freq2(1,:), real(log10(Spectro_Hs2)'));  %
figure(3);pcolor(doy_all+1,Spectro_freq2(2,:), real(log10(Spectro_Hs2)'));  %
set(gca,'CLim',[-5 0.5])
shading(gca,'interp')
h = colorbar;
set(get(h,'ylabel'),'String', 'log_{10} \psi (m^2/s)');
% colormap('Gray')
% axis([326 330 0 1])
cmap = b2rcolormap(real(log10(Spectro_Hs2)'));
xlabel('Year (2011)');ylabel('Frequency (Hz)')
set(gca(gcf),'XTick',stdt:xs:endt);datetick('x','mm/dd','keepticks','keeplimits');
set(gca(gcf),'TickDir','out') 

figure;pcolor(doy_all(oo)+1,Spectro_freq1(2,:), real(log10(Spectro_Hs1(oo,:))'));  %
set(gca,'CLim',[-5 0.5])
shading(gca,'interp')
h = colorbar;

% oo=find((Hs_stdispline-Hs_std>1) | lHs_all<6000);
% test=Spectro_Hs1;test(oo,:)=NaN;
oo=find((Hs_stdiinterp-Hs_std>.1));
test=Spectro_Hs2;test(oo,:)=NaN;
dtest=doy_all;%dtest(oo)=[];

figure;pcolor(dtest+1,Spectro_freq1(2,:), real(log10(test)'));  %
set(gca,'CLim',[-5 0.5])
shading(gca,'interp')
h = colorbar;
cmap = b2rcolormap(real(log10(Spectro_Hs2)'));




Hs_heave2wavespline_leg3=load('C:\Data\CRUISES\DYNAMO_2011\data\Riegl_corrected\Riegl_heavesplineinterpolatedtowave_v012012_leg3\Hs_hourly_leg3.prn');
doy_leg3=Hs_heave2wavespline_leg3(:,1)+Hs_heave2wavespline_leg3(:,2)/24;
ii=find(abs(Hs_heave2wavespline_leg3(:,7))<0.06  & abs(Hs_heave2wavespline_leg3(:,5))<1 & Hs_heave2wavespline_leg3(:,6) >10);

badpts=find(abs(doy_leg3-324.4)==min(abs(doy_leg3-324.4)))-1;
Hs_heave2wavespline_leg3(badpts,:)=NaN;test(badpts,:)=ones(1,62)*NaN;
badpts=find(abs(doy_leg3-315.3)==min(abs(doy_leg3-315.3)));
Hs_heave2wavespline_leg3(badpts,:)=NaN;test(badpts,:)=ones(1,62)*NaN;
figure;plot(doy_leg3(ii),Hs_heave2wavespline_leg3(ii,3),'b.--',doy_leg3(ii),Hs_heave2wavespline_leg3(ii,4),'r-','LineWidth',2);grid;

doy_allo=repmat(doy_all,1,62);


% doy=Hs_raw(:,1)+Hs_raw(:,2)/24+Hs_raw(:,3)/3600/24;
% 
% ii=find(~isnan(Hs));
% Hs=Hs_raw(:,end);
% Hsi=interp1(doy(ii),Hs(ii),doy,'spline'); Hsi(end-20:end)=(Hsi(end-21)); Hsi(1:10)=Hsi(11);
% figure;plot(doy,Hs,'.',doy,Hsi,'r--');
% 
% [SHs,FHs]=psd2(detrend(Hs(~isnan(Hs))),length(Hs(~isnan(Hs))),10,hamming(length(Hs(~isnan(Hs)))));[SHs1,FHs1]=specsmoo(SHs,10); 
% [SHsi,FHsi]=psd2(detrend(Hsi(~isnan(Hsi))),length(Hsi(~isnan(Hsi))),10,hamming(length(Hsi(~isnan(Hsi)))));[SHsi1,FHsi1]=specsmoo(SHsi,10); 
% 
% figure;clf;loglog(FHs1,FHs1.*SHs1,'b',FHs1,FHs1.*SHs1,'r')

Spectro_freqi=Spectro_freq(38,:);
Spectro_Hsi=ones(31*24,62)*NaN;
for oo=1:length(Spectro_freq)
    try
        Spectro_Hsi(oo,:)=interp1(Spectro_freq(oo,~isnan(Spectro_freq(oo,:))),Spectro_Hs(oo,~isnan(Spectro_freq(oo,:))),Spectro_freqi,'spline');
    catch
        Spectro_Hsi(oo,:)=ones(1,62)*NaN;
    end    
end    

% doy_all(isnan(Spectro_Hsi(:,end)))=[];
Spectro_Hsi(isnan(Spectro_Hsi(:,end)),:)=[];



stdt=floor(doy_all(1)+datenum(2011,1,0));
endt=ceil(doy_all(end)+datenum(2011,1,0)); 
xs=2;
% figure;imagesc(doy_all,Spectro_freqi, real(log10(Spectro_Hsi)),[-6 -1]);  %
figure;imagesc(doy_all+datenum(2011,1,0),Spectro_freqi, real(log10(Spectro_Hsi)),[-3 1]);  %
axis xy;ylim([0 .5]);
xlabel('November (2011)');ylabel('log_{10}Hs(m^2/s)')
set(gca(gcf),'XTick',stdt:xs:endt);datetick('x','mm/dd','keepticks','keeplimits');
h = colorbar;
axes(h)


lud=find(floor(doy_all)==328);
figure;loglog(Spectro_freq(1,:), Spectro_Hs(lud(1),:),Spectro_freqi,Spectro_Hsi(lud(1),:),'r')

        