% loads & plots sp files
% 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.

% sp file format:
%{
columns     variable        description
1           doy             decimal day of year
2:49        Fson_u          frequency bins
50:97       Su              U spectrum, Edson correction
98:145      Sv              V spectrum, Edson correction
146:193     Sw              W spectrum, Edson correction
194:241     St              T spectrum
242:289     Sq              sp.hum. spectrum, LI7500
290:337     Cuw             UW cospectrum
338:385     Cvw             VW cospectrum
386:433     Ctw             TW cospectrum
434:481     Cqw             qW cospectrum
%}

%% initialize
clear;
close all;

%% load data

% insert path to desired sp file here and edit file name as required
% SA = load 'D:\Data\Calwater_2015\RHB\flux\Processed\procturb\daily\Calwater_2015_sp_v4_021.txt'
% DA = load 'D:\Data\Calwater_2015\RHB\flux\Processed\procturb\daily\Calwater_2015_da_v4_021.txt'

% or... load sp & da data for entire cruise
% SA = load ('/Users/byronb/Raw_Data_Files.noindex/CalWater_2015/RHB/flux/Processed_v4/CalWater_2015_sp_v4_decorr.txt');
SA = load('D:\Data\CAPRICORN_2016\Investigator\flux\Processed\CAPRICORN_2016_sp_decorr.txt');

% DA = load ('/Users/byronb/Raw_Data_Files.noindex/CalWater_2015/RHB/flux/Processed_v4/CalWater_2015_da_v4_decorr.txt');
DA = load('D:\Data\CAPRICORN_2016\Investigator\flux\Processed\CAPRICORN_2016_da_decorr.txt');

lic = load('D:\Data\CAPRICORN_2016\Investigator\flux\Processed\CAPRICORN_2016_cr_decorr.txt');

%% define criteria and load selected spectra
jd = DA(:,1);
u = DA(:,16);
reldir = DA(:,42);     % rel wind dir
mean_sog = DA(:,97);
sig_h = DA(:,113);     % sdev heading
sig_u = DA(:,98);      % sdev ship speed
org = DA(:,24);        % rain rate
sp2 = sqrt(DA(:,130)); % vplat sdev
missing = DA(:,158);   % missing lines of sonic data in each 10-min interval
bad = DA(:,159);       % # bad sonic pnts per 10-min
cql = DA(:,140);       % cq by method a, after fixit
agc = lic(:,17);

clear Su Sv Sw St Sq Cuw Cvw Ctw jj

% apply criteria from da file to select spectra
jj = find(jd<76 & (reldir<60 | reldir>-60) & sig_h<5 & sig_u<1.4 & missing<50 & org<5 & sp2<1.2);
% criteria for high winds, on station,
% jj = find((reldir<60 | reldir>-60) & sig_h<10 & sig_u<0.4 & sp2<1.2 & u>10 & missing<50 & org<1 & agc<60);
% criteria for low winds, underway,
% jj = find((reldir<60 | reldir>-60) & sig_h<5 & sig_u<1.4 & sp2<1.2 & u<8 & missing<50);

% read desired spectra
kk = find(isfinite(SA(:,2)),1,'first');
fson = SA(kk,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);

%% plot mean spectra
close all

figure(1);
loglog(fson,fson.*nanmean1(Su,1),'b-');
title('U spectra');
grid;

figure(2);
loglog(fson,fson.*nanmean1(Sw,1),'b-');
title('W spectra');
grid;

figure(3);
loglog(fson,fson.*nanmean1(St,1),'b-');
title('T spectra');
grid;

figure(4);
loglog(fson,fson.*nanmean1(Sq,1),'b-');
title('q spectra');
grid;

figure(5);
semilogx(fson,fson.*nanmean1(Cuw,1),'b-');
title('UW Cospectra');
grid;

figure(6);
semilogx(fson,fson.*nanmean1(Ctw,1),'b-');
title('WT Cospectra');
grid;

figure(7);
semilogx(fson,fson.*nanmean1(Cqw,1),'b-');
title('Wq Cospectra');
grid;