% loads & plots sp files for WFIP2 flux tower PS01

% sp file format:
%{
	30-min sp array of spectra & cospectra: N rows x 1179 columns

	column designations as follows, all spectra are 62 bin

    COLUMNS     VARIABLE        DESCRIPTION
    1           doy             decimal day of year
    2:63        Fx              frequency array for all spectra/cospectra
    64:125      Pu03            U03 spectrum
    126:187     Pv03            V03 spectrum
    188:249     Pw03            W03 spectrum
    250:311     Pt03            T03 spectrum
    312:373     Pq03            q03 spectrum, LI7500
    374:435     Cuw03           wu03 cospectrum
    436:497     Cvw03           wv03 cospectrum
    498:559     Ctw03           wt03 cospectrum
    560:621     Cqw03           wq03 cospectrum
    622:683     Pu10            U10 spectrum
    684:745     Pv10            V10 spectrum
    746:807     Pw10            W10 spectrum
    808:869     Pt10            T10 spectrum
    870:931     Cuw10           wu10 cospectrum
    932:993     Cvw10           wv10 cospectrum
    994:1055    Ctw10           wt10 cospectrum
    1056:1117   Pco203          co2 spectrum (licor)
    1118:1179   Cco2w03         wco2 cospectrum (licor)
%}

%% initialize
clear;
close all;

% matlab script path
path_prog = fullfile('/Users/byronb/Documents/MATLAB/2016_WFIP2');
cd(path_prog);                  % cd to matlab directory
restoredefaultpath;             % reset matlab path
addpath(genpath(path_prog));    % add scripts to path
rehash;

%% load data

% load giant spectra array: 14976 x 1179
% SP = load('/Users/byronb/Raw_Data_Files.noindex/WFIP2_2016/PS01/processed_final/sp_PS01_all_PF.txt');
load ('/Users/byronb/Raw_Data_Files.noindex/WFIP2_2016/PS01/processed_final/sp_PS01_all_PF.mat');
% load all other variables: 139 arrays of 14976x1
load ('/Users/byronb/Raw_Data_Files.noindex/WFIP2_2016/PS01/processed_final/da_PS01_all_PF_rad_soil.mat');


%% define individual spectra arrays, all are 14976 x 62

fson = SP(1,2:63);

Su03 = SP(:,64:125);
Sv03 = SP(:,126:187);
Sw03 = SP(:,188:249);
St03 = SP(:,250:311);
Sq03 = SP(:,312:373);
Cuw03 = SP(:,374:435);
Cvw03 = SP(:,436:497);
Ctw03 = SP(:,498:559);
Cqw03 = SP(:,560:621);
Pco203 = SP(:,1056:1117);
Cco2w03 = SP(:,1118:1179);

Su10 = SP(:,622:683);
Sv10 = SP(:,684:745);
Sw10 = SP(:,746:807);
St10 = SP(:,808:869);
Cuw10 = SP(:,870:931);
Cvw10 = SP(:,932:993);
Ctw10 = SP(:,994:1055);

clear SP  % clear up a bit of memory

%% plot mean spectra based on following criteria

% define a time range of interest
kk = find(dnum>=dnum(1) & dnum<=dnum(end));     % time range

% define various criteria to filter spectra before averaging
jj = find(agc_max<60);                  % licor criteria
ll = find(wdir10>0 & wdir10<360);       %  & other limits

yy = intersect(ll,kk);                  % general filter
zz = intersect(yy,jj);                  % general + licor filter

% or, use predefined logical arrays
yy = good03 & good10;
zz = yy & good_licor;

% Date String
Dstr_yy = [datestr(dnum(yy(1)),'mm/dd/yy, HHMM'),' to ',datestr(dnum(yy(end)),'mm/dd/yy, HHMM')];
Dstr_zz = [datestr(dnum(zz(1)),'mm/dd/yy, HHMM'),' to ',datestr(dnum(zz(end)),'mm/dd/yy, HHMM')];

close all

% compute mean spectra/cospectra and plot
figure;
loglog(fson,fson.*nanmean1(Su03(yy,:),1),'b-',fson,fson.*nanmean1(Su10(yy,:),1),'r-','linewidth',2); grid;
title(['WFIP2 PS01 U Spectrum, ',Dstr_yy]);
ylabel('S_{uu}(f) x f'); xlabel('f (Hz)');
legend('3m','10m','location','best');

figure;
loglog(fson,fson.*nanmean1(Sw03(yy,:),1),'b-',fson,fson.*nanmean1(Sw10(yy,:),1),'r-','linewidth',2); grid;
title(['WFIP2 PS01 W Spectrum, ',Dstr_yy]);
ylabel('S_{ww}(f) x f'); xlabel('f (Hz)');
legend('3m','10m','location','best');

figure;
loglog(fson,fson.*nanmean1(St03(yy,:),1),'b-',fson,fson.*nanmean1(St10(yy,:),1),'r-','linewidth',2); grid;
title(['WFIP2 T Spectrum, ',Dstr_yy]);
ylabel('S_{tt}(f) x f'); xlabel('f (Hz)');
legend('3m','10m','location','best');

figure;
loglog(fson,fson.*nanmean1(Sq03(zz,:),1),'b-','linewidth',2); grid;
title(['WFIP2 PS01 q Spectrum, ',Dstr_zz]);
ylabel('S_{qq}(f) x f'); xlabel('f (Hz)');

figure;
loglog(fson,fson.*nanmean1(Pco203(zz,:),1),'b-','linewidth',2); grid;
title(['WFIP2 PS01 CO_2 Spectrum, ',Dstr_zz]);
ylabel('S_{cc}(f) x f'); xlabel('f (Hz)');

figure;
semilogx(fson,fson.*nanmean1(Cuw03(yy,:),1),'b-',fson,fson.*nanmean1(Cuw10(yy,:),1),'r-','linewidth',2); hold on;
semilogx([1e-4,10],[0,0],'k--'); grid;
title(['WFIP2 PS01 UW Cospectrum, ',Dstr_yy]);
ylabel('C_{uw}(f) x f'); xlabel('f (Hz)');
legend('3m','10m','location','best');

figure;
semilogx(fson,fson.*nanmean1(Ctw03(yy,:),1),'b-',fson,fson.*nanmean1(Ctw10(yy,:),1),'r-','linewidth',2); hold on;
semilogx([1e-4,10],[0,0],'k--'); grid;
title(['WFIP2 PS01 WT Cospectrum, ',Dstr_yy]);
ylabel('C_{wt}(f) x f'); xlabel('f (Hz)');
legend('3m','10m','location','best');

figure;
semilogx(fson,fson.*nanmean1(Cqw03(zz,:),1),'b-','linewidth',2); hold on;
semilogx([1e-4,10],[0,0],'k--'); grid;
title(['WFIP2 PS01 Wq Cospectrum, ',Dstr_zz]);
ylabel('C_{wq}(f) x f'); xlabel('f (Hz)');



