% clear;
warning off;

% spectral moments table
final_table = [];

% loading data for sonic wind speed
flux = load('\\shares\psdshares\psd3\cruises\CALWATER_2015\RHB\flux\Allcruise_met_flux_v3\Calwater_2015_PSD_flux_5min_allcruise3.txt');
U = flux(:,2);
uday = flux(:,1);
k = 1;
for i = 15:41
    for j = 0:23
        U_hr(k,1) = mean(U(uday >= (i + j/24) & uday <= (i + (j+1)/24)));
        U_ind(k,1:3) = [i+j/24 i j];
        k = k+1;
    end
end

% bin calculation on one file example
day = '013';
hour = '04';

fid = fopen(['\\shares\psdshares\psd3\cruises\CALWATER_2015\RHB\Riegl\raw\rgl15' day hour '_raw.txt']);
header = fgets(fid);
temp = textscan(fid, ['%2f%2f%3f' repmat('%f',1,333*4)],'Delimiter',' ,');
riegl_data = cell2mat(temp);
fclose(fid);

timeind = str2num(day) + str2num(hour)/24 + riegl_data(:,1)/60/24 + riegl_data(:,2)/60/60/24 +riegl_data(:,3)/1000/60/60/24;

angle = [timeind riegl_data(:,6:4:end)];

angle(:,2:end) = angle(:,2:end)./40*36;
angle(:,2:end) = 180 - (angle(:,2:end) + 45);

angle_edge = angle(1,2:13:end);
angle_diff = diff(angle_edge);

angle_bin = angle_edge(2:end) - abs(mean(angle_diff))/2;
dist = 16.3.*tan(angle_bin*pi/180);
rangex_bin = 16.3./cos(angle_bin.*pi/180);

fu = 9.81/2/pi./U_hr;
for i = 16:27 
    for j = 0:23
        day = sprintf('%03d',i)
        hour = sprintf('%02d',j)
        % read Riegl data
        fid = fopen(['\\shares\psdshares\psd3\cruises\CALWATER_2015\RHB\Riegl\raw\rgl15' day hour '_raw.txt']);
        if fid > 0
            header = fgets(fid);
            temp = textscan(fid, ['%2f%2f%3f' repmat('%f',1,333*4)],'Delimiter',' ,');
            riegl_data = cell2mat(temp);
            if isempty(riegl_data)
                continue
            end
            timeind = str2num(day) + str2num(hour)/24 + riegl_data(:,1)/60/24 + riegl_data(:,2)/60/60/24 + riegl_data(:,3)/1000/60/60/24;
            secind = (timeind-timeind(1))*3600*24;
            
            range = [timeind riegl_data(:,4:4:end)];
            amplitude = [timeind riegl_data(:,5:4:end)];
            angle = [timeind riegl_data(:,6:4:end)];
            syncTimer = [timeind riegl_data(:,7:4:end)];

            [~,c] = unique(timeind);
            timeind = timeind(c);
            range = range(c,:);
            angle = angle(c,:);
            range(range==0) = NaN;
            angle(:,2:end) = angle(:,2:end)./40*36;
            angle(:,2:end) = 180 - (angle(:,2:end) + 45);
            
            nn = size(range,1);
            if nn < 35000
                continue
            end
            
            % bin averaging
            idx = 2:13:334;
            range_bin = [];
            for k = 1:25
               range_bin(:,k) = nanmean(range(:,idx(k):idx(k+1)),2); 
            end
            
            
            % read motion data
            flag = true;
            try
                motdata=load(['\\shares\psdshares\psd3\cruises\CALWATER_2015\RHB\flux\Processed\motion_decorr_v2\motion_v2_' day '_' hour '_decorr.txt']);
            catch
               flag = false; 
            end
            
            if flag 
                tmotion=motdata(:,1);
                heave=motdata(:,7);
                roll=motdata(:,8);
                pitch=motdata(:,9);

                kk = sum(tmotion<=timeind(end));
                tmotion = tmotion(1:kk);
                heave = heave(1:kk);
                roll = roll(1:kk);
                pitch = pitch(1:kk);

                waveint = [];
                eta = [];

                for k = 1:25
                    jj = find(~isnan(range_bin(:,k)));
                    waveint(:,k) = interp1(timeind(jj), range_bin(jj,k), tmotion);
                end
                % motion correction
                z = heave + 16.3;
                L = waveint;
                theta = angle_bin.*pi/180;
                wave_spec = [];
                wave_frq = [];
                for k = 1:25
                    eta(:,k) = z - L(:,k).*(sin(theta(k)).*sin(pitch) + cos(theta(k)).*real(sqrt(1-sin(pitch).^2 - sin(roll).^2)));
                end
                % frequency of sampling
                freq = 10;
                if (i <= 18)
                    freq = 20;
                end
                % wave spectra calculation fo bins #21-23
                for k = 21:23
                    jj = find(~isnan(eta(:,k)));
                    wv = eta(jj,k);
                    n = length(wv);
                    [Sw,F] = psd(detrend(wv), n, freq, hamming(n),'linear');
                    Sw = Sw*2/freq;
                    [sws,fs,dfs]=spec_smooth(Sw,F,1/freq);
                    % bad hour correction
                    if k==23 & i==27 & j==6
                        sws = sws(1:end-1);
                        fs = fs(1:end-1);
                        dfs = dfs(1:end-1);
                    end
                    % wave spectra moments
                    wave_spec(:,k-20) = sws;
                    wave_frq(:,k-20) = fs; 
                    m0(k-20) = sqrt(sum(sws.*dfs));
                    m1(k-20) = sqrt(sum(fs.*sws.*dfs));
                    m2(k-20) = sqrt(sum(fs.*fs.*sws.*dfs));
                    m4(k-20) = sqrt(sum(fs.^4.*sws.*dfs));
                end
                % median over bins
                m0 = nanmedian(m0);
                m1 = nanmedian(m1);
                m2 = nanmedian(m2);
                m4 = nanmedian(m4);
                
                % Frequency of separation between swell and wind waves
                Uz = U_hr(U_ind(:,2)==i & U_ind(:,3)==j);
                Fu = fu(U_ind(:,2)==i & U_ind(:,3)==j);
                
                % only wind waves zero moment
                m0_hi = sqrt(sum(sws(wave_frq(:,3)>Fu).*dfs(wave_frq(:,3)>Fu)));
                
                % slope spectra calculation
                slope_eta = (eta(:,22)-eta(:,20))/(dist(22)-dist(20));
                jj = find(~isnan(slope_eta));
                sl = slope_eta(jj);
                n = length(sl);
                [Ssl,Fsl] = psd(detrend(sl), n, freq, hamming(n),'linear');
                Ssl = Ssl*2/freq;
                [Ssls,Fsls,dfss] = spec_smooth(Ssl,Fsl,1/freq);
                
                m0_sl = sqrt(sum(Ssls.*dfss));
                
                slope_frq = Fsls';
                slope_spec = Ssls';
                
                % wave spectra amplitude sigma
                wave_spec1 = interp1(wave_frq(:,1),wave_spec(:,1),wave_frq(:,3));
                wave_spec2 = interp1(wave_frq(:,2),wave_spec(:,2),wave_frq(:,3));
                wave_spec3 = wave_spec(:,3);
                wave_freq = wave_frq(:,3);
                
                wave_med = nanmedian([wave_spec1 wave_spec2 wave_spec3],2);
                sigma = nanmedian(std([wave_spec1 wave_spec2 wave_spec3],0,2));
                
                % Fpeak calculation
                [y,x] = max(wave_med(10:end-10).*wave_freq(10:end-10));
                 
                Fpeak = wave_freq(x+9);
                Cpeak = 9.81/2/pi/Fpeak;
                Sfp = y/Fpeak;
                
                % wave spectra at 0.4 Hz
                [~,idf] = min(abs(0.4-wave_freq));
                S04 = wave_med(idf);

                figure
                loglog(wave_frq, wave_spec.*wave_frq)
                set(gca, 'LineWidth',2, 'FontWeight','bold','FontSize',10)
                hold on
                hs(1) = loglog([Fpeak Fpeak],[1e-7 1e+3],'--r','linewidth',3);
                hs(2) = loglog([Fu Fu],[1e-7 1e+3],'--g','linewidth',3);
                ylabel({'Wave Spectrum f*S(f)'})
                title(['\fontsize{12pt}\bf Wave Spectrum - ' day ' ' hour]);
                xlabel('Frequency [Hz]');
                legend(hs, 'Fpeak','Fu')
                grid on

                saveas(gcf,['C:\Users\ekonopleva\Desktop\Riegl\Processed\New\wavespec_' day hour '.jpg'])
                close all

                figure
                loglog(slope_frq, slope_spec.*slope_frq)
                set(gca, 'LineWidth',2, 'FontWeight','bold','FontSize',10)
                ylabel({'Slope Spectrum f*S(f)'})
                title(['\fontsize{12pt}\bf Slope Spectrum - ' day ' ' hour]);
                xlabel('Frequency [Hz]');
                grid on

                saveas(gcf,['C:\Users\ekonopleva\Desktop\Riegl\Processed\New\slopespec_' day hour '.jpg'])
                close all
                
                final_table = [final_table; i+j/24 Uz Fpeak Sfp Cpeak S04 Fu m0 m1 m2 m4 sigma m0_hi m0_sl];
                
                wave_spec_table =[wave_frq(:,1) wave_spec(:,1) wave_frq(:,2) wave_spec(:,2) wave_frq(:,3) wave_spec(:,3) slope_frq slope_spec];
                save(['C:\Users\ekonopleva\Desktop\Riegl\Processed\New\TExt\wavespec_' day hour '.txt'],'wave_spec_table','-ascii','-tabs')                

            end
        end
        fclose all;
    end 
end
save(['C:\Users\ekonopleva\Desktop\Riegl\Processed\New\wave_moments.txt'],'final_table','-ascii','-tabs')