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

rge=283:301;

%try #2 - interpolate to 36000 and then computes Spectra
Spectro_freq1=ones(length(rge)+1*24,62)*NaN;
Spectro_Hs1=ones(length(rge)+1*24,62)*NaN;
Spectro_freq2=ones(length(rge)+1*24,62)*NaN;
Spectro_Hs2=ones(length(rge)+1*24,62)*NaN;
doy_all=ones(length(rge)+1*24,1)*NaN;
lHs_all=ones(length(rge)+1*24,1)*NaN;
Hs_stdispline=ones(length(rge)+1*24,1)*NaN;
Hs_stdiinterp=ones(length(rge)+1*24,1)*NaN;
Hs_std=ones(length(rge)+1*24,1)*NaN;
lFHs1=ones(length(rge)+1*24,1)*NaN;
lFHs2=ones(length(rge)+1*24,1)*NaN;
for ddd=283:301
    jd=sprintf('%03i',ddd);
    for hhh=00:23
        ic=(ddd-283)*24+hhh+1;
        [ddd hhh ic]
        hr=sprintf('%02i',hhh);
        if exist(['C:\data\CRUISES\HiWinGS_2013\flux_data\Riegl_corrected\Riegl_heavesplineinterpolatedtowave_minusriegldplusheave_vcfrp2\Wavecorrected_',jd hr,'.prn'])
            Hs_raw=load(['C:\data\CRUISES\HiWinGS_2013\flux_data\Riegl_corrected\Riegl_heavesplineinterpolatedtowave_minusriegldplusheave_vcfrp2\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(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;



%%%% do subplots
% for ddd=283:295
%     jd=sprintf('%03i',ddd);
%     filecat(['C:\Data\CRUISES\HiWinGS_2013\flux_data\Riegl_corrected\Riegl_waveplineinterpolatedtoheave_minusriegldplusheave_vcfrp\Hs_' jd '*.prn'],['Hs_doy' jd '.txt']);
% end
% filecat(['C:\Data\CRUISES\HiWinGS_2013\flux_data\Riegl_corrected\Riegl_waveplineinterpolatedtoheave_minusriegldplusheave_vcfrp\Hs_doy*.txt'],['Hs_hourly_all.prn']);
Hs_heave2wavespline_all=load('C:\data\CRUISES\HiWinGS_2013\flux_data\Riegl_corrected\Riegl_heavesplineinterpolatedtowave_minusriegldplusheave_vcfrp2\Hs_hourly_all.prn');
doy_all=Hs_heave2wavespline_all(:,1)+Hs_heave2wavespline_all(:,2)/24;
% % ii=find(abs(Hs_heave2wavespline_leg3(:,7))<0.1);   %finds nanmean(wave)>0.1
% ii=find(abs(Hs_heave2wavespline_all(:,7))<0.06  & abs(Hs_heave2wavespline_leg3(:,5))<1 & Hs_heave2wavespline_leg3(:,6) >10);
ii=find(abs(Hs_heave2wavespline_all(:,4))<10 );

x=load(['C:\Data\CRUISES\HiWinGS_2013\Knorr\flux\Processed\HiWinGS_2013_PSD_flux_5min_all.txt']);
jdy=x(:,1);
% jdy=x(:,1)+datenum(2013,1,0);%julian day to Matlab format
% jdy=jdy+datenum(Vdate(1),1,0);  %converts time to Serial Date Number
U=x(:,2);%true wind,PSD sonic (m/s)
% stdt=floor(jdy(1));
stdt=283;
endt=ceil(doy_all(end)); 
xs=2;  %defines spacing for plot on x-axis
cruise='HiWinGS_2013';

figure(2);clf;
set(gcf,'defaultaxesfonts',12)
set(gca,'fontname','time','fontsize',12);

    subplot(3,1,2);pcolor(doy_all,Spectro_freq1(1,:), real(log10(Spectro_Hs1)'));  %
    ylim([0 0.5]);xlim([stdt endt]);
    title(sprintf('Wave spectrogram'),'FontWeight','Bold','Interpreter','none')
set(gca,'CLim',[-5 2])
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_Hs1)'));
% xlabel('Year (2011)');
ylabel('Frequency (Hz)')
set(gca(gcf),'XTick',stdt:xs:endt);datetick('x','mm/dd','keepticks','keeplimits');
set(gca(gcf),'TickDir','out') 
 h2=get(gca,'position'); 
 
 hold on;plot([284 284],[0 0.5],'k--',[286 286],[0 0.5],'k--','LineWidth',2)
 hold on;plot([287.2 287.2],[0 0.5],'k--',[289.85 289.85],[0 0.5],'k--','LineWidth',2)
 hold on;plot([290.7 290.7],[0 0.5],'k--',[294.4 294.4],[0 0.5],'k--','LineWidth',2)

    subplot(3,1,1);plot(jdy,U,'b','LineWidth',2);
xlim([stdt endt]);ylabel('WS (ms^{-1})');ylim([0 30])
title(sprintf('5-min wind speed'),'FontWeight','Bold','Interpreter','none')
set(gca(gcf),'XTick',stdt:xs:endt);datetick('x','mm/dd','keepticks','keeplimits');
grid on;%grid minor
        h1=get(gca,'position');h1(3)=h2(3); %sets same width for figure
 set(gca,'position',h1)
 
  hold on;plot([284 284],[0 30],'k--',[286 286],[0 30],'k--','LineWidth',2)
 hold on;plot([287.2 287.2],[0 30],'k--',[289.85 289.85],[0 30],'k--','LineWidth',2)
 hold on;plot([290.7 290.7],[0 30],'k--',[294.4 294.4],[0 30],'k--','LineWidth',2)


%     subplot(2,1,2);plot(doy_leg3(ii),Hs_heave2wavespline_leg3(ii,3),'b.--',doy_leg3(ii),Hs_heave2wavespline_leg3(ii,4),'r-','LineWidth',2)
    subplot(3,1,3);plot(doy_all(:),Hs_heave2wavespline_all(:,3),'b.--',doy_all(ii),Hs_heave2wavespline_all(ii,4),'r-','LineWidth',2)
legend('Uncorrected', 'Corrected','location','South');
ylabel('Hs (m)');xlabel('Year (2013)');
title(sprintf('Hourly significant wave height'),'FontWeight','Bold','Interpreter','none')
axis([stdt endt 0 15]);
set(gca(gcf),'XTick',stdt:xs:endt);datetick('x','mm/dd','keepticks','keeplimits');
grid on;%grid minor
        h3=get(gca,'position');h3(3)=h2(3); %sets same width for figure
 set(gca,'position',h3)
  hold on;plot([284 284],[0 15],'k--',[286 286],[0 15],'k--','LineWidth',2)
 hold on;plot([287.2 287.2],[0 15],'k--',[289.85 289.85],[0 15],'k--','LineWidth',2)
 hold on;plot([290.7 290.7],[0 15],'k--',[294.4 294.4],[0 15],'k--','LineWidth',2)

 orient tall
% if printit==1
%     print('-djpeg90 ',['C:\data\CRUISES\HiWinGS_2013\flux_data\Riegl_corrected\Riegl_heavesplineinterpolatedtowave_minusriegldplusheave_vcfrp2\HiWinGS2013_spectrogram_to1Hz_withHsandWS_v2.jpg']);
% end; 






waveriderHW1=load('C:\Data\CRUISES\HiWinGS_2013\flux_data\Riegl_corrected\waverider\HW1a_buoy_data');
waveriderHW2=load('C:\Data\CRUISES\HiWinGS_2013\flux_data\Riegl_corrected\waverider\HW2_buoy_data');
waveriderHW3=load('C:\Data\CRUISES\HiWinGS_2013\flux_data\Riegl_corrected\waverider\HW3_buoy_data');
figure;plot(doy_all(:),Hs_heave2wavespline_all(:,3),'b.--',doy_all(ii),Hs_heave2wavespline_all(ii,4),'r-','LineWidth',2)
hold on;plot(waveriderHW1.STATS.doy,waveriderHW1.STATS.Hs,'k');
hold on;plot(waveriderHW2.STATS.doy,waveriderHW2.STATS.Hs,'k');
hold on;plot(waveriderHW3.STATS.doy,waveriderHW3.STATS.Hs,'k');
legend('Riegl Uncorrected', 'Riegl Corrected','WaveRider','location','South');
ylabel('Hs (m)');xlabel('Year (2013)');
title(sprintf('Hourly significant wave height'),'FontWeight','Bold','Interpreter','none')
axis([stdt endt 0 10]);
set(gca(gcf),'XTick',stdt:xs:endt);datetick('x','mm/dd','keepticks','keeplimits');
grid on;%grid minor
        h3=get(gca,'position');h3(3)=h2(3); %sets same width for figure
 set(gca,'position',h3)
  hold on;plot([284 284],[0 10],'k--',[286 286],[0 10],'k--','LineWidth',2)
 hold on;plot([287.2 287.2],[0 10],'k--',[289.85 289.85],[0 10],'k--','LineWidth',2)
 hold on;plot([290.7 290.7],[0 10],'k--',[294.4 294.4],[0 10],'k--','LineWidth',2)

