% 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

    da 10-min variable designations as follows:  144 x 222 array

    Note:  historical da format is retained, but Hare corrections not computed
    in this version of the flux code & therefore set to NaN in the da output.

    Also note: COARE input format is for vectorized version 3.5.

    No odec on HiWinGS, so these also set to NaN.

    COL     VARIABLE        DESCRIPTION
    1       doy             decimal day of year
    2:10    xship           COARE 3.5 input, scs:[Utrue,Tair,RH,P_mb,Tsea,
                            Rs,Rl,lat,rain]
    11:15   no_data         To maintain historical column numbering...
    16:24   xpsd            COARE 3.5 input, PSD, as above
    25:29   no_data         To maintain historical column numbering...
    30      u_ship          ship mean true wind speed (m/s)
    31      u_ship_min      ship minimum true wind speed (m/s)
    32      u_ship_max      ship maximum true wind speed (m/s)
    33      u_ship_std      ship true wind speed stdev (m/s)
    34      wdir_ship       ship mean true wind direction
    35      wdir_ship_min   ship minimum true wind direction
    36      wdir_ship_max   ship maximum true wind direction
    37      wdir_ship_std   ship true wind direction stdev
    38      rwspd_psd       streamwise relative wind speed (m/s)from PSD sonic
    39      Un              N component, corrected true windspeed, PSD sonic
    40      Ue              E component, corrected true windspeed, PSD sonic
    41      wdir_psd        true wind direction from N,E components, PSD sonic
    42      rwdir_psd       relative wind direction from PSD sonic
    43      azm             azimuth angle from first rotation, radians
    44      tilt            tilt angle from first rotation, radians
    45      Tsea_psd        seasnake SST, degC
    46      Tair_psd        air temperature, PSD vaisala, degC
    47      RH_psd          RH, PSD vaisala
    48      RH_ship         RH, ship sensor
    49      Q_psd           specific humidity, PSD sensors, g/kg
    50      Q_ship          specific humidity, ship sensors, g/kg
    51      Qsea_psd        sea surface specific humidity from seasnake SST
    52      rhoa_psd        moist air density from PSD sensors (kg/m3)
    53:73   yship           COARE 3.5 output, ship data:[usr,tau,hsb,hlb,hlwebb,
                         tsr,qsr,zot,zoq,Cd,Ch,Ce,L,zet,dter,dqer,tkt,RF,
                         Cdn_10,Chn_10,Cen_10]
    74:94   ypsd            COARE 3.5 output, PSD data, as above
    95      Rl_std_psd      std PIR flux, PSD
    96      Rs_std_psd      std PSP flux, PSD
    97      sog_ship        mean sog, ship gps, m/s
    98      sog_ship_std    std sog, ship gps
    99      sog_ship_min    minimum sog, ship gps
    100     sog_ship_max    maximum sog, ship gps
    101     cog_ship        mean cog, deg, ship gps
    102     cog_ship_std    std cog, ship gps
    103     cog_ship_min    minimum cog, ship gps
    104     cog_ship_max    maximum cog, ship gps
    105     cog_psd         mean cog, deg, PSD gps
    106     cog_psd_std     std cog, PSD gps
    107     cog_psd_min     minimum cog, PSD gps
    108     cog_psd_max     maximum cog, PSD gps
    109     sog_psd         mean sog, m/s, PSD gps
    110     sog_psd_min     minimum sog, PSD gps
    111     sog_psd_max     maximum sog, PSD gps
    112     lrg             mean heading, deg, ship gyro
    113     lrg_std         std heading, ship gyro
    114     lrg_min         minimum heading, ship gyro
    115     lrg_max         maximum heading, ship gyro
    116     no_data         wu, Hare correction, not computed
    117     wu              wu covariance, Edson correction
    118     no_data         wv, Hare correction, not computed
    119     wv              wv covariance, Edson correction
    120     no_data         wT, Hare correction, not computed
    121     wt              wT covariance, Edson correction
    122     no_data         uvar, Hare correction, not computed
    123     no_data         vvar, Hare correction, not computed
    124     no_data         wvar, Hare correction, not computed
    125     uvar            u variance, Edson correction
    126     vvar            v variance, Edson correction
    127     wvar            w variance, Edson correction
    128     tvar            T variance
    129     Sp1             U platform velocity std, Edson correction
    130     Sp2             V platform velocity std, Edson correction
    131     Sp3             W platform velocity std, Edson correction
    132     hd_sin_std      stdev sine lrg heading
    133     hd_cos_std      stdev cosine lrg heading
    134     Cu2a            Cu structure func parameter, method a
    135     Cu2b            Cu structure func parameter, method b
    136     Cw2a            Cw structure func parameter, method a
    137     Cw2b            Cw structure func parameter, method b
    138     Ct2a            Ct structure func parameter, method a
    139     Ct2b            Ct structure func parameter, method b
    140     Cq2a            Ch2o structure func parameter, method a
    141     Cq2b            Ch2o structure func parameter, method b
    142     usid_a          ustar from idiss, method a
    143     usid_b          ustar from idiss, method b
    144     tsid_a          tstar from idiss, method a
    145     tsid_b          tstar from idiss, method b
    146     qsid_a          qstar from idiss, method a, kg/kg
    147     qsid_b          qstar from idiss, method b, kg/kg
    148     Lid_a           L from idiss, method a
    149     Lid_b           L from idiss, method b
    150     lon             longitude, scs gps
    151     lat             latitiude, scs gps
    152     Ubar            streamwise true windspeed from Edson correction
    153     odec            odec (ship speed through water), m/s
    154     odec_std        stdev odec
    155     odec_min        minimum odec
    156     odec_max        maximum odec
    157     P_mb            atmospheric pressure, PSD sensor, mb
    158     missing         missing data lines in son file per 10 min
    159     badPnts         # bad data points (spikes, etc.) in sonic data per 10 min
    160     hs_cov          sensible heat flux from covariance
    161     hl_cov          latent heat flux from covariance, uncorrected
    162     hlc_cov         webb corrected hl_cov
    163     sal_scs         salinity, ship system, ppt
    164:184 yuh             COARE 3.5 output, UH data, as above
    185     wu_uh           wu covariance, UH sonic
    186     wv_uh           wv covariance, UH sonic
    187     wt_uh           wt covariance, UH sonic
    188     uvar_uh         u variance, UH sonic
    189     vvar_uh         v variance, UH sonic
    190     wvar_uh         w variance, UH sonic
    191     tvar_uh         T variance, UH sonic
    192     Sp1_uh          U platform velocity std, UH motpak
    193     Sp2_uh          V platform velocity std, UH motpak
    194     Sp3_uh          W platform velocity std, UH motpak
    195     Cu2a_uh         Cu structure func parameter, method a, UH sonic
    196     Cu2b_uh         Cu structure func parameter, method b, UH sonic
    197     Cw2a_uh         Cw structure func parameter, method a, UH sonic
    198     Cw2b_uh         Cw structure func parameter, method b, UH sonic
    199     Ct2a_uh         Ct structure func parameter, method a, UH sonic
    200     Ct2b_uh         Ct structure func parameter, method b, UH sonic
    201     usid_a_uh       ustar from idiss, method a, UH sonic
    202     usid_b_uh       ustar from idiss, method b, UH sonic
    203     tsid_a_uh       tstar from idiss, method a, UH sonic
    204     tsid_b_uh       tstar from idiss, method b, UH sonic
    205     Lid_a_uh        L from idiss, method a, UH sonic
    206     Lid_b_uh        L from idiss, method b, UH sonic
    207     rwspd_uh        streamwise relative wind speed (m/s), UH sonic
    208     rwdir_uh        relative wind direction, UH sonic
    209     Un_uh           N component, corrected true windspeed, UH sonic
    210     Ue_uh           E component, corrected true windspeed, UH sonic
    211     wspd_uh         true wind speed from N,E components, UH sonic
    212     wdir_uh         true wind direction from N,E components, UH sonic
    213     azm_uh          azimuth angle from first rotation, radians, UH sonic
    214     tilt_uh         tilt angle from first rotation, radians, UH sonic
    215     hs_cov_uh       covariance sensible heat flux, UH sonic
    216     rs_psd2         psp2 solar radiation
    217     rl_psd2         pir1 IR radiation
    218     Ts_noise        psd WindMasterPro Tsonic noise variance
    219     Ts_noise_uh     uh R2A Tsonic noise variance
    220     missingLic      missing licor points per 10min interval
    221     badLic          bad licor H2O points per 10min interval
    222     missingMot      missing motion points per 10min interval
%}

%% 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
missing = DA(:,158); % missing sonic points per 10min interval
bad = DA(:,159);     % # bad sonic pnts per 10-min interval
missingLic = DA(:,220); % missing licor points per 10min interval
badLic = DA(:,221);     % bad licor H2O points per 10min interval
missingMot = DA(:,222); % missing motion points per 10min interval

clear Su Sv Sw St Sq Cuw Cvw Ctw jj

% apply criteria from da file to select spectra...
%    all data
% jj = find((reldir<60 & reldir>-60) & sig_h<5 & sig_u<1.4 & missing<50 & org<5 & sp2<1.2);
%    for Davis Strait Leg
% jj = find((reldir<60 & reldir>-60) & sig_h<5 & sig_u<1.4 & sp2<1.2 & missing<50 & jd<277 & jd>257);
%    for HiWinGS Leg
jj = find((jd<317 & jd>283 & reldir<60 & reldir>-60) & sig_h<5 & sig_u<1.4 & sp2<1.2); % for uh
kk = intersect(jj, find(missing<50)); % for psd
mm = intersect(kk,find(missingLic<50)); % for psd licor

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

Su = SA(kk,50:97);
Sv = SA(kk,98:145);
Sw = SA(kk,146:193);
St = SA(kk,194:241);
Sq = SA(mm,242:289);
Cuw = SA(kk,290:337);
Cvw = SA(kk,338:385);
Ctw = SA(kk,386:433);
Cqw = SA(mm,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.*nanmedian1(Su,1),'r-',fson,fson.*nanmedian1(Su_uh,1),'b-');
title('U spectra: psd:red, uh:blue');
grid;

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

figure(7);
loglog(fson,fson.*nanmedian1(St,1),'r-',fson,fson.*nanmedian1(St_uh,1),'b-');
title('Ts spectra: psd:red, uh:blue');
grid;

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

figure(4);
semilogx(fson,fson.*nanmedian1(Cvw,1),'r-',fson,fson.*nanmedian1(Cvw_uh,1),'b-');
title('VW Cospectra: psd:red, uh:blue');
grid;

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

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