
%
% Loads the MVP data from the .mat files generated with L0toL1.m and
% processes them in order to align T with cond/SV due to the different 
% response times
% 
% The optimum delay is found through a cross-correlation technique and the
% results are fitted through a model depending on the vertical velocity of
% the fish: c0 + c1/vSpeed
% 
% Ullmand and Hebert, JAOT, 31, 2014
%
%
% v0
%   Lucio Bellomo, 23 June 2015
%


clear all; clc;

% Add the subfolders needed here to the Matlab path
addpath('../functions/',genpath('../GSW_Toolbox/'));


%% Parameters to be set by the user
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Group of transects to be processed
% - LDA_Cross, LDA_Before, LDA_After, LDB_Cross, LDB_Before, LDB_After, LDC_Cross
% groupOfTransects = {'LDA_Cross','LDA_Before'};
% groupOfTransects = {'LDA_After'};
% groupOfTransects = {'LDB','LDC'};
groupOfTransects = {'LDB'};

% Minimum pressure value to be reached to consider a profile as valid
minValidPres = 10;

% Low-pass filter of the data (s)
tauPres   = 2;    % pres needs a lot of filtering because its accuracy is 3 m and the vertical velocity of the MSFFF >= 3 m/s
tauOthers = 0.25;

% Minimum pressure value to consider a profile valid
presMin = 10;

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

% Folder where the ProfilesDescription.info file is
pathInfo = '/media/NeNe/DatiVari/AtalanteDisk/OUTPACE_MVP/MVP_DATA/All_Data/';

% Folder where the .mat are
pathL1 = '/media/NeNe/DatiVari/AtalanteDisk/OUTPACE_MVP/MVP_DATA/All_Data/L1/';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% Read the transect filenames
fileNames = dir([pathL1 filesep '*.mat']);
nProf = length(fileNames);


%% Read the ProfilesDescription.info file and fetch the time intervals corresponding to the considered group(s) of transect
for iGrp = 1 : length(groupOfTransects)
    fid = fopen([pathInfo filesep 'ProfilesDescription.info']);
    iTrans = 1;
    nCharGrp = length(groupOfTransects{iGrp});
    while ~feof(fid)
        ln = fgetl(fid);
        if isempty(ln)
            continue;
        end
        if strcmp(ln(1:nCharGrp),groupOfTransects{iGrp})
            ind = strfind(ln,char(9));
            group(iGrp).transectName{iTrans} = ln(1:ind(1)-1);
            group(iGrp).dateStart(iTrans,1)    = datenum(ln(ind(1)+1:ind(2)-1),'yyyy-mm-dd_HHMMSS');
            group(iGrp).dateStop(iTrans,1)     = datenum(ln(ind(2)+1:end),     'yyyy-mm-dd_HHMMSS');
            iTrans = iTrans+1;
        end
    end
    fclose(fid);
end


%% Process the profiles
vSpeedAll = [];
delayAll  = [];
coefsAll  = [];
% Group by group
for iGrp = 1 : length(groupOfTransects)
    
    
    vSpeedGrp = [];
    delayGrp  = [];
    % Transect by transect
    for iTrans = 1 : length(group(iGrp).transectName)
        
        hWait = waitbar(0,['Processing .mat files - '  strrep(group(iGrp).transectName{iTrans},'_','\_')]);   % waitbar for each transect
        
        
        % Profile by profile
        for iProf = 1 : nProf
            
            % Check the date and process the file only if it lies within
            % the interval of the dates of the current group
            fileDate = datenum(fileNames(iProf).name(5:end-4),'yyyy-mm-dd_HHMMSS');
            if fileDate >= group(iGrp).dateStart(iTrans)-1/86400 && fileDate <= group(iGrp).dateStop(iTrans)+1/86400    % to avoid finite computer resolution errors, add/subtract one second (1/86400 days)
                fileName = [pathL1 filesep fileNames(iProf).name];
                load(fileName); % this gives three structures: info, metadata and data
                
                % - check whether it is empty
                if isempty(data)
                    continue;
                end
                
                % - remove unuseful fields
                data = rmfield(data,{'LOPCcounts','analogOutput'});  % remove LOPCcounts as we do not know yet how to investigare the assessment
                
                % - low-pass filter data
                if params.fSampl.value >= 16
                    % -- 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,'SV')
                        data.soundVel.value = filter_updn(b,a,data.soundVel.value,1);
                    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
                
                % - 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');
                indEnd   = indStart + nChunk-1;
                iDel = 1;
                while indEnd <= length(data.pres.value)
                    dTemp = diff(data.temp.value(indStart:indEnd));
                    if strcmp(params.instrType,'uCTD+fluo')
                        dOther = diff(data.cond.value(indStart:indEnd));
                    else
                        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 % date inside the transect date range
            
            waitbar(iProf/nProf);
            
        end % profiles
        
        % Close waitbar
        close(hWait);
        pause(0.5);
        
    end % transects
    
    
    % 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 strcmp(params.instrType,'SVPT+fluo')
        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(iGrp); 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 ' groupOfTransects{iGrp}],'Interpreter','none');
    
    vSpeedAll = [vSpeedAll; vSpeedGrp];
    delayAll  = [delayAll;  delayGrp];
    coefsAll  = [coefsAll,  coefsGrp];
    
end % groups


%% Evaluate the optimum delay model parameters (c0 + c1/vSpeed)
% Sort values and remove outliers
[vSpeedAll, indSort] = sort(vSpeedAll);
delayAll = delayAll(indSort);

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

% Plot all the delays
figure(iGrp+1); hold on; grid on;
plot(vSpeedAll,delayAll,'b*');
plot(vSpeedAll,coefsOverall(1) + coefsOverall(2)./vSpeedAll,'r');
xlabel('Vertical speed (m/s)');
ylabel('Delay (s)');
title('ALL profiles');
