% loads & plots sp files for HiWinGS 2013
% v1 data is for std motion correction
% v2 data are for modified tilt correction scheme - correct tilt in
%   ship frame before adjusting for ship speed, then rotate azimuth.
% v3 data: decorrelate q-licor and Tsonic for motion

% sp file format:
%{
        10-min spec output arrays for spectra and cospectra: 144x817

        column designations as follows, all spectra are 48 bin

        COLUMNS     VARIABLE        DESCRIPTION
        1           doy             decimal day of year
        2:49        Fxx             frequency bins
        50:97       Puu             U spectrum
        98:145      Pvv             V spectrum
        146:193     Pww             W spectrum
        194:241     Ptt             T spectrum
        242:289     Pqq             q spectrum, LI7500
        290:337     Cuw             wu cospectrum
        338:385     Cvw             wv cospectrum
        386:433     Ctw             wt cospectrum
        434:481     Cqw             wq cospectrum
        482:529     Puu_uh          U spectrum, UH sonic
        530:577     Pvv_uh          V spectrum, UH sonic
        578:625     Pww_uh          W spectrum, UH sonic
        626:673     Ptt_uh          T spectrum, UH sonic
        674:721     Cuw_uh          wu cospectrum, UH sonic
        722:769     Cvw_uh          wv cospectrum, UH sonic
        770:817     Ctw_uh          wt cospectrum, UH sonic
%}

%% initialize
clear
close all

%% load data
% load sp & da data for entire cruise
load '/Users/byronb/Data_Working.noindex/HiWinGS_2013/Knorr/flux/Processed_v3/HiWinGS_2013_sp_v3_decorr.txt';
SA = HiWinGS_2013_sp_v3_decorr;
clear HiWinGS_2013_sp_v3_decorr

load '/Users/byronb/Data_Working.noindex/HiWinGS_2013/Knorr/flux/Processed_v3/HiWinGS_2013_da_v3_decorr.txt';
DA = HiWinGS_2013_da_v3_decorr;
clear HiWinGS_2013_da_v3_decorr

%% define criteria and load selected spectra

% edit and rerun these last two sections as desired to refine plots

jd = DA(:,1);
u = DA(:,38);        % streamwise relative wind speed
reldir = DA(:,42);   % relative wind direction
mean_sog = DA(:,97);
sig_h = DA(:,113);   % stdev heading
sig_u = DA(:,98);    % stdev ship speed
org = DA(:,24);      % psd org rain rate
sp2 = DA(:,130);     % v-motion std
bad = DA(:,159);     % # bad sonic pnts per 10-min interval

clear Su Sv Sw St Sq Cuw Cvw Ctw jj

% apply criteria from da file to select spectra...
%    all data
% jj = find((reldir<90 & reldir>-90) & sig_h<5 & sig_u<1.4 & bad<300 & org<5 & sp2<1.2);
%    for Davis Strait Leg
% jj = find((reldir<90 & reldir>-90) & sig_h<5 & sig_u<1.4 & sp2<1.2 & bad<300 & jd<277 & jd>257);
%    for HiWinGS Leg
jj = find((reldir<90 & reldir>-90) & sig_h<5 & sig_u<1.4 & sp2<1.2 & bad<300 & jd<317 & jd>283);

% read desired spectra
fson = SA(1,2:49);
date = SA(jj,1);

Su = SA(jj,50:97);
Sv = SA(jj,98:145);
Sw = SA(jj,146:193);
St = SA(jj,194:241);
Sq = SA(jj,242:289);
Cuw = SA(jj,290:337);
Cvw = SA(jj,338:385);
Ctw = SA(jj,386:433);
Cqw = SA(jj,434:481);
Su_uh = SA(jj,482:529);
Sv_uh = SA(jj,530:577);
Sw_uh = SA(jj,578:625);
St_uh = SA(jj,626:673);
Cuw_uh = SA(jj,674:721);
Cvw_uh = SA(jj,722:769);
Ctw_uh = SA(jj,770:817);

%% plot mean spectra
figure(1);
loglog(fson,fson.*nanmean1(Su,1),'r-',fson,fson.*nanmean1(Su_uh,1),'b-');
title('U spectra: psd:red, uh:blue');
grid;

figure(2);
loglog(fson,fson.*nanmean1(Sw,1),'r-',fson,fson.*nanmean1(Sw_uh,1),'b-');
title('W spectra: psd:red, uh:blue');
grid;

figure(3);
semilogx(fson,fson.*nanmean1(Cuw,1),'r-',fson,fson.*nanmean1(Cuw_uh,1),'b-');
title('UW Cospectra: psd:red, uh:blue');
grid;

figure(4);
semilogx(fson,fson.*nanmean1(Ctw,1),'r-',fson,fson.*nanmean1(Ctw_uh,1),'b-');
title('WT Cospectra: psd:red, uh:blue');
grid;

figure(5);
semilogx(fson,fson.*nanmean1(Cqw,1),'r-');
title('Wq Cospectra: psd:red');
grid;