% HiWinGS2013_specreader.m - loads spectra from concatenated project
% data file HiWinGS_2013_sp_v3_.txt.  These are 10-min spectra
% produced by nau_process_HiWinGS_2013.m script and concatenated
% by run_HiWinGS_2013_v080714.m

% Assumes da file has been loaded and the following variables defined:
% reldir, sig_h, sig_u, SP2, org, ww, wwj, ugw, U, vv, vvj, ugu

% data file structure:
%{
columns     variable        description
1           doy             decimal day of year
2:49        Fson_u          frequency bins
50:97       Su              U spectrum, Hare correction
98:145      Sv              V spectrum, Hare correction
146:193     Sw              W spectrum, Hare correction
194:241     St              T spectrum
242:289     So              H2O spectrum, LI7500
290:337     Suj             U spectrum, Edson correction
338:385     Svj             V spectrum, Edson correction
386:433     Swj             W spectrum, Edson correction
434:481     Cuw             UW cospectrum, Hare correction
482:529     Cvw             VW cospectrum, Hare correction
530:577     Ctw             TW cospectrum, Hare correction
578:625     Cow             H2O:W cospectrum, Hare correction
626:673     Cuwj            UW cospectrum, Edson correction
674:721     Cvwj            VW cospectrum, Edson correction
%}

%% initialize
datasetv='_v3';
decorr=0;
dcname=[];
data_path='/Volumes/User/bblomquist/Documents/Cruises/HiWinGS_2013/Analysis_BWB/procturb/total/';
if decorr==1
    dcname='_vdecorr';
end;

%% load data 
SA=load([data_path,'HiWinGS_2013_sp' datasetv dcname '.txt']); %size(SA) =  [5220  721]

lll=48;

DOY=SA(:,1);
fson=SA(:,2:lll+1);
Su=SA(:,lll+2:lll+1+48);
Sv=SA(:,lll*2+2:lll*2+1+48);
Sw=SA(:,lll*3+2:lll*3+1+48);
St=SA(:,lll*4+2:lll*4+1+48);
So=SA(:,lll*5+2:lll*5+1+48);
Suj=SA(:,lll*6+2:lll*6+1+48);
Svj=SA(:,lll*7+2:lll*7+1+48);
Swj=SA(:,lll*8+2:lll*8+1+48);
Cuw=SA(:,lll*9+2:lll*9+1+48);
Cvw=SA(:,lll*10+2:lll*10+1+48);
Ctw=SA(:,lll*11+2:lll*11+1+48);
Cow=SA(:,lll*12+2:lll*12+1+48);
Cuwj=SA(:,lll*13+2:lll*13+1+48);
Cvwj=SA(:,lll*14+2:lll*14+1+48);

%% filter using criteria and compute medians
% clear prior variables, if present
clear Summ Sumx Svmm Svmx Swmm Swmx Stmm Stmx Somm Somx 
clear Sujmm Sujmx Svjmm Svjmx Swjmm Swjmx nm dum
clear Cuwmm Cuwmx Cvwmm Cvwmx Ctwmm Ctwmx Cowmm Cowmx Cuwjmm Cuwjmx Cvwjmm Cvwjmx
clear uu

% define windspeed bins
uu(1:13)=0:12;
uu(14:21)=14:2:28;

% compute medians in windspeed bins
for i=1:length(uu)-1
    ii=find(U>=uu(i) & U<uu(i+1));  % data within windspeed bin
    % apply selection criteria
%    jj=find((reldir(ii)>300 | reldir(ii)<60) & ushp_im(ii) <5 & sig_h(ii)<4 & sig_u(ii)<0.2 &sp2(ii)<0.8 & sqrt(ww(ii))./ugw(ii)<.7+.0015*U(ii).^2 & sqrt(vv(ii))./ugu(ii)<1.4);
%    jj=find((reldir(ii)>270 | reldir(ii)<90) & sig_h(ii)<4 & sig_u(ii)<0.6 &sp2(ii)<0.8 & sqrt(ww(ii))./ugw(ii)<.8+.0015*U(ii).^2 & sqrt(vv(ii))./ugu(ii)<1.5);
%    jj=find((reldir(ii)>270 | reldir(ii)<90) & sig_h(ii)<4 & sig_u(ii)<0.6 &sp2(ii)<0.8);
    jj=find((reldir(ii)<90 | reldir(ii)>270) & sig_h(ii)<5 & sig_u(ii)<0.6 & sp2(ii)<0.8  & org(ii)<5 & sqrt(wwj(ii))./ugw(ii)<.8+.0015*U(ii).^2 & sqrt(vvj(ii))./ugu(ii)<1.5); %& isledist>3.5 );

   if isempty(jj)  | length(jj)==1  % one or fewer valid measurements
      Summ(i,:)=NaN;      Svmm(i,:)=NaN;      Swmm(i,:)=NaN;
      Stmm(i,:)=NaN;      Somm(i,:)=NaN;
      Sujmm(i,:)=NaN;     Svjmm(i,:)=NaN;     Swjmm(i,:)=NaN;
      Cuwmm(i,:)=NaN;     Cvwmm(i,:)=NaN;     Ctwmm(i,:)=NaN;
      Cowmm(i,:)=NaN;     Cuwjmm(i,:)=NaN;    Cvwjmm(i,:)=NaN;
      Sumx(i,:)=NaN;      Svmx(i,:)=NaN;      Swmx(i,:)=NaN;
      Stmx(i,:)=NaN;      Somx(i,:)=NaN;
      Sujmx(i,:)=NaN;     Svjmx(i,:)=NaN;     Swjmx(i,:)=NaN;
      Cuwmx(i,:)=NaN;     Cvwmx(i,:)=NaN;     Ctwmx(i,:)=NaN;
      Cowmx(i,:)=NaN;     Cuwjmx(i,:)=NaN;    Cvwjmx(i,:)=NaN;
            
      nm(i)=0;
      dum(i)=NaN;
   else     % compute medians
      nm(i)=length(jj);
      dum(i)=median(U(ii(jj)));
      Summ(i,:)=median(Su(ii(jj),:));
      Sumx(i,:)=nanmean(Su(ii(jj),:));
      Svmm(i,:)=median(Sv(ii(jj),:));
      Svmx(i,:)=nanmean(Sv(ii(jj),:));
      Swmm(i,:)=median(Sw(ii(jj),:));
      Swmx(i,:)=nanmean(Sw(ii(jj),:));
      Stmm(i,:)=median(St(ii(jj),:));
      Stmx(i,:)=nanmean(St(ii(jj),:));
      Somm(i,:)=median(So(ii(jj),:));
      Somx(i,:)=nanmean(So(ii(jj),:));
      Sujmm(i,:)=median(Suj(ii(jj),:));
      Sujmx(i,:)=nanmean(Suj(ii(jj),:));
      Svjmm(i,:)=median(Svj(ii(jj),:));
      Svjmx(i,:)=nanmean(Svj(ii(jj),:));
      Swjmm(i,:)=median(Swj(ii(jj),:));
      Swjmx(i,:)=nanmean(Swj(ii(jj),:));
      
      Cuwmm(i,:)=median(Cuw(ii(jj),:));   
      Cuwmx(i,:)=nanmean(Cuw(ii(jj),:));
      Cvwmm(i,:)=median(Cvw(ii(jj),:));  
      Cvwmx(i,:)=nanmean(Cvw(ii(jj),:));
      Ctwmm(i,:)=median(Ctw(ii(jj),:));
      Ctwmx(i,:)=nanmean(Ctw(ii(jj),:));
      Cowmm(i,:)=median(Cow(ii(jj),:)); 
      Cowmx(i,:)=nanmean(Cow(ii(jj),:));
      Cuwjmm(i,:)=median(Cuwj(ii(jj),:));   
      Cuwjmx(i,:)=nanmean(Cuwj(ii(jj),:));  
      Cvwjmm(i,:)=median(Cvwj(ii(jj),:));
      Cvwjmx(i,:)=nanmean(Cvwj(ii(jj),:));
   end;
end;

%% further filtering??? 
ii=find(nm>5); % more than 5 results passing criteria
% what is ia???
% jj=find((reldir(ia)>270 | reldir(ia)<90) & sig_h(ia)<5 & sig_u(ia)<1.4 &sp2(ia)<1.2 & org(ia)<5 & sqrt(ww(ia))./ugw(ia)<2+.0015*U(ia).^2 & sqrt(vv(ia))./ugu(ia)<2);

%% covariance plots
%------ wu
figure;
subplot(2,1,1)
for ll=1:length(ii)
    col=rand(3,1);
    semilogx(fson(1,:),fson(1,:).*(Cuwjmm(ii(ll),:)),'Color',col,'LineWidth',1.5);
%     legend(sprintf('%f',dum(ii(ll))))
    hold on;
end
grid;%ylim([0 0.02])
legendCell = cellstr(num2str(dum(ii)', '%.1f'));
legend(legendCell);
title(['HIWINGS2013 - Cuw averaged by wind bins']);xlabel('freqency (Hz)');ylabel('f*Cxx(f)');
subplot(2,1,2);
title(['HIWINGS2013 - Cuw averaged']);xlabel('freqency (Hz)');ylabel('f*Cxx(f)');
semilogx(fson(1,:),fson(1,:).*nanmean(Cuwj(jj,:)),'m--',fson(1,:),fson(1,:).*median(Cuwj(jj,:)),'g--','LineWidth',3);grid
legend('mean','median')
orient tall

%------ wv
figure;
subplot(2,1,1)
for ll=1:length(ii)
    col=rand(3,1);
    semilogx(fson(1,:),fson(1,:).*(Cvwjmm(ii(ll),:)),'Color',col,'LineWidth',1.5);
%     legend(sprintf('%f',dum(ii(ll))))
    hold on;
end
grid;%ylim([0 0.02])
legendCell = cellstr(num2str(dum(ii)', '%.1f'));
legend(legendCell);
title(['HIWINGS2013 - Cvw averaged by wind bins']);xlabel('freqency (Hz)');ylabel('f*Cxx(f)');
subplot(2,1,2);
title(['HIWINGS2013 - Cvw averaged']);xlabel('freqency (Hz)');ylabel('f*Cxx(f)');
semilogx(fson(1,:),fson(1,:).*nanmean(Cvwj(jj,:)),'m--',fson(1,:),fson(1,:).*median(Cvwj(jj,:)),'g--','LineWidth',3);grid
legend('mean','median')

%------ wT
figure;
subplot(2,1,1)
for ll=1:length(ii)
    col=rand(3,1);
    semilogx(fson(1,:),fson(1,:).*(Ctwmm(ii(ll),:)),'Color',col,'LineWidth',1.5);
%     legend(sprintf('%f',dum(ii(ll))))
    hold on;
end
grid;ylim([0 0.02])
legendCell = cellstr(num2str(dum(ii)', '%.1f'));
legend(legendCell);
title(['HIWINGS2013 - Cwt averaged by wind bins']);xlabel('freqency (Hz)');ylabel('f*Cxx(f)');
subplot(2,1,2);
title(['HIWINGS2013 - Cwt averaged']);xlabel('freqency (Hz)');ylabel('f*Cxx(f)');
semilogx(fson(1,:),fson(1,:).*nanmean(Ctw(jj,:)),'m--',fson(1,:),fson(1,:).*median(Ctw(jj,:)),'g--','LineWidth',3);grid
legend('mean','median')
orient tall

%------ wT
figure;
set(gcf,'defaultaxesfonts',12)
set(gca,'fontname','time','fontsize',12);
loglog(fson(1,:),fson(1,:).*median(Ctw(jj,:)),fson(1,:),fson(1,:).^(-4/3)*median(Ctw(jj,36)),'--','LineWidth',2);grid
title(['HIWINGS2013 - C_{wT} averaged']);xlabel('freqency (Hz)');ylabel('f*C_{wT}(f)');
legend('C_{wT}','-4/3');ylim([1e-5 1e-2])

%------ T spectrum
figure;
set(gcf,'defaultaxesfonts',12)
set(gca,'fontname','time','fontsize',12);
loglog(fson(1,:),fson(1,:).*median(St(jj,:)),fson(1,:),fson(1,:).^(-2/3)*median(St(jj,36)),'--','LineWidth',2);grid
title(['HIWINGS2013 - S_T averaged']);xlabel('freqency (Hz)');ylabel('f*S_T(f)');
legend('S_T','-2/3');ylim([8e-4 1e-1])
ylim([1e-4 2e-3])

%------ U,V,W spectra
figure;
set(gcf,'defaultaxesfonts',12)
set(gca,'fontname','time','fontsize',12);
loglog(fson(1,:),fson(1,:).*median(Suj(jj,:)),fson(1,:),fson(1,:).*median(Svj(jj,:)),fson(1,:),fson(1,:).*median(Swj(jj,:)),fson(1,:),fson(1,:).^(-2/3)*median(Sw(jj,36)),'--','LineWidth',2);grid
legend('S_U','S_V','S_W','-2/3');%ylim([8e-4 1e-1])
title(['HIWINGS2013 - S_{x} averaged']);xlabel('freqency (Hz)');ylabel('f*S_{x}(f)');

