%{
    Hs_spectrogram_CAPRICORN_2016.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)

    Raw Riegl data are processed to remove as much interference from
        bow spray as possible.

    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 where data
    loss is less than 50% (<3000 points out of 6000).

    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.

    Nothing is computed for hours where all data are bad or missing.

    SCRIPT EXECUTION IS SLOW DUE TO INTERP & NAN FILLING OPERATIONS...

    Feb 2016 BWB
%}

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

% matlab script path
path_prog = fullfile('c:\Users\Sergio\Documents\MATLAB\2016_CAPRICORN\Riegl_Hs_fixed');
cd(path_prog);                  % cd to matlab directory
restoredefaultpath              % reset matlab path
addpath(genpath(path_prog));	% add scripts to path
savepath; rehash;

% set start/end day of year
stjd = 106;
endjd = 106;

% cruise ID and paths
cruise = 'CAPRICORN_2016';    % Cruise name
ship = 'Investigator';             % Research vessel
yyyy = '2016';
plotit = true;

% cruise data root directory - edit as required for your file system and OS
% note: windows uses the '\' separator and mac-unix uses '/'
% data_drive = '/Users/byronb/Raw_Data_Files.noindex/';
data_drive = 'D:\Data';

% 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,'missing'));
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(:);
    WvStats = NaN(144,17); WvStats(:,1) = jd_default;
    WvSpec = NaN(144,1+3*78); WvSpec(:,1) = jd_default;
    WvSpecRaw = NaN(144,1+1025*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,78);
        Swraw = NaN(6,1025);

        % define paths to data files
        dfl1 = fullfile(path_raw_data_flux,[yyyy(3:4),jd],['las0',yyyy(3:4),jd,hr,'_raw.txt']);
        dfl2 = fullfile(path_motionData,['motion_v4_',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,~,~] = Riegl_reader(dfl1,ddd,hhh);
            f_riegl = 1/(mean(diff(jd_wave))*86400);    % riegl sampling frequency, Hz
            if all(isnan(wave))
                continue
            end;
            % crudely despike with std limit
             wave(abs(wave)>4*nanstd1(wave)) = NaN;
            % OR despike by rejecting points where d/dt > 7 m/s
%             zz = find(isfinite(wave));
%             wdiff = diff(wave(zz))*f_riegl; wdiff(end+1) = 0;   % d/dt in m/s
%             wdiff2 = NaN(size(wave)); wdiff2(zz) = wdiff;
%             zzz = find(abs(wdiff2)>7); wave(zzz) = NaN;
%             wave(abs(wave)>4*nanstd1(wave)) = NaN; % then reject std outliers

            wave2 = interp1(jd_wave,wave,jd_wave,'pchip','extrap'); % clean array without NaNs

            if plotit && (hhh==0 || hhh==6|| hhh==12|| hhh==18)   % save timeseries plot of raw riegl range data
                figure; plot(jd_wave, wave,'b.'); grid; xlabel('DOY'); ylabel('Riegl range (m)');
                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
            end;

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

            %% Interp to common 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');
            wave2int = interp1(jd_wave,wave2,tmotion,'linear','extrap');

            %% divide into 10-min segments and process
            seconds = ((jd_wave*24) - floor(jd_wave(1)*24))*3600;
            aa = find(seconds<600,1,'last');
            bb = find(seconds<1200,1,'last');
            cc = find(seconds<1800,1,'last');
            dd = find(seconds<2400,1,'last');
            ee = find(seconds<3000,1,'last');
            ii = [0,aa,bb,cc,dd,ee,length(wave)];
            kk = 0:6000:36000;

            for jj = 1:6
                tmot = tmotion(kk(jj)+1:kk(jj+1));
                seg2 = wave2int(kk(jj)+1:kk(jj+1));
                hve2 = heave(kk(jj)+1:kk(jj+1));
                pit2 = pitch(kk(jj)+1:kk(jj+1));
                rol2 = roll(kk(jj)+1:kk(jj+1));

                jdseg = jd_wave(ii(jj)+1:ii(jj+1));
                seg = wave(ii(jj)+1:ii(jj+1));
                hve = heaveint(ii(jj)+1:ii(jj+1));
                pit = pitchint(ii(jj)+1:ii(jj+1));
                rol = rollint(ii(jj)+1:ii(jj+1));
                N = length(seg);
                missing(jj) = 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(jj) = 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 + 14.96;
                L = seg + 14.96;
                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(jj) = 4*std_h;

                if (missing(jj)>3000 || N<5900)    % skip spectra if >50% missing
                    continue
                else
                    % corrected, clean wave height with NaNs filled in
                    z2 = hve2 + 14.96;
                    L2 = seg2 + 14.96;
                    N2 = length(seg2);
                    h2 = z2 - L2.*(sin(theta).*sin(pit2)+cos(theta).*sqrt(1-(sin(pit2).^2-(sin(rol2)).^2)));
                    h2 = h2 - nanmean1(h2);
                    % 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,h2);
                    % spectrum of clean array with psd2 function
                    win = hamming(2048);
                    [Sw,F] = psd2(detrend(h2),2048,f_mot,win);
                    Swraw(jj,:) = Sw;
                    [Sws,Fsw,Df] = specsmoo4wave(Sw,f_mot);
                    Sws1(jj,:) = Sws;
                    % compute moments and associated stats from smoothed spectrum
                    % output moms = [Hs,Tz,Te,Tpc,m0,m1,m2,m3,m4,mn1,mn2]
                    moms(jj,:) = 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(jj) = Fsw(idx);
                        Tp(jj) = 1/Fsw(idx);
                    else
                        Fp(jj) = NaN;
                        Tp(jj) = NaN;
                    end;

                    % compute spectrum from array with missing data - autocovariance method
                    % first, detrend with linear fit
%                     win = hamming(N);
%                     xx = find(isfinite(h));
%                     coef = polyfit(jdseg(xx),h(xx),1);
%                     trend = polyval(coef,jdseg);
%                     h = h - trend;
%                     h_win = h.*win;
%                     acw = nanacv(h_win,N);      % compute autocovariance w/NaNs in data - slow
%                     Sww = ifft(acw);            % inverse fft of acv
%                     f2 = 0:N;
%                     ff2 = f2'*f_riegl/2/length(f2);  % frequency array
%                     ff20 = ff2(2);                  % fundamental frequency
%                     % compute spectrum from ifft
%                     Syw = 2*abs(Sww(1:N+1))/ff20*2.52;      % 2.52 for hamming normalization
% %                    [Sws_acv,Fsw_acv] = specsmoo(Syw,f_riegl);    % smoothed spectrum
%                     [Sws_acv,Fsw_acv] = logSmoothSpec_Nbins(Syw,ff2,78);
%                     Sws0(jj,:) = Sws_acv';
                end;
            end;

            if plotit && (hhh==0 || hhh==6|| hhh==12|| hhh==18)    % save specrta plot
                Sws1m = nanmean1(Sws1,1);
%                 Sws0m = nanmean1(Sws0,1);
                if ~(all(isnan(Sws1m))) % && all(isnan(Sws0m)))
                    ff = find(Fsw>0.01,1,'first'); % gg = find(Fsw_acv>0.01,1,'first');
                    figure; loglog(Fsw(ff:end)',Sws1m(ff:end),'ro-'); %,Fsw_acv(gg:end)',Sws0m(gg:end),'bo-');
                    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; %legend('psd2','nanacv');
                    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;

        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];
%             Facv6 = [Fsw_acv';Fsw_acv';Fsw_acv';Fsw_acv';Fsw_acv';Fsw_acv'];
%             WvSpec(zzz,:) = [jd_10,Facv6,Sws0,Fsw6,Df6,Sws1];
            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
        % 10 min missing data
        figure; plot(WvStats(:,1),WvStats(:,2),'bo-'); grid;
        ylabel('missing points per 10 min'); xlabel('DOY');
        title(['Missing Data: jd ',sprintf('%03i',ddd)]);
        ppath = fullfile(path_wave_plots,'missing',['missing_riegl_',jd,'.png']);
        print('-dpng',ppath);

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