%look at PSD for when the Rielg agree and disagree with WR and spar...

% rge=292;
rge=288;

Spectro_freq2=ones(length(rge)+1*13*6,62)*NaN;
Spectro_Hs2=ones(length(rge)+1*13*6,62)*NaN;
Spectro_freq1=ones(length(rge)+1*13*6,62)*NaN;
Spectro_Hs1=ones(length(rge)+1*13*6,62)*NaN;
Hs_fromspectra=ones(length(rge)+1*13*6,1)*NaN;
doy_all=ones(length(rge)+1*13*6,1)*NaN;
lHs_all=ones(length(rge)+1*13*6,1)*NaN;
Hs_stdiinterp=ones(length(rge)+1*13*6,1)*NaN;
Hs_std=ones(length(rge)+1*13*6,1)*NaN;
lFHs2=ones(length(rge)+1*13*6,1)*NaN;
for ddd=rge
    jd=sprintf('%03i',ddd);
    for hhh=12%:24;%00:12
        ic=1;
        [ddd hhh ic]
        hr=sprintf('%02i',hhh);
        if exist(['D:\Desktop\Data\CRUISES\HiWinGS_2013\flux_data\Riegl_corrected\Riegl_waveplineinterpolatedtoheave_minusriegldplusheave_vcfrp\Wavecorrected_',jd hr,'.prn'])
            Hs_raw=load(['D:\Desktop\Data\CRUISES\HiWinGS_2013\flux_data\Riegl_corrected\Riegl_waveplineinterpolatedtoheave_minusriegldplusheave_vcfrp\Wavecorrected_',jd hr,'.prn']);
            Hs=Hs_raw(:,end);
            
        %removes first and last 30s for motcorr filter tapering issues
        Hs(1:300)=NaN;        Hs(end-299:end)=NaN;
             
            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 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_raw=reshape(Hs_interp,6000,6);
            
            for lu=1:6  %computes on 10min basis                          
                Hs_stdiinterp(ic)=nanstd(Hs_raw(:,lu));
                [SHs,FHs]=psd2(detrend(Hs_raw(:,lu)),length(Hs_raw(:,lu)),10,hamming(length(Hs_raw(:,lu))));[SHs2,FHs2]=specsmoo(SHs,10); 
                Spectro_freq2(ic,1:length(SHs2))=FHs2;  
                Spectro_Hs2(ic,1:length(SHs2))=SHs2;
                lFHs2(ic)=length(FHs2); 
            
                [Sl, Fl, Df]=specsmoo4wave(SHs,10);
                Spectro_freq1(ic,1:length(Sl))=Fl;  
                Spectro_Hs1(ic,1:length(Sl))=Sl;
%                 [Hs,Tz,Te,Tpc]=Wavestatsfromspectra(Sl,Fl,Df);
                [Hs,Tz,Te,Tpc]=Wavestatsfromspectra(SHs,FHs,diff(FHs(1:2)))
                Hs_fromspectra(ic)=Hs;
                                
    %         doy_all(ic)=Hs_raw(:,1)+Hs_raw(:,2)/24+Hs_raw(:,3)/3600/24;
                doy_all(ic)=ddd+hhh/24+10*(lu-1)/24/60;
                ic=ic+1;
            end    
        end    
    end
end    

Spectro_freq2(isnan(doy_all),:)=[];Spectro_Hs2(isnan(doy_all),:)=[];
lHs_all(isnan(doy_all))=[];
Hs_stdiinterp(isnan(doy_all))=[];
Hs_std(isnan(doy_all))=[];
Hs_fromspectra(isnan(doy_all))=[];
doy_all(isnan(doy_all))=[];

% Spectro_Hs2(:,Spectro_freq2(1,:)>1.2)=NaN;Spectro_freq2(Spectro_freq2(1,:)>1.2,:)=NaN;


[l,c]=size(Spectro_Hs2);
figure(1);clf;
set(gcf,'defaultaxesfonts',12)
set(gca,'fontname','time','fontsize',12);
for ll=1:l
    loglog(Spectro_freq2(ll,:),Spectro_freq2(ll,:).*Spectro_Hs2(ll,:),'b');hold on
    loglog(Spectro_freq1(ll,:),Spectro_freq1(ll,:).*Spectro_Hs2(ll,:),'c');hold on
end;
hold on;loglog(Spectro_freq2(ll,:),Spectro_freq2(ll,:).*nanmedianx(Spectro_Hs2),'r','LineWidth',2);
hold on;loglog(Spectro_freq2(1,end-25:end),Spectro_freq2(1,end-25:end).^(-4)*Spectro_Hs2(1,end-11),'k--','LineWidth',2);
xlabel('Frequency (Hz)');ylabel('f*Sw(f)  (m^2)')
grid
title(sprintf('Wave Height PSD (%03i)',ddd))

waverider=load('D:\Desktop\Data\CRUISES\HiWinGS_2013\flux_data\Riegl_corrected\waverider\WR.stats2.mat');
ii=find(waverider.stats2.doy>=doy_all(1) & waverider.stats2.doy<=doy_all(end));
for ll=1:length(ii)
    loglog(waverider.stats2.f(:,1),waverider.stats2.f(:,1).*nanmedianx(waverider.stats2.S(:,ii(ll))')','k');
end;
hold on;loglog(waverider.stats2.f(:,1),waverider.stats2.f(:,1).*nanmedianx(waverider.stats2.S(:,ii)')','g','LineWidth',2);


figure(2);clf;
plot(doy_all,4*Hs_stdiinterp,'b',doy_all,Hs_fromspectra,'g')
hold on;plot(waverider.stats2.doy(ii),waverider.stats2.Hs(ii),'m')


% % 
% % 
% % 
% % 
% % 
% % [l,c]=size(Spectro_Hs2);
% figure(2);clf;
% set(gcf,'defaultaxesfonts',12)
% set(gca,'fontname','time','fontsize',12);
% for ll=1:l
%     loglog(Spectro_freq2(ll,:),Spectro_Hs2(ll,:),'b');hold on
% end;
% hold on;loglog(Spectro_freq2(ll,:),nanmedianx(Spectro_Hs2),'r','LineWidth',2);
% hold on;loglog(Spectro_freq2(1,end-25:end),Spectro_freq2(1,end-25:end).^(-5)*Spectro_Hs2(1,end-11),'k--','LineWidth',2);
% xlabel('Frequency (Hz)');ylabel('Sw(f) (m^2 s)')
% grid
% title(sprintf('Wave Height PSD (%03i)',ddd))
% % 
% % % RIEGLmat=[doy_all Spectro_freq2 Spectro_Hs2];
% % % sdg=['save ', 'D:\Desktop\Data\CRUISES\HiWinGS_2013\flux_data\Riegl_corrected\Wavecorrected_',jd hr,'.prn zapmat -ascii'];
% % % eval(sdg)
% % save(['D:\Desktop\Data\CRUISES\HiWinGS_2013\flux_data\Riegl_corrected\HiWinGS2013_RIEGL_doy' jd '.mat'], 'doy_all','Spectro_freq2','Spectro_Hs2');
% % 
% % f=Spectro_freq2(1,:);
% % SHs=nanmedianx(Spectro_Hs2);
% % save(['D:\Desktop\Data\CRUISES\HiWinGS_2013\flux_data\Riegl_corrected\HiWinGS2013_RIEGL_doy' jd '.mat'], 'doy_all','Spectro_freq2','SHs');
% % 
% % % figure;
% % % waveriderHW3=load('D:\Desktop\Data\CRUISES\HiWinGS_2013\flux_data\Riegl_corrected\waverider\HW3_buoy_data');
% % % ii=find(waveriderHW3.STATS.doy>=doy_all(1) & waveriderHW3.STATS.doy<=doy_all(end));
% % % for ll=1:length(ii)
% % %     loglog(waveriderHW2.SPEC.f,waveriderHW2.SPEC.f.*waveriderHW2.SPEC.S(ll,:),'b');hold on
% % % end;
% % % hold on;loglog(waveriderHW2.SPEC.f,waveriderHW2.SPEC.f.*nanmedianx(waveriderHW2.SPEC.S(ll,:)),'g','LineWidth',2);
