% 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
% load 'D:\Data\Calwater_2015\RHB\flux\Processed\procturb\daily\Calwater_2015_sp_v1_021.txt'
% load 'D:\Data\Calwater_2015\RHB\flux\Processed\procturb\daily\Calwater_2015_da_v1_021.txt'
% SA = Calwater_2015_sp_v1_021;
% clear Calwater_2015_sp_v1_021;
% DA = Calwater_2015_da_v1_021;
% clear Calwater_2015_da_v1_021;

% or... load sp & da data for entire cruise
load '/Users/byronb/Raw_Data_Files.noindex/CalWater_2015/RHB/flux/Processed/CalWater_2015_sp_v2_decorr.txt';
% load 'D:\Data\Calwater_2015\RHB\flux\Processed\procturb\total\CalWater_2015_sp_v2_decorr.txt'
SA = CalWater_2015_sp_v2_decorr;
clear CalWater_2015_sp_v2_decorr

load '/Users/byronb/Raw_Data_Files.noindex/CalWater_2015/RHB/flux/Processed/CalWater_2015_da_v3_decorr.txt';
% load 'D:\Data\Calwater_2015\RHB\flux\Processed\procturb\total\CalWater_2015_da_v3_decorr.txt'
DA = CalWater_2015_da_v3_decorr;
clear CalWater_2015_da_v3_decorr


%% 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(:,10);        % rain rate
sp2 = sqrt(DA(:,130)); % vplat sdev
bad = DA(:,159);       % # bad sonic pnts per 10-min

clear Su Sv Sw St Sq Cuw Cvw Ctw jj

% apply criteria from da file to select spectra
jj = find(jd>25 & jd<27 & (reldir<90 | reldir>-90) & sig_h<5 & sig_u<1.4 & bad<300 & org<5 & sp2<1.2);
% criteria for high winds, on station,
% jj = find((reldir<90 | reldir>270) & sig_h<5 & sig_u<1.4 & sp2<1.2 & u>10 & bad<300 & jd<25.58333 & jd>19.50);
% criteria for low winds, underway,
% jj = find((reldir<90 | reldir>270) & sig_h<5 & sig_u<1.4 & sp2<1.2 & u<8 & bad<300 & jd<19.25 & jd>16.0);

% 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);

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

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

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

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