
% 
% Loads the MVP data from the .mat files generated with L0toL1.m and
% processes them to generate a .mat and a .nc file for each group of
% transects (LDA_Cross, LDA_Before, LDA_After, ..., Test_1, Test_2, etc.)
% 
% N.B.: If 'LDA' or 'Test' are chosen, then the .nc file will not be 
%       written, since within LDA the set of processing parameters changes 
%       (uCTD, then SVPT) and the mat2nc.m is not ready for this yet
% 
% 
% Data processing consists in the following:
% 1) remove profiles where a minimum theshold pressure has not been reached
% 2) retrieve chla concentration from the analog output of the MVP
% 3) low-pass filtering
% 4) T, C, and SV post-calibration (if available)
% 5) alignement of T with cond/SV due to the different response times
% 6) down/up casts separation -> only the downcast is recorded
% 7) rebinning onto a regular pressure grid
% 8) alignment of chla with hydrological variables due to the different
%    position of the sensors in the MSFFF (fluo is in the nose, whereas
%    the others are in the back or in the middle) (only for the downcast, 
%    since for the upcast the MSFFF is supposed to be horizontal)
% 
% *** Some of these processing steps have been inspired by 
%     Ullmand and Hebert, JAOT, 31, 2014
% 
% 
% v0
%   Lucio Bellomo, 29 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, LDA_Cross, LDA_Before, LDA_After, LDB, LDB_Cross, LDB_Before,
%   LDB_After, LDC_Cross (or LDC), Test, Test_1..7
groupOfTransects = {'Test_1'};

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

% AnalogOutput -> Chlorophyll-a concentration conversion
convertChla = [0.041 4.587]; % [darkCounts scaleFactor]

% 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; % these values come from Ullman and Hebert, JAOT, 31, 2014

% Post-calibration coefs (NaN if no post-calibration)
calibT_uCTD    = [NaN NaN NaN]; % [offset coefTemp coefPres]
calibT_SVPT_1  = [NaN NaN NaN]; % [offset coefTemp coefPres]
calibT_SVPT_2  = [NaN NaN NaN]; % [offset coefTemp coefPres]
calibC_uCTD    = [13.9002 0.6578 -0.0051]; % [offset coefCond coefPres]
calibSV_SVPT_1 = [NaN NaN NaN]; % [offset coefSV coefPres]
calibSV_SVPT_2 = [NaN NaN NaN]; % [offset coefSV coefPres]

% Alignment: delay of T wrt cond or SV (NaN if no delay)
delayTparams_uCTD   = [-0.13528 0.72684];
delayTparams_SVPT_1 = [1.03965 0.72826];
delayTparams_SVPT_2 = [0.41884 0.84381];

% Resampling on a regular pressure grid
presMin        = 10;
presMax        = 350;
presDelta_uCTD = 5; % with the uCTD the sampling rate was 1 Hz!!!
presDelta_SVPT = 2;

% Alignment: distance between fluorometer and CTD/SVPT (m)
deltaPresFluo_uCTD   = 1.05;    % LDA cross and before
deltaPresFluo_SVPT_1 = 1.05;    % LDA after
deltaPresFluo_SVPT_2 = 0.52;    % LDB and LDC

% 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/';

% Folder where the .mat (final) and .nc are to be written
pathL2 = '/media/NeNe/DatiVari/AtalanteDisk/OUTPACE_MVP/MVP_DATA/All_Data/L2/';

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



%% 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
% Group by group
for iGrp = 1 : length(groupOfTransects)
    
    
    % Initialization    
    iProfOK = 0;
    
    
    % 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 LOPCcounts as we do not know yet how to investigare the assessment
                data = rmfield(data,'LOPCcounts');
                
                % - check that the minimum pressure value has been reached
                if max(data.pres.value) >= minValidPres
                    iProfOK = iProfOK + 1;
                else
                    continue;
                end
                
                % - remove salinity, if present, since it will be computed
                %   more properly later
                if isfield(data,'psal')
                    data = rmfield(data,'psal');
                end
    
                % - transform the ADC voltage into chla concentration and
                %   remove analogOutput from the data structure
                if ~all(isnan(convertChla))
                    data.chla.value = convertChla(2) * (data.analogOutput.value - convertChla(1));
                    data.chla.unit  = 'mg/m^3';
                else
                    data.chla = data.analogOutput;
                end
                data = rmfield(data,'analogOutput');
                
                % - 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,'soundVel')
                        data.soundVel.value = filter_updn(b,a,data.soundVel.value,1);
                    end
                    data.chla.value = filter_updn(b,a,data.chla.value,1);
                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';
                
                % - post-calibration of T, cond, and soundVel
                %   For the parameter 'P', the model is of the kind:
                %   P <- a0 + a1*P + a2*pres
                % -- first choose the right set of parameters
                if ~isempty(strfind(groupOfTransects{iGrp},'Test')) % if it's a post-cal test profile, do not calibrate
                    calibT  = [NaN NaN NaN];
                    calibC  = [NaN NaN NaN];
                    calibSV = [NaN NaN NaN];
                elseif strcmp(params.instrType,'uCTD+fluo')                                         % LDA_Cross and LDA_Before
                    calibT = calibT_uCTD;
                    calibC = calibC_uCTD;
                elseif metadata.date.value < datenum('06 Mar 2015 12:00:00','dd mmm yyyy HH:MM:SS') % LDA_After
                    calibT  = calibT_SVPT_1;
                    calibSV = calibSV_SVPT_1;
                else                                                                                % from LDB_Cross till the end
                    calibT  = calibT_SVPT_2;
                    calibSV = calibSV_SVPT_2;
                end
                % -- then apply it
                if ~all(isnan(calibT))
                    data.temp.value = calibT(1) + ...
                                      calibT(2)*data.temp.value + ...
                                      calibT(3)*data.pres.value;
                end
                if strcmp(params.instrType,'uCTD+fluo')
                    if ~all(isnan(calibC))
                        data.cond.value = calibC(1) + ...
                                          calibC(2)*data.cond.value + ...
                                          calibC(3)*data.pres.value;
                    end
                else
                    if ~all(isnan(calibSV))
                        data.soundVel.value = calibSV(1) + ...
                                              calibSV(2)*data.soundVel.value + ...
                                              calibSV(3)*data.pres.value;
                    end
                end                
    
                % - alignment of hydrological parameters
                % -- first choose the right set of parameters
                if ~isempty(strfind(groupOfTransects{iGrp},'Test')) % if it's a post-cal test profile, align only if the cast has been done in freewheel mode
                    if strcmp(metadata.fileName,'MVP_2015-02-25_110647')                                            % Test_1
                        delayT = delayTparams_uCTD;
                    elseif (metadata.date.value >= datenum('06 Mar 2015 23:08:00','dd mmm yyyy HH:MM:SS') && ...    % Test_5
                            metadata.date.value <= datenum('07 Mar 2015 00:01:59','dd mmm yyyy HH:MM:SS')) || ...
                            metadata.date.value >= datenum('14 Mar 2015 22:52:00','dd mmm yyyy HH:MM:SS')           % Test_7
                        delayT = delayTparams_SVPT_2;
                    else
                        delayT = [NaN NaN];                                                                         % all other tests
                    end
                elseif strcmp(params.instrType,'uCTD+fluo')                                         % LDA_Cross and LDA_Before
                    delayT = delayTparams_uCTD;
                elseif metadata.date.value < datenum('06 Mar 2015 12:00:00','dd mmm yyyy HH:MM:SS') % LDA_After
                    delayT = delayTparams_SVPT_1;
                else                                                                                % from LDB_Cross till the end
                    delayT = delayTparams_SVPT_2;
                end
                % -- then apply it
                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 = delayCTD(data.temp.value,params.fSampl.value,-timeDelayT);
                end
        
                % - compute derived quantities
                % -- practical salinity and sound velocity (if uCTD) or
                %    conductivity (if SVPT)
                if strcmp(params.instrType,'uCTD+fluo')
                    data.psal.value     = gsw_SP_from_C(data.cond.value,data.temp.value, ...
                                                        data.pres.value);
                    data.soundVel.value = soundspeed_from_temppsal(data.temp.value,data.psal.value, ...
                                                                   data.pres.value,metadata.lat,1);
                    data.soundVel.unit  = 'm/s';
                else
                    data.psal.value = 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';
    
                % - 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);
                end
                
                % - for the fluorometer align the measurement with pressure
                % -- first choose the right set of parameters
                if ~isempty(strfind(groupOfTransects{iGrp},'Test')) % if it's a post-cal test profile, align only if the cast has been done in freewheel mode
                    if strcmp(metadata.fileName,'MVP_2015-02-25_110647')                                            % Test_1
                        deltaPresFluo = deltaPresFluo_uCTD;
                    elseif (metadata.date.value >= datenum('06 Mar 2015 23:08:00','dd mmm yyyy HH:MM:SS') && ...    % Test_5
                            metadata.date.value <= datenum('07 Mar 2015 00:01:59','dd mmm yyyy HH:MM:SS')) || ...
                            metadata.date.value >= datenum('14 Mar 2015 22:52:00','dd mmm yyyy HH:MM:SS')           % Test_7
                        deltaPresFluo = deltaPresFluo_SVPT_2;
                    else
                        deltaPresFluo = NaN;                                                                        % all other tests
                    end
                elseif strcmp(params.instrType,'uCTD+fluo')                                         % LDA_Cross and LDA_Before
                    deltaPresFluo = deltaPresFluo_uCTD;
                elseif metadata.date.value < datenum('06 Mar 2015 12:00:00','dd mmm yyyy HH:MM:SS') % LDA_After
                    deltaPresFluo = deltaPresFluo_SVPT_1;
                else                                                                                % from LDB_Cross till the end
                    deltaPresFluo = deltaPresFluo_SVPT_2;
                end
                % -- then apply it
                if ~isnan(deltaPresFluo)
                    presFluo = data.pres.value + deltaPresFluo;    % actual pressure where the fluorometer is placed
                    data.chla.value = interp1q(presFluo,data.chla.value,data.pres.value);
                end

                % - resample data onto the regular pressure grid (median filter)
                % -- first choose the right pressure step
                if strcmp(params.instrType,'uCTD+fluo')
                    presDelta = presDelta_uCTD;
                else
                    presDelta = presDelta_SVPT;
                end
                % -- build the regular pressure grid
                presGrid = [presMin : presDelta : presMax]';
                % -- resample
                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'});
                
                % - update the data structure comprising of all the profiles
                dataAll(iProfOK) = data;
                
                % - update the metadata structure comprising of all the profiles
                metadata.transectName = group(iGrp).transectName{iTrans};
                metadataAll(iProfOK) = metadata;
                
            end % date inside the transect date range
            
            waitbar(iProf/nProf);
            
        end % profiles
        
        
        % Close waitbar
        close(hWait);
        pause(0.5);
        
    end % transects
    
    
    % Save .mat and .nc for each group of transects
    fileNameMat = [pathL2 filesep groupOfTransects{iGrp} '.mat'];
    
    
    % Processing parameters
    params.minValidPres.value = minValidPres;
    params.minValidPres.unit  = 'dbar';
    params.convertChla.value = convertChla;
    params.convertChla.unit  = {'V','mg/m^3/V'};
    params.tauPres.value   = tauPres;
    params.tauPres.unit    = 's';
    params.tauOthers.value = tauOthers;
    params.tauOthers.unit  = 's';
    params.calibT.value    = calibT;
    params.calibT.unit     = {'degC','','degC/dbar'};
    if strcmp(params.instrType,'uCTD+fluo')
        params.calibC.value = calibC_uCTD;
        params.calibC.unit  = {'mS/cm','','mS/cm/dbar'};
    else
        params.calibSV.value = calibSV;
        params.calibSV.unit  = {'m/s','','m/s/dbar'};
    end
    params.delayT.value        = delayT;
    params.delayT.unit         = {'m','s'};
    params.deltaPresFluo.value = deltaPresFluo;
    params.deltaPresFluo.unit  = 'm';
    
    
    % Metadata and data
    metadata = metadataAll;
    data     = dataAll;
    
    
    % Save .mat
    save(fileNameMat,'params','metadata','data');
    disp(blanks(1)');
    disp(['*** Done - Created ' groupOfTransects{iGrp} '.mat ***']);
    
    
    % Save NetCDF
    if strcmp(groupOfTransects{iGrp},'LDA') || strcmp(groupOfTransects{iGrp},'Test')
        msg = ['Couldn''t create ' groupOfTransects{iGrp} '.nc because the ' ...
               'processing parameters change for each test and the .mat->NetCDF ' ...
               'conversion isn''t\nready yet to cope with this (it''s in the ' ...
               'global attributes that the parameters should be made profile' ...
               '-dependent)!\n'];
        fprintf(msg);
    else
        mat2nc(fileNameMat);
        disp(['*** Done - Created ' groupOfTransects{iGrp} '.nc ***']);
    end
    disp(blanks(1)');

end % groups

