%{
    Hs_spectrogram_PISTON_MISOBOB_2018.m

    Computing Hs from 4*std(h), where std(h) is computed 3 ways:
    1) From std of the crudely corrected surface displacement
        time series (-lidar_range + heave).
    2) From std of the motion corrected surface displacement
        time series (CF equations).
    3) From the integral of the psd2 variance spectrum of the interpolated
        surface displacement timeseries (i.e. h with no NaNs)

    Standard deviations are estimated from the median & 16th/84th percentiles
        of wave height to further discriminate against spikes and outliers.

    Computes spectra, moments, and stats for every 10 min interval with
    sufficient data.

    Fp and Tp are computed when a variance peak between 0.03 and 0.3 Hz is
    detected with magnitude > 0.001.

    Hs by methods 1 & 2 are always computed.

    Feb 2018 BWB
%}

%% Initialize run parameters
close all;
fclose all;
clear;
warning ('off','MATLAB:interp1:NaNstrip');
warning ('off','MATLAB:MKDIR:DirectoryExists');

% cruise ID and paths
% data_drive = '/Users/byronb/Raw_Data_Files.noindex';
data_drive = '/Users/eliz/DATA/';
cruise = 'ATOMIC_2020';    % Cruise name
ship = 'Brown';                      % Research vessel
yyyy = '2020';
plotit = true;

% matlab script path
restoredefaultpath;
path_prog = fullfile(data_drive,cruise,ship,'scientific_analysis','programs','Riegl_Hs_fixed');
cd(path_prog);                  % cd to matlab directory
addpath(genpath(path_prog));	% add scripts to path
rehash toolboxcache

% set leg number manually, 2-5
% select dates and location based on leg number
% *** leg 1, Taiwan to Taiwan was cancelled ***
% leg 2 jd 155:171, MISOBOB1 Chennai to Colombo
% leg 3 jd 181:201, MISOBOB2 Colombo to Colombo
% leg 4 jd 219:254, PISTON2 Taiwan to Palau
% leg 5 jd 255:290, PISTON3 Palau to Taiwan

stjd = 251; endjd = 251;

% build paths relative to data_drive
% these and subsequent path constructs are OS-independent due to use of fullfile()
% these first 4 paths should be preexisting
path_raw_data_flux = fullfile(data_drive,cruise,ship,'flux','Raw');
path_proc_data_flux = fullfile(data_drive,cruise,ship,'flux','Processed');
path_proc_images_flux = fullfile(data_drive,cruise,ship,'flux','Processed_Images');
path_motionData = fullfile(path_proc_data_flux,['motion_decorr']);

path_wave_plots = fullfile(path_proc_images_flux,'riegl_wave_plots');
mkdir(path_wave_plots);
mkdir(fullfile(path_wave_plots,'hourly_raw_riegl'));
mkdir(fullfile(path_wave_plots,'hourly_spectra'));
mkdir(fullfile(path_wave_plots,'Hs'));
mkdir(fullfile(path_wave_plots,'Tp'));
path_wave_data = fullfile(path_proc_data_flux,'riegl_wave_data');
mkdir(path_wave_data);
mkdir(fullfile(path_wave_data,'Spectra'));
mkdir(fullfile(path_wave_data,'Stats'));

% various constants
theta = 5*pi/180;      % tilt angle of riegl, radians

%% begin daily loop here

for ddd = stjd:endjd
    jd = sprintf('%03i',ddd);
    jd_default = ddd + (0:10:1430)/1440; jd_default = jd_default(:); % 10 min timestamp
    WvStats = NaN(144,17); WvStats(:,1) = jd_default;
    WvSpec = NaN(144,1+3*65); WvSpec(:,1) = jd_default;
    WvSpecRaw = NaN(144,1+513*2); WvSpecRaw(:,1) = jd_default;
    
    % begin hourly loop here
    for hhh = 0:23
        hr = sprintf('%02i',hhh);
        jd_10 = ddd + hhh/24 + (0:10:50)/1440; jd_10 = jd_10(:);
        missing = NaN(6,1);
        Hs = NaN(6,1);
        Hs_cf = NaN(6,1);
        moms = NaN(6,11);
        Tp = NaN(6,1);
        Fp = NaN(6,1);
        Sws1 = NaN(6,65);
        Swraw = NaN(6,513);

        % define paths to data files
        dfl1 = fullfile(path_raw_data_flux,[yyyy(3:4),jd],['lsr0',yyyy(3:4),jd,hr,'_raw.txt']);
        dfl2 = fullfile(path_motionData,['motion_',jd,'_',hr,'_decorr.txt']);

        % if wave or motion file missing, continue to next hour
        if (exist(dfl1,'file')==0 || exist(dfl2,'file')==0)
            continue
        else
            disp(['Reading jd ',jd,' hour ',hr]);
            %% read riegl data
            [~,wave,jd_wave,med_range,raw_missing] = Riegl_reader(dfl1,ddd,hhh);
            f_riegl = 1/(mean(diff(jd_wave))*86400);    % riegl sampling frequency, Hz
            if raw_missing/length(wave)>0.1
                continue
            end
            % crudely despike with std limit
            wave(abs(wave)>4*nanstd1(wave)) = NaN;

            wave2 = interp1(jd_wave,wave,jd_wave,'pchip','extrap'); % clean array without NaNs
            
            figure; plot(jd_wave(1000:1400), wave2(1000:1400),'r-',jd_wave(1000:1400), wave(1000:1400),'b.'); grid;
            xlabel('DOY'); ylabel('Riegl range (m)'); xlim([jd_wave(1000),jd_wave(1400)]);
            title(['Riegl Wave Height: jd ',sprintf('%03i',ddd),'; hr ',sprintf('%02i',hhh)]);
            ppath = fullfile(path_wave_plots,'hourly_raw_riegl',['riegl_range_',jd,'_',hr,'.png']);
            print('-dpng',ppath);
            close all

            %% read motion
            motdata = load(dfl2);
            tmotion = ddd + hhh/24 + (0:35999)/864000; tmotion = tmotion(:);
            heave = motdata(:,7);
            roll = motdata(:,8);
            pitch = motdata(:,9);

            %% Interp to wave timebase
            heaveint = interp1(tmotion,heave,jd_wave,'linear','extrap');
            pitchint = interp1(tmotion,pitch,jd_wave,'linear','extrap');
            rollint = interp1(tmotion,roll,jd_wave,'linear','extrap');
            f_mot = f_riegl;

            %% divide into 10-min segments and process
            
            N = length(wave2);
            kk = 1;
            for jj = 1:N/6:N
                jdseg = jd_wave(jj:jj+N/6-1);
                seg = wave2(jj:jj+N/6-1);
                hve = heaveint(jj:jj+N/6-1);
                pit = pitchint(jj:jj+N/6-1);
                rol = rollint(jj:jj+N/6-1);
                M = length(seg);
                missing(kk) = length(find(isnan(seg)));    % count missing points in raw riegl
                % compute Hs from raw ht timeseries with crude heave correction
                % use 16th & 84th percentiles for best estimate of std(h)
                pct = prctilex((-seg+hve),[16,50,84]); std_seg = mean(diff(pct));
                Hs(kk) = 4*std_seg;

                % Correction from CF equation
                % eta = z-L*[sin(theta)sin(pit)+cos(theta)*SQRT(1-(sin(pit)^2-(sin(roll))^2)];
                % where theta is the mean tilt forward of the riegl
                % pit is pitch euler angle (bow down)
                % roll is roll euler angle (port up).
                % z is obtained from our motion system that is co-located with the riegl.
                z = hve + med_range;
                L = seg + med_range;
                h = z - L.*(sin(theta).*sin(pit)+cos(theta).*sqrt(1-(sin(pit).^2-(sin(rol)).^2)));
                h = h - nanmean1(h);
                % compute Hs from raw ht timeseries using 16th & 84th
                % percentiles for best extimate of std deviation
                pct = prctilex(h,[16,50,84]); std_h = mean(diff(pct));
                Hs_cf(kk) = 4*std_h;

                % High pass filter clean wave array: 0.05 Hz passband / 0.04 Hz stopband
                wp = 0.05/(f_mot/2);  % high pass filter coefficients for later use
                ws = 0.04/(f_mot/2);
                [n,wn] = buttord(wp,ws,3,7);
                [bhigh,ahigh] = butter(n,wn,'high');
                h2 = filtfilt(bhigh,ahigh,h);
                % spectrum of clean array with psd2 function
                win = hamming(1024);
                [Sw,F] = psd2(detrend(h2),1024,f_mot,win);
                Swraw(kk,:) = Sw;
                [Sws,Fsw,Df] = specsmoo4wave(Sw,f_mot);
                Sws1(kk,:) = Sws;
                % compute moments and associated stats from smoothed spectrum
                % output moms = [Hs,Tz,Te,Tpc,m0,m1,m2,m3,m4,mn1,mn2]
                moms(kk,:) = WaveStatsFromSpectra(Sws,Fsw,Df);

                % look for a peak in the spectrum between 0.04 and 0.3 Hz
                rr = find(Fsw>0.04,1,'first');
                ss = find(Fsw<0.3,1,'last');
                [maxSws,idx] = max(Sws(rr:ss));
                idx = idx + rr -1;
                % compute Fp & Tp from psd2 spectrum
                % set some limits for an acceptable peak
                if (maxSws > 0.001)
                    Fp(kk) = Fsw(idx);
                    Tp(kk) = 1/Fsw(idx);
                else
                    Fp(kk) = NaN;
                    Tp(kk) = NaN;
                end
            kk = kk+1;
            end
            
            %% Plot mean hourly wave spectrum
            Sws1m = nanmean1(Sws1,1);
            if ~(all(isnan(Sws1m)))
                ff = find(Fsw>0.01,1,'first');
                figure; loglog(Fsw(ff:end)',Sws1m(ff:end),'ro-');
                ylim([max(min(Sws1m(ff:end)),10^-6),10^(ceil(max(log10(Sws1m(ff:end)))))]);
                xlabel('Frequency (Hz)'); ylabel('S_{hh} (m^2/Hz)'); grid;
                title(['Wave Spectrum jd ',sprintf('%03i',ddd),'; hr ',sprintf('%02i',hhh)]);
                ppath = fullfile(path_wave_plots,'hourly_spectra',['wave_spectrum_',jd,'_',hr,'.png']);
                print('-dpng',ppath);
                close all
            end
        end
        
        %% save hourly results and go to next hour
        zzz = hhh*6 + (1:6);    % Indicies for this hour in daily arrays
        WvStats(zzz,:) = [jd_10,missing,Hs,Hs_cf,moms,Fp,Tp];
        if ~all(isnan(Sws1))
            Fw6 = [F';F';F';F';F';F'];
            Fsw6 = [Fsw;Fsw;Fsw;Fsw;Fsw;Fsw];
            Df6 = [Df;Df;Df;Df;Df;Df];
            WvSpec(zzz,:) = [jd_10,Fsw6,Df6,Sws1];
            WvSpecRaw(zzz,:) = [jd_10,Fw6,Swraw];
        end
    end
    
    %% save results for this day and go to next day
    WvStats_fname = ['WvStats_',jd,'.txt'];
    fl = fullfile(path_wave_data,'Stats',WvStats_fname);
    save(fl,'WvStats','-ascii');
    WvSpec_fname = ['WvSpec_',jd,'.txt'];
    fl = fullfile(path_wave_data,'Spectra',WvSpec_fname);
    save(fl,'WvSpec','-ascii');
    WvSpecRaw_fname = ['WvSpecRaw_',jd,'.mat'];
    fl = fullfile(path_wave_data,'Spectra',WvSpecRaw_fname);
    save(fl,'WvSpecRaw');
    
    % daily plots
    if plotit && any(isfinite(WvStats(:,3)))
        % 10 min Hs
        figure; plot(WvStats(:,1),WvStats(:,3),'bo',WvStats(:,1),WvStats(:,4),'ro',WvStats(:,1),WvStats(:,5),'go'); grid;
        ylabel('Hs (m)'); xlabel('DOY'); ylim([0,max(0.25,max(WvStats(:,3))+1)]);
        xlim([WvStats(1,1),WvStats(end,1)]);
        title(['Significant Wave Height: jd ',sprintf('%03i',ddd)]);
        legend('4*std(-range+heave)','4*std(hcf)','4*sqrt(IntSpec)','location','best');
        ppath = fullfile(path_wave_plots,'Hs',['Hs_riegl_',jd,'.png']);
        print('-dpng',ppath);
        
        % 10 min Tp
        figure; plot(WvStats(:,1),WvStats(:,17),'bo'); grid;
        ylabel('Tp (s)'); xlabel('DOY'); ylim([0,max(1,max(WvStats(:,17))+1)]);
        xlim([WvStats(1,1),WvStats(end,1)]);
        title(['Wave Period: jd ',sprintf('%03i',ddd)]);
        ppath = fullfile(path_wave_plots,'Tp',['Tp_riegl_',jd,'.png']);
        print('-dpng',ppath);
        close all
    end
end

% concatenate daily result files into single whole-cruise file
filecat(fullfile(path_wave_data,'Stats','WvStats_*.txt'), fullfile(path_wave_data,[cruise,'_WaveStats.txt']),'Overwrite');
filecat(fullfile(path_wave_data,'Spectra','WvSpec_*.txt'), fullfile(path_wave_data,[cruise,'_WaveSpectra.txt']),'Overwrite');
