
% 
% Loads a MVP .mat file generated with L1toL2.m and a RosetteCTD .mat file
% generated with functions/CNVtoMAT.m.
% 
% Chlorophyll-a profiles are then compared and fitting coefficients computed
% through a Least Squares Error minimization.
% 
% The correction model is the following:
% 
%              chla = c0 + c1*chla,
% 
% where actually c0 is the Scale Factor and -c0/c1 are the Dark Counts.
% 
% To find the correspondance between MVP and RosetteCTD casts, see the
% document Report_Utilization_MVP_OUTPACE_2015.pdf
% 
% 
% v0
%   Lucio Bellomo, 25 June 2015
% 


clear all; clc;


%% Parameters to be set by the user
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MVP file name  (.mat)
MVPfile = 'Test_1';
% MVPfile = 'Test_5';
% MVPfile = 'Test_7';

% CTD file name
CTDfile = 'out_c_021';
% CTDfile = 'out_c_080';
% CTDfile = 'out_c_104';

% Minimum pressure at which the fitting should start (to avoid using the
% upper part of the water column, where variability is higher due to e.g.
% internal waves)
presMin = 0;

% Maximum pressure at which the fitting should stop (if no max pres, put Inf)
presMax = Inf;

% Folder where the MVP .mat are
pathMVP = '/media/NeNe/DatiVari/AtalanteDisk/OUTPACE_MVP/MVP_DATA/All_Data/L2/';

% Folder where the RosetteCDT .mat are
pathCTD = '/media/NeNe/DatiVari/AtalanteDisk/OUTPACE_MVP/MVP_DATA/All_Data/PostCal/RosetteCTD/';

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


%% Load CTD data and fetch the variables to be used for the fitting, plus some metadata
load([pathCTD filesep CTDfile '.mat']);

% Extract the data
CTDdate = data(1).value;
CTDlon  = data(2).value;
CTDlat  = data(3).value;
CTDpres = data(4).value(:);
for iField = 5 : length(data)
    if ~isempty(strfind(data(iField).name,'Fluorescence'))
        CTDchla = data(iField).value;
    end
end

% Only downcast
[~, indBottom] = max(CTDpres);
CTDpres = CTDpres(1:indBottom);
CTDchla = CTDchla(1:indBottom);


%% Load MVP data and fetch the variables to be used for the fitting, plus some metadata
load([pathMVP filesep MVPfile '.mat']);

nProf = length(metadata);
for iProf = 1 : nProf
    MVPdate(iProf) = metadata(iProf).date.value;
    MVPlon(iProf)  = metadata(iProf).lon.value;
    MVPlat(iProf)  = metadata(iProf).lat.value;
    MVPpres(:,iProf) = data(iProf).pres.value;
    MVPchla(:,iProf) = data(iProf).chla.value;
    chlaUnit{iProf} = data(iProf).chla.unit;
end

% Check on the units
for iProf = 1 : nProf
    if ~strcmp(chlaUnit{iProf},'V')
        errMsg = ['The unit of chla for the MVP is not Volt! The way the ' ...
                  'calibration is handled here and in L1toL2 requires chla ' ...
                  'to be in Volt. To have this, you have to run again L1toL2 ' ...
                  'with the parameter convertChla = [NaN NaN].'];
        error(errMsg);
    end
end


%% Some statistics
timeDiff = (CTDdate - MVPdate)*24*60;   % min
disp(['The difference in time between the two profiles is ' num2str(round(timeDiff)) ' min']);

[dx, dy] = lonlat2km(CTDlon,CTDlat,MVPlon,MVPlat);
dist = sqrt(dx.^2 + dy.^2)*1000;
disp(['The distance between the two profiles is ' num2str(round(dist)) ' m']);
disp(blanks(1)');


%% First plot
hFig = figure; hold on; grid on;
plot(CTDchla,CTDpres,'-b.');
for iProf = 1 : nProf
    plot(MVPchla,MVPpres(:,iProf),'-rx');
end
xlabel('Chlorophyll-a concentration (ug/l)');
ylabel('Pressure (dbar)');
set(gca,'YDir','rev');


%% Do you want to use all the MVP profiles (if more than 1) or want to select one or more?
if nProf > 1
    msg = ['Type 0 if you want to work with all MVP profiles;\ntype the profile ' ...
           'number(s) if you only a subset (ex.: [1 3])\n'];
    iProf = input(sprintf(msg));
    nProf = length(iProf);

    % Fetch the profiles
    MVPdate = MVPdate(iProf);
    MVPlon  = MVPlon(iProf);
    MVPlat  = MVPlat(iProf);
    MVPpres = MVPpres(:,iProf);
    MVPchla = MVPchla(:,iProf);

    % Merge the profiles
    [MVPpres, ind] = sort(MVPpres(:));
    MVPchla = MVPchla(ind);
end


%% Interpolate the CTD data onto the MVP ones, since the latter are less precise and resolved
CTDpres_interp = MVPpres;
CTDchla_interp = interp1q(CTDpres,CTDchla,MVPpres);


%% Plot
close(hFig);

% Chlorophyll-a
hChla = figure; hold on; grid on;
plot(CTDchla_interp,CTDpres_interp,'-b.');
plot(MVPchla,MVPpres,' rx');
set(gca,'YDir','rev');
xlabel('Chlorophyll-a concentration (ug/l)');
ylabel('Pressure (dbar)');


%% Post-calibration
% Compute the coefficients
indOK = ~isnan(CTDchla_interp) & ~isnan(MVPchla);
A = [sum(indOK),              sum(MVPchla(indOK));
     sum(MVPchla(indOK)), sum(MVPchla(indOK).^2)];
b = [sum(CTDchla_interp(indOK)); ...
     sum(CTDchla_interp(indOK).*MVPchla(indOK))];
coefs = A \ b;

% Apply the coefficients
MVPchlaCorr = coefs(1) + coefs(2)*MVPchla;
disp('The coefficients are:');
for iCoef = 1 : length(coefs)
    disp(['c' num2str(iCoef-1) ' = ' num2str(coefs(iCoef),'%.4f')]);
end
disp('Or, equivalently:');
darkCounts = -coefs(1)/coefs(2);
scaleFactor = coefs(2);
disp(['darkCounts = '  num2str(darkCounts,'%.4f')]);
disp(['scaleFactor = ' num2str(scaleFactor,'%.4f')]);

% Update the plot
figure(hChla);
plot(MVPchlaCorr,MVPpres,'m');
legend('CTD','MVP','MVP post-cal');
