%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%                    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_L1toL2 (MATLAB SCRIPT)
%%
%% 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 
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all; close all;
% Load relevant cruise parameters.
MVP_cruise_params;

fileNames = dir([pathL1 filesep '*.mat']);
nFiles = length(fileNames);

transectName = [];
% Set up some metadata parameters.
if length(transchoice) == 1
    transectName = Transect_Names(transchoice);
else
    for iName = 1:length(transchoice)
        transectName = [transectName Transect_Names(transchoice(iName))];
    end
end
dateStart    = datenum(Transect_Start{transchoice(1)},'yyyy-mm-dd_HHMMSS');
dateStop     = datenum(Transect_Stop{transchoice(end)} ,'yyyy-mm-dd_HHMMSS');

%savename = [pathL2 filesep sprintf('Transect_%02d_qc_%d.mat',transchoice,qc_user)];

% Since each cast may have multiple QC states, we need to get the unique
% list of cast names, and then the index of each QC status. For the L2
% treatement, we will extract the highest QC flag available for
% processing...

% Get the indices of all the profiles selected...
d1 = dateStart; d2 = dateStop;
minInd = nFiles; maxInd = 0;

CastList = [];
Struc_List = {};

for iFile = 1:nFiles % Cycle through all files in L1 folder...
    tempname = fileNames(iFile).name;
    tempname = tempname(5:end-6); % Ignore QC flag, just get the name...
    CastList = [CastList; tempname];
    try
        tempdate = datenum(tempname,'yyyy-mm-dd_HHMMSS');
    catch
        fprintf('There seems to be a bad filename here: %s\n',fileNames(iFile).name)        
        disp('Skipping cast');
        continue
    end
    if tempdate >= (d1-60/86400) && tempdate <= (d2+60/86400)
        %CastList = [CastList; tempname];
        Struc_List = [Struc_List {tempname}];
    end
end

nCast = size(CastList); nCast = nCast(1); % Number of files belonging to casts, but some are copies prob...
UniqueCastList = unique(Struc_List); % List of unique cast files...
nUnique = length(UniqueCastList);

Index_Flag0 = NaN(nUnique,1);
Index_Flag1 = NaN(nUnique,1);
Index_Flag2 = NaN(nUnique,1);
Index_Flag3 = NaN(nUnique,1);
Index_Bad   = NaN(nUnique,1);
Index_Good  = NaN(nUnique,1);

% Cycle through each unique cast name, find where it exists, and then
% extract the QC flag index that matches with qc_user...
for iCast = 1:nUnique
    castname = UniqueCastList{iCast};
    for iFile = 1:nFiles
        tempcompare = CastList(iFile,:);
        if strcmp(castname,tempcompare)
            qc_extract = str2double(fileNames(iFile).name(end-4));
            if qc_extract == qc_flag
                Index_Good(iCast) = iFile;
            end
            
            switch qc_extract
                case 0
                    Index_Flag0(iCast) = iFile;
                case 1
                    Index_Flag1(iCast) = iFile;
                case 2
                    Index_Flag2(iCast) = iFile;
                case 3
                    Index_Flag3(iCast) = iFile;
                case 9
                    Index_Bad(iCast) = iFile;
            end
        end
    end    
end

Unique_Good = find(~isnan(Index_Good));
nGood = length(Unique_Good);


%%%%%%%%%%%%% NOW TIME TO RUN THE LAG QC %%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch qc_user
    
    %%%%%%%% RUN MATCHING OF T/C GRADIENTS METHOD %%%%%%%%%%%%%%
    case 1
        % Run Lucio's method
        
        Lag_coefs = []; % Find the best lag using lucio's method.
        %%%%%%%%%%%%% VERY IMPORTANT LINE %%%%%%%%%%%%%%%%%%%%%%%%%
        delayT = MVP_lags_lucio(fileNames,Index_Good(Unique_Good));
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        disp('Lag coefficients found!');
        
        
        disp('Now applying correction to all the chosen casts...');
        drawnow
        
        iProfOK = 0; % initialize for good casts...
        
        hWait = waitbar(0,'Processing L1 casts');
        
        for iProf = 1:nGood %Trans_All(1):Trans_All(2)
            load([pathL1 filesep fileNames(Index_Good(Unique_Good(iProf))).name]);
            if params.QCcast >=1 && params.QCcast ~= 9
                
                % For all the good casts...
                % First remove LOPC, if it exists...
                if isfield(data,{'LOPCcounts'})
                    data = rmfield(data,{'LOPCcounts'});
                end
                
                % Check to see if ADCcounts is present, and if so, convert
                % it to Voltage (Conversion parameter specified in
                % MVP_cruise_params.m)
                if isfield(data,{'ADCcounts'})
                    data.analog.value = (data.ADCcounts.value).*volt_input / 2^(bit_input);
                    data.analog.unit = 'V';
                    data = rmfield(data,{'ADCcounts'});
                end

                iProfOK = iProfOK+1;                
                
                
                % - low-pass filter data, if possible
                if filter_flag 
                    if exist('butter')==2
                        % -- pressure
                        fCutPres = (1/tauPres)/(params.fSampl.value/2);
                        [b,a] = butter(3,fCutPres,'low');
                        data.pres.value = MVP_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 = MVP_filter_updn(b,a,data.temp.value,1);
                        if isfield(data,{'cond'})
                            data.cond.value = MVP_filter_updn(b,a,data.cond.value,1);
                        end
                        if isfield(data,{'soundVel'})
                            data.soundVel.value = MVP_filter_updn(b,a,data.soundVel.value,1);
                        end
                    else
                        disp('The function butter does not exist...');
                        pause(2);
                    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);
                    % Below keeps only upcast...
                    %data.(fields{iField}).value = data.(fields{iField}).value(ind+1:end);
                end
                
                % Delay the Temperature
                if ~all(isnan(delayT))
                    timeDelayT                                = delayT(1) + delayT(2)./data.vertSpeed.value;
                    timeDelayT(abs(data.vertSpeed.value) < 1) = delayT(1) + delayT(2)/1;        % if vertSpeed is lower than 1 m/s, the second term of (1) starts to explode. Hence vertSpeed is set to 1.
                    timeDelayT(isnan(timeDelayT))             = 0;                              % for NaN values in vertSpeed do not add any delay
                    data.temp.value = MVP_delayCTD(data.temp.value,params.fSampl.value,-timeDelayT);
                end
                
                % - compute derived quantities
                % -- practical salinity and sound velocity (if uCTD) or
                %    conductivity (if SVPT)
                if isfield(data,{'cond'})
                    data.psal.value     = gsw_SP_from_C(data.cond.value,data.temp.value, ...
                        data.pres.value);
                    data.soundVel.value = MVP_soundspeed_from_temppsal(data.temp.value,data.psal.value, ...
                        data.pres.value,metadata.lat,1);
                    data.soundVel.unit  = 'm/s';
                else % if isfield(data,'soundVel')
                    data.psal.value = MVP_salinity_from_soundspeed(data.soundVel.value,data.temp.value, ...
                        data.pres.value,metadata.lat.value,1);
                    data.cond.value = gsw_C_from_SP(data.psal.value,data.temp.value, ...
                        data.pres.value);
                    data.cond.unit  = 'mS/cm';
                end
                data.psal.unit = 'psu';
                
                % -- absolute salinity
                asal = gsw_SA_from_SP(data.psal.value,data.pres.value, ...
                    metadata.lon.value,metadata.lat.value);
                
                % -- potential temperature
                data.ptemp.value = gsw_pt0_from_t(asal,data.temp.value,data.pres.value);
                data.ptemp.unit = 'degC';
                
                % -- conservative temperature
                ctemp = gsw_CT_from_t(asal,data.temp.value,data.pres.value);
                
                % -- density
                data.dens.value = gsw_rho(asal,ctemp,data.pres.value);
                data.dens.unit  = 'kg/m^3';
                
                % -- potential density anomaly (-1000)
                data.pdens.value = gsw_rho(asal,ctemp,0)-1000;
                data.pdens.unit  = 'kg/m^3';
                
                % - resample data onto the regular pressure grid (median filter)
                % -- build the regular pressure grid
                presGrid = (presMin : presDelta : presMax)';
                % -- resample
                fields = fieldnames(data);
                for iPres = 1 : length(presGrid)
                    ind = find(abs(data.pres.value - presGrid(iPres)) < presDelta/2);
                    for iField = 2 : length(fields)
                        dataTmp.(fields{iField}).value(iPres,1) = median(data.(fields{iField}).value(ind));
                        dataTmp.(fields{iField}).unit           = data.(fields{iField}).unit;
                    end
                end
                dataTmp.pres.value = presGrid;
                dataTmp.pres.unit  = data.pres.unit;
                data = dataTmp; clear dataTmp;
                
                % - update the data structure comprising of all the profiles
                dataAll(iProfOK) = data;
                
                % - update the metadata structure comprising of all the profiles
                metadata.transectName = sprintf('Transect_%02d',transchoice);
                metadataAll(iProfOK) = metadata;
                
                % Let's save the cast in L2 folder...
                save([pathL2 filesep fileNames(Index_Good(Unique_Good(iProf))).name(1:end-4) '_lag_1.mat'],'data','metadata','params');
                waitbar(iProfOK/nGood);
            end % Outside of "good cast" block
        end  % End of cast loop
        
        close(hWait)
        
        data = dataAll;
        metadata = metadataAll;
        params.tauPres.value   = tauPres;
        params.tauPres.unit    = 's';
        params.tauOthers.value = tauOthers;
        params.tauOthers.unit  = 's';
        params.delayT.value        = delayT;
        params.delayT.unit         = {'m','s'};
        
        savename = [pathL2 filesep sprintf('%s.mat',Transect_Names{transchoice})];
        save(savename,'metadata','data','params');
        
        disp('Done with converting the files.');
        disp('There are .mat files saved for each cast as well as one for the Transect in pathL2...');
        
        pause(3)
        
        %%%%% END OF QC METHOD NUMBER 1 %%%%%%%%%%
        
        
        %%%%% RUN THE VARIANCE METHOD %%%%%%%%
    case 2
        
        % These warnings come up during interpolation. This is to prevent
        % cluttering up the screen.
        warning('off','MATLAB:chckxy:IgnoreNaN')
        % Run variance method
        
        Best_lag_vec = [];

        iProfOK = 0;
        hWait = waitbar(0,'Calculating lag ...');
          
        fprintf('Calculating lags for all %03d chosen casts...\n',nGood);
        %disp(sprintf('Calculating lags for all %03d chosen casts...',nGood));
        %disp(char(10));
        for iProf = 1:nGood  % Trans_All(1):Trans_All(2)
            fprintf('Calculating lag for cast %03d \n',iProf);
            %disp(sprintf('Calculating lag for cast %03d ',iProf));
            load([pathL1 filesep fileNames(Index_Good(Unique_Good(iProf))).name]);
            if params.QCcast >= 0
                iProfOK = iProfOK + 1;
                % Extract downcast, run the find_lags script.
                indmaxp = find(data.pres.value == max(data.pres.value),1,'first');
                if isfield(data,{'cond'})
                    if fast_flag == 0
                        best_lag = MVP_find_lags(data.pres.value(1:indmaxp),data.cond.value(1:indmaxp),...
                            data.temp.value(1:indmaxp),metadata.lon.value,metadata.lat.value,...
                            10,0,[1 0.1],1,dens_flag);
                    else
                        best_lag = MVP_find_lags_fast(data.pres.value(1:indmaxp),data.cond.value(1:indmaxp),...
                            data.temp.value(1:indmaxp),metadata.lon.value,metadata.lat.value,...
                            10,0,1,1,dens_flag);
                    end
                else
                    if fast_flag == 0
                        best_lag = MVP_find_lags(data.pres.value(1:indmaxp),data.soundVel.value(1:indmaxp),...
                            data.temp.value(1:indmaxp),metadata.lon.value,metadata.latvalue,...
                            10,0,[1 0.1],0,dens_flag);
                    else
                        best_lag = MVP_find_lags_fast(data.pres.value(1:indmaxp),data.soundVel.value(1:indmaxp),...
                            data.temp.value(1:indmaxp),metadata.lon.value,metadata.latvalue,...
                            10,0,1,0,dens_flag);
                    end
                end

                Best_lag_vec = [Best_lag_vec nanmean(best_lag)];
                waitbar(iProfOK/nGood);
            end
            
        end
        
        close(hWait);
        
        fprintf('The best lag for these data is: %f \n',nanmean(Best_lag_vec));
        %disp(sprintf('The best lag for these data is: %f ',nanmean(Best_lag_vec)));
        hist(Best_lag_vec,50);
        xlabel('Lag (scans) ');
        ylabel('Number of occurrences ');
        title('Histogram of best lags for this transect.');
        drawnow
        pause(2)
        disp('Now making the correction for the casts...');
        
        Best_lag = nanmean(Best_lag_vec); % We'll use the nearest 0.1 scan...
        Best_lag = round(Best_lag*10)/10;
        
        iProfOK = 0;
        hWait = waitbar(0,'Processing L1 files ...');
                
        for iProf = 1:nGood  % Trans_All(1):Trans_All(2)
            load([pathL1 filesep fileNames(Index_Good(Unique_Good(iProf))).name]);
            if params.QCcast >= 0
                
                % For all the good casts...
                %%% AD condition added to analyze data from only uCTD (when the SSFFF fish is used)
                if isfield(data,{'LOPCcounts'})
%                if ~strcmp(instrType,'uCTD') %% LR change Cruise_Setup.instrType in instrType from MVP_cruise_params
                    data = rmfield(data,{'LOPCcounts'});
                end
                iProfOK = iProfOK+1;
                
                % Check to see if ADCcounts is present, and if so, convert
                % it to Voltage (Conversion parameter specified in
                % MVP_cruise_params.m)
                if isfield(data,{'ADCcounts'})
                    data.analog.value = (data.ADCcounts.value).*volt_input / 2^(bit_input);
                    data.analog.unit = 'V';
                    data = rmfield(data,{'ADCcounts'});
                end

                % - low-pass filter data, if possible
                if filter_flag %AD 9999 instead of 16
                    if exist('butter')== 2
                        % -- pressure
                        fCutPres = (1/tauPres)/(params.fSampl.value/2);
                        [b,a] = butter(3,fCutPres,'low');
                        data.pres.value = MVP_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 = MVP_filter_updn(b,a,data.temp.value,1);
                        
                        if isfield(data,{'cond'})
                            data.cond.value = MVP_filter_updn(b,a,data.cond.value,1);
                        end
                        if isfield(data,{'soundVel'})
                            data.soundVel.value = MVP_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);
                    % Below keeps the upcast instead of the downcast...
                    %data.(fields{iField}).value = data.(fields{iField}).value(ind+1:end);
                end
                
                % Delay the Temperature
                % Interpolate the profile, to the 0.1 scan precision...
                prof_len = length(data.pres.value);
                interp_space = linspace(1,prof_len,prof_len*10);
                interp_short = interp_space(1:end-floor(Best_lag*10));
                % Interpolate the data fields, perform the lag on
                % temperatures, and then interpolate back.
                
                % In case there is now analog instead of ADCcounts...
                fields = fieldnames(data);
                
                for iField = 1 : length(fields)
                    dataInt.(fields{iField}).value = pchip((1:prof_len),data.(fields{iField}).value,interp_space);
                    if strcmp((fields{iField}),'temp')
                        dataInt.(fields{iField}).value = dataInt.(fields{iField}).value(1+floor(Best_lag*10):end);
                    else
                        dataInt.(fields{iField}).value = dataInt.(fields{iField}).value(1:end-floor(Best_lag*10));
                    end
                    data.(fields{iField}).value = pchip(interp_short,dataInt.(fields{iField}).value,(1:prof_len-ceil(Best_lag)));
                end
                
                % - compute derived quantities
                % -- practical salinity and sound velocity (if uCTD) or
                %    conductivity (if SVPT)
                if isfield(data,{'cond'})
                    data.psal.value     = gsw_SP_from_C(data.cond.value,data.temp.value, ...
                        data.pres.value);
                    data.soundVel.value = MVP_soundspeed_from_temppsal(data.temp.value,data.psal.value, ...
                        data.pres.value,metadata.lat,1);
                    data.soundVel.unit  = 'm/s';
                else
                    data.psal.value = MVP_salinity_from_soundspeed(data.soundVel.value,data.temp.value, ...
                        data.pres.value,metadata.lat.value,1);
                    data.cond.value = gsw_C_from_SP(data.psal.value,data.temp.value, ...
                        data.pres.value);
                    data.cond.unit  = 'mS/cm';
                end
                data.psal.unit = 'psu';
                
                % -- absolute salinity
                asal = gsw_SA_from_SP(data.psal.value,data.pres.value, ...
                    metadata.lon.value,metadata.lat.value);
                
                % -- potential temperature
                data.ptemp.value = gsw_pt0_from_t(asal,data.temp.value,data.pres.value);
                data.ptemp.unit = 'degC';
                
                % -- conservative temperature
                ctemp = gsw_CT_from_t(asal,data.temp.value,data.pres.value);
                
                % -- density
                data.dens.value = gsw_rho(asal,ctemp,data.pres.value);
                data.dens.unit  = 'kg/m^3';
                
                % -- potential density anomaly (-1000)
                data.pdens.value = gsw_rho(asal,ctemp,0)-1000;
                data.pdens.unit  = 'kg/m^3';
                
                % - resample data onto the regular pressure grid (median filter)
                % -- first choose the right pressure step
                presGrid = (presMin : presDelta : presMax)';
                % -- resample
                fields = fieldnames(data);
                for iPres = 1 : length(presGrid)
                    ind = find(abs(data.pres.value - presGrid(iPres)) < presDelta/2);
                    for iField = 2 : length(fields)
                        dataTmp.(fields{iField}).value(iPres,1) = median(data.(fields{iField}).value(ind));
                        dataTmp.(fields{iField}).unit           = data.(fields{iField}).unit;
                    end
                end
                dataTmp.pres.value = presGrid;
                dataTmp.pres.unit  = data.pres.unit;
                data = dataTmp; clear dataTmp;
                %data = orderfields(data,{'pres','temp','ptemp','cond','soundVel','psal','dens','pdens','chla','vertSpeed'});
                
                % Incorporate the lag into the params structure...
                params.delayT.value  = Best_lag;
                params.delayT.unit   = {'scan'};
                
                % - update the data structure comprising of all the profiles
                dataAll(iProfOK) = data;
                
                % - update the metadata structure comprising of all the profiles
                metadata.transectName = sprintf('Transect_%02d',transchoice);
                metadataAll(iProfOK) = metadata;
                

                % Let's save the cast in L2 folder...
                save([pathL2 filesep fileNames(Index_Good(Unique_Good(iProf))).name(1:end-4) sprintf('_lag2_dens_%01d_fast_%01d.mat',dens_flag,fast_flag)],'data','metadata','params');
                waitbar(iProfOK/nGood);
                
            end % Outside of "good cast" block
        end  % End of cast loop
        
        close(hWait)
        data = dataAll;
        metadata = metadataAll;
        
        params.tauPres.value   = tauPres;
        params.tauPres.unit    = 's';
        params.tauOthers.value = tauOthers;
        params.tauOthers.unit  = 's';
        params.delayT.value    = Best_lag_vec;
        params.delayT.unit         = {'scan'};

        
        savename = [pathL2 filesep sprintf('%s.mat',Transect_Names{transchoice})];
        save(savename,'metadata','data','params');
        
        disp('Done with converting the files.');
        pause(3)
        
        disp('There are .mat files saved for each cast as well as one for the Transect in pathL2...');
        
        %%%%%%%%%%%% END OF VARIATIONAL LAG METHOD %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        
end

%%%%%%%%%%% END OF LAG QC %%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%% PLOTTING! :D %%%%%%%
% Now to plot the transect just processed

disp('Now we are ready to plot the data...');
pause(2)

% Plotting parameters...
MVP_plot_transect(pathL2,savename(length(pathL2)+2:end-4),color_axes);

disp('When finished with the plot, please press any key to exit...')
pause

disp('End of L1toL2. ');
disp('Exiting....');
pause(1)
clear all; clc; close all;


% End of L1toL2.m

% BTW, to run the MVP_plot_transect function from nothing... do:
% MVP_cruise_params
% MVP_plot_transect(pathL2,sprintf('Transect_%02d_%01d',transchoice,qc_user),color_axes);
%
