%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%                    LATEXtools                         %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  This  file  is  free software ;  it  is distributed  in  %
%%  the hope that  it  will  be  useful,  but  without  any  %
%%  warranty.  You  can  redistribute  it and/or modify  it  %
%%  under  the  terms  of  the  GNU  General Public License  %
%%  as  published  by  the   Free  Software  Foundation  at  %
%%       http://www.gnu.org/copyleft/gpl.html                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% MVP_lags_lucio (MATLAB FUNCTION)
%%
%% Calculates optimal lags for L1 MVP data, and plots it.
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 21/04/16 : louise.rousselet@mio.osupytheas.fr :
%%            andrea.doglioli@mio.osupytheas.fr  : 
%%            alain.de-verneil@mio.osupytheas.fr : creation 
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Created on Matlab R2014b
% v0.0 Alain de Verneil, borrowing code from align_CTD_SVP.m written by
% Lucio Bellomo
% 6 Apr 2016

function coefsGrp = MVP_lags_lucio(fileNames,Trans_ind)

MVP_cruise_params;

% Correlation computation parameters
timeChunk    = 10;  % correlation window length (s)
timeDelayMax = 2;   % maximum allowable delay between T and SV/cond (s)

vSpeedAll = [];
delayAll  = [];
coefsAll  = [];

% Group by group
vSpeedGrp = [];
delayGrp  = [];

% Profile by profile

profCnt = length(Trans_ind); % Profile Counter

for iProf = 1:profCnt
    
    fileName = [pathL1 filesep fileNames(Trans_ind(iProf)).name];
    load(fileName); % this gives three structures: info, metadata and data
    
    if params.QCcast < 1 | params.QCcast==9   % Bad cast... empty or ID'ed by user as bad..
        continue;
    end
        
    % Smooth the data, if possible...
    if fSample >= 16 %AD 9999 instead of 16
        if exist('butter')==2 % If the user can use the butterworth low-pass filter...
            % -- pressure
            fCutPres = (1/tauPres)/(params.fSampl.value/2);
            [b,a] = butter(3,fCutPres,'low');
            data.pres.value = filter_updn(b,a,data.pres.value,1);
            % -- other variables
            fCutOthers = (1/tauOthers)/(params.fSampl.value/2);
            [b,a] = butter(3,fCutOthers,'low');
            data.temp.value = filter_updn(b,a,data.temp.value,1);
            if isfield(data,'cond')
                data.cond.value = filter_updn(b,a,data.cond.value,1);
            end
            if isfield(data,'soundVel')
                data.soundVel.value = filter_updn(b,a,data.soundVel.value,1);
            end
        end
    end
    
    % - compute the vertical speed (central finite differences)
    timeGrid = metadata.date.value + [0:length(data.temp.value)-1]'/params.fSampl.value/24/60/60;  % uniform time grid (days)
    data.vertSpeed.value = [NaN; diff(data.pres.value(2:end))+diff(data.pres.value(1:end-1)); NaN]*params.fSampl.value/2;
    data.vertSpeed.unit  = 'm/s';
    
    % - separate down- and up-casts and keep only the downcast
    [~, ind] = max(data.pres.value); % maximum pressure value
    fields = fieldnames(data);
    for iField = 1 : length(fields)
        data.(fields{iField}).value = data.(fields{iField}).value(1:ind);
        %data.(fields{iField}).value = data.(fields{iField}).value(ind+1:end);
    end
    
    presMin = 10;
    % - compute correlation between temperature and SV gradients
    nChunk = timeChunk*params.fSampl.value;
    nDel   = timeDelayMax*params.fSampl.value;
    indStart = find(data.pres.value >= presMin,1,'first');
    indStart2 = find(data.vertSpeed.value> 1,1,'first');
    indStart = max([indStart indStart2]);
    indEnd   = indStart + nChunk-1;
    iDel = 1;
    while indEnd <= length(data.pres.value)
        dTemp = diff(data.temp.value(indStart:indEnd));
        if isfield(data,'cond')
            dOther = diff(data.cond.value(indStart:indEnd));
        elseif isfield(data,'soundVel')
            dOther = diff(data.soundVel.value(indStart:indEnd));
        end
        [res, lags] = xcorr(dTemp,dOther,'unbiased');
        res  = res(nChunk-nDel:nChunk+nDel);
        lags = lags(nChunk-nDel:nChunk+nDel);
        [~, ind] = max(res);
        if ind == 1 || ind == length(res)
            delayChunk = NaN;
        else
            coefs = polyfit(lags(ind-1:ind+1),res(ind-1:ind+1)',2); % 2nd degree fitting
            lagsOpt = -coefs(2)/(2*coefs(1));
            delayChunk = lagsOpt/params.fSampl.value;
        end
        vertSpeedChunk = median(data.vertSpeed.value(indStart:indEnd));
        indStart = indEnd+1;
        indEnd   = indStart + nChunk-1;
        iDel = iDel+1;
        
        vSpeedGrp = [vSpeedGrp; vertSpeedChunk];
        delayGrp  = [delayGrp;  delayChunk];
    end
    
end % profiles

pause(0.5);


% Remove NaN values
indOK = ~isnan(vSpeedGrp) & ~isnan(delayGrp);
vSpeedGrp = vSpeedGrp(indOK);
delayGrp  = delayGrp(indOK);


% Sort values and remove outliers
[vSpeedGrp, indSort] = sort(vSpeedGrp);
delayGrp = delayGrp(indSort);

if fSample < 16
    indOK = abs(delayGrp) < 1 & abs(vSpeedGrp) > 2;
else
    indOK = delayGrp > 0 & abs(vSpeedGrp) > 2;
end

% These are the original outpace conditions... I'm guessing if Hz at 1,
% like w/ the uCTD for OUTPACE, then use the abs(delayGrp) condition...

% if strcmp(params.instrType,'SVPT')
%     indOK = delayGrp > 0 & abs(vSpeedGrp) > 2;
% else
%     indOK = abs(delayGrp) < 1 & abs(vSpeedGrp) > 2;
% end

vSpeedGrp = vSpeedGrp(indOK);
delayGrp  = delayGrp(indOK);

% Evaluate the optimum delay model parameters (c0 + c1/vSpeed)
A = [length(vSpeedGrp) sum(1./vSpeedGrp); ...
    sum(1./vSpeedGrp) sum(1./vSpeedGrp.^2)];
b = [sum(delayGrp); ...
    sum(delayGrp./vSpeedGrp)];
coefsGrp = A\b;

% Plot the delays per group
figure(); hold on; grid on;
plot(vSpeedGrp,delayGrp,'b*');
plot(vSpeedGrp,coefsGrp(1) + coefsGrp(2)./vSpeedGrp,'r');
xlabel('Vertical speed (m/s)');
ylabel('Delay (s)');
title(['Delays for ' ],'Interpreter','none');

