% Function find_lags
%
% This function takes in one MVP profile, and searches iteratively for the
% best constant lag to apply that minimizes the RMS of both salinity and
% density. 
% The function does this by interpolating the raw data onto successively
% higher resolution lag grids, and searching for the smooth (relatively)
% minimum in the calculated RMS.
%
% Note: while this will work on an entire cast, preferably the
% upcast has already been removed prior to calling this function.
%
% Inputs:
%  - press:       The pressure timeseries [dBar]
%  - cond:        Conductivity timeseries [mS/cm]
%  - temp:        Temperature timeseries  [deg C]
%  - lon:         Longitude (used in GSW) [decimal deg. [0:360] OR
%                 [-180:180]
%  - lat:         Latitude (used in GSW toolbox) [XXX.XX decimal degrees N
%                 (-90:90)]
%  - interp_hi:   Upper threshold of lags to be used.
%  - interp_lo:   Lower threshold of lags to be used.
%  - delt_levels: The gap between tested lags (in Scans)
%  - Cond_flag:   Flag that tells us if this is conductivity or sound
%                 speed.. 0 for SVP, 1 for Cond
%  - dens_flag:   User says whether they want to calculate the density RMS.
%
% Outputs:
%  - best_lag:    The empirically found lag that minimizes salinity(dens).
%                 If the dens_flag is set to True, then a vector of the two
%                 values is returned.
%
% Example of use:
% Given pressure, conductivity, and temperature vectors, along with the
% latitude, the following will begin by looking at 0 lag to 10 scan lags,
% with a first loop that has a resolution of 1 scan, and a second loop that
% uses a resolution of 0.1 interpolate scans. Density is calculated also.
%
% lags = find_lags(P,C,T,lon,lat,10,0,[1 0.1],1);
%
% Created on Matlab R2014b
% v0.0 Alain de Verneil
% 6 Apr 2016

function best_lag = MVP_find_lags(press,cond,temp,lon,lat,interp_hi,interp_lo,delt_levels,Cond_flag,dens_flag)

nLevs = length(delt_levels);    % User-prescribed number of loops to run.
index_min = [];                 % variable to find minimum RMS
best_lag = [];                   % Initializes the returned variable.  
len_prof = length(press);       % Length of dataset

for iLev = 1:nLevs  % Cycle through levels of analysis.
        
    delt = delt_levels(iLev); % Sets the lag spacing.
    
    % Interpolate the raw data onto the lag-spaced grid.
    interp_space = 0:delt:len_prof;
    intCond = pchip(1:len_prof,cond,interp_space);
    intTemp = pchip(1:len_prof,temp,interp_space);
    intPres = pchip(1:len_prof,press,interp_space);
    
    delt_space = 0:delt:len_prof;
    
    
    find_new_hi = find(delt_space>=interp_hi,1,'first');
    find_new_lo = find(delt_space<=interp_lo,1,'last'); % so... this resets the range of lags based on the interp_space?
    interp_hi = delt_space(find_new_hi);
    interp_lo = delt_space(find_new_lo);
    
    rms_std_dev = zeros(1,length(interp_lo:delt:interp_hi)); % rms of salinity
    if dens_flag
        rms_std_dev2 = zeros(1,length(interp_lo:delt:interp_hi)); % rms of density
    end
    lag_vec = interp_lo:delt:interp_hi;
    
    % The Heavy Lifting!
    for qw = 1:length(lag_vec) % cycle through lag vector
        
        % Perform the lags...
        newcond2 = intCond(1:end-floor((lag_vec(qw)/delt)));
        newpress2 = intPres(1:end-floor((lag_vec(qw)/delt)));
        newtemp2 = intTemp(1+floor((lag_vec(qw)/delt)):end);
        
        if Cond_flag
            % Calculate derived variables (salinity, density) using
            % conductivity.
            newsal2 = gsw_SA_from_SP(gsw_SP_from_C(newcond2,newtemp2,newpress2),newpress2,lon,lat);
        else
            % Calculate derived vars w/ soundspeed.
            newsal2 = salinity_from_soundspeed(newcond2,newtemp2,newpress2,lat,1);
            newsal2 = gsw_SA_from_SP(newsal2,newpress2,lon,lat);
        end
        if dens_flag
            newCT2  = gsw_CT_from_t(newsal2,newtemp2,newpress2);
            newdens2 = gsw_rho_CT(newsal2,newCT2,newpress2);
        end
        
        chunk_min = []; % reset the min/max variables... this is NOT index_min used later...
        chunk_max = [];
        if length(newsal2)<200
            dchun = 10;
        else
            dchun = 100;
        end
        for m = 1:length(newsal2) % goes through every datapoint
            chunk_min = ceil(max(1,m-length(newsal2)/dchun)); % Finds a nearby chunk (1/100th) of the profile
            chunk_max = floor(min(length(newsal2),m+length(newsal2)/dchun)); % to calculate RMS from...
            
            % Add up values for RMS
            rms_std_dev(qw) = rms_std_dev(qw) + (newsal2(m)-(mean(newsal2(chunk_min:chunk_max))))^2;
            if dens_flag
                rms_std_dev2(qw) = rms_std_dev2(qw) + (newdens2(m)-(mean(newdens2(chunk_min:chunk_max))))^2;
            end
            
        end
        
        % Calculate the RMS
        rms_std_dev(qw) = rms_std_dev(qw)/length(newsal2);
        rms_std_dev(qw) = sqrt(rms_std_dev(qw));
        if dens_flag
            rms_std_dev2(qw) = rms_std_dev2(qw)/length(newdens2);
            rms_std_dev2(qw) = sqrt(rms_std_dev2(qw));
        end
    end
    
    index_min = find(rms_std_dev == min(rms_std_dev)); % THIS index_min is min of rms...
    if dens_flag
        index_min2 = find(rms_std_dev2 == min(rms_std_dev2));
    end
    
    if iLev == 1
        interp_hi = lag_vec(index_min(1))+1.5*delt;
        interp_lo = lag_vec(index_min(1))-1.5*delt;
        if interp_lo<0
            interp_lo = 0;
        end
    else
        best_lag = nanmean(lag_vec(index_min));
        if dens_flag
            best_lag = [best_lag nanmean(lag_vec(index_min2))];
        end
    end
    
end

if isempty(best_lag)
    disp('Strangely, no index_min found! Setting return value to zero.');
    best_lag = 0;
end