
% 
% Loads a MVP .mat file generated with L1toL2.m and a RosetteCTD .mat file
% generated with functions/CNVtoMAT.m.
% 
% Conductivity profiles are then compared and fitting coefficients computed
% through a Least Squares Error minimization.
% 
% The correction model is the following:
% 
%                   tempCorr = c0 + c1*temp + c2*pres
% 
% 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';

% CTD file name
CTDfile = 'out_c_021';

% 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,'Conductivity '))  % space needed at the end!
        CTDcond = data(iField).value;
    end
    if ~isempty(strfind(data(iField).name,'Salinity, Practical '))  % space needed at the end!
        CTDpsal = data(iField).value;
    end
end

% Only downcast
[~, indBottom] = max(CTDpres);
CTDpres = CTDpres(1:indBottom);
CTDcond = CTDcond(1:indBottom);
CTDpsal = CTDpsal(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;
    MVPtemp(:,iProf) = data(iProf).temp.value;
    MVPcond(:,iProf) = data(iProf).cond.value;
    MVPpsal(:,iProf) = data(iProf).psal.value;
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(CTDcond,CTDpres,'-b.');
for iProf = 1 : nProf
    plot(MVPcond,MVPpres(:,iProf),'-rx');
end
xlabel('Conductivity (mS/cm)');
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);
    MVPtemp = MVPtemp(:,iProf);
    MVPcond = MVPcond(:,iProf);
    MVPpsal = MVPpsal(:,iProf);

    % Merge the profiles
    [MVPpres, ind] = sort(MVPpres(:));
    MVPtemp = MVPtemp(ind);
    MVPcond = MVPcond(ind);
    MVPpsal = MVPpsal(ind);
    MVPother = MVPother(ind);
end


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


%% Plot
close(hFig);

% Conductivity
hTemp = figure;
h1 = subplot(1,2,1); hold on; grid on;
plot(CTDcond_interp,CTDpres_interp,'-b.');
plot(MVPcond,MVPpres,' rx');
set(gca,'YDir','rev');
xlabel('Conductivity (mS/cm)');
ylabel('Pressure (dbar)');
h2 = subplot(1,2,2); hold on; grid on;
plot(MVPcond./CTDcond_interp,MVPpres);
set(gca,'YDir','rev');
xlabel('Conductivity ratio MVP / CTD');
ylabel('Pressure (dbar)');
linkaxes([h1 h2],'y');

% Salinity
hPsal = figure;
h1 = subplot(1,2,1); hold on; grid on;
plot(CTDpsal_interp,CTDpres_interp,'-b.');
plot(MVPpsal,MVPpres,' rx');
set(gca,'YDir','rev');
ylabel('Pressure (dbar)');
xlabel('Practical salinity (psu)');
h2 = subplot(1,2,2); hold on; grid on;
plot(MVPpsal./CTDpsal_interp,MVPpres);
set(gca,'YDir','rev');
xlabel('Salinity ratio MVP / CTD');
ylabel('Pressure (dbar)');
linkaxes([h1 h2],'y');


%% Post-calibration
% Compute the coefficients
indOK = ~isnan(CTDcond_interp) & ~isnan(MVPcond);
A = [sum(indOK),          sum(MVPcond(indOK)),                 sum(MVPpres(indOK));
     sum(MVPcond(indOK)), sum(MVPcond(indOK).^2),              sum(MVPcond(indOK).*MVPpres(indOK));
     sum(MVPpres(indOK)), sum(MVPcond(indOK).*MVPpres(indOK)), sum(MVPpres(indOK).^2)];
b = [sum(CTDcond_interp(indOK)); ...
     sum(CTDcond_interp(indOK).*MVPcond(indOK));...
     sum(CTDcond_interp(indOK).*MVPpres(indOK))];
% coefs = A(1:2,1:2) \ b(1:2);
coefs = A \ b;

% Apply the coefficients
% MVPcondCorr = coefs(1) + coefs(2)*MVPcond;
MVPcondCorr = coefs(1) + coefs(2)*MVPcond + coefs(3)*MVPpres;
disp('The coefficients are:');
for iCoef = 1 : length(coefs)
    disp(['c' num2str(iCoef-1) ' = ' num2str(coefs(iCoef),'%.4f')]);
end

% Update the plot
figure(hTemp);
subplot(1,2,1);
plot(MVPcondCorr,MVPpres,'m');
legend('CTD','MVP','MVP post-cal');
subplot(1,2,2);
plot(MVPcondCorr./CTDcond_interp,MVPpres,'m');
legend('MVP/CTD','MVP post-cal/CTD');

% Compute the salinity
if all(MVPdate < datenum('28 Feb 2015 00:00:00','dd mmm yyyy HH:MM:SS'))
    MVPpsalCorr = gsw_SP_from_C(MVPcondCorr,MVPtemp,MVPpres);
elseif all(MVPdate >= datenum('28 Feb 2015 00:00:00','dd mmm yyyy HH:MM:SS'))
    MVPpsalCorr = salinity_from_soundspeed(MVPother,MVPcondCorr,MVPpres,MVPlat,1);
else
    error('The profiles are partly from the uCTD and partly from the SVPT - Impossible to compute the salinity!');
end

% Update the salinity plot
figure(hPsal);
subplot(1,2,1);
plot(MVPpsalCorr,MVPpres,'m');
legend('CTD','MVP','MVP post-cal');
subplot(1,2,2);
plot(MVPpsalCorr./CTDpsal_interp,MVPpres,'m');
legend('MVP/CTD','MVP post-cal/CTD');

