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

% CTD file name
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,'Sound Velocity'))
        CTDsoundVel = 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);
CTDsoundVel = CTDsoundVel(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;
    MVPsoundVel(:,iProf) = data(iProf).soundVel.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(CTDsoundVel,CTDpres,'-b.');
for iProf = 1 : nProf
    plot(MVPsoundVel,MVPpres(:,iProf),'-rx');
end
xlabel('Sound Velocity (m/s)');
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);
    MVPsoundVel = MVPsoundVel(:,iProf);
    MVPpsal     = MVPpsal(:,iProf);

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


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


%% Plot
close(hFig);

% Sound Velocity
hSV = figure;
h1 = subplot(1,2,1); hold on; grid on;
plot(CTDsoundVel_interp,CTDpres_interp,'-b.');
plot(MVPsoundVel,MVPpres,' rx');
set(gca,'YDir','rev');
xlabel('Sound Velocity (m/s)');
ylabel('Pressure (dbar)');
h2 = subplot(1,2,2); hold on; grid on;
plot(MVPsoundVel./CTDsoundVel_interp,MVPpres);
set(gca,'YDir','rev');
xlabel('Sound Velocity 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(CTDsoundVel_interp) & ~isnan(MVPsoundVel);
A = [sum(indOK),              sum(MVPsoundVel(indOK)),                 sum(MVPpres(indOK));
     sum(MVPsoundVel(indOK)), sum(MVPsoundVel(indOK).^2),              sum(MVPsoundVel(indOK).*MVPpres(indOK));
     sum(MVPpres(indOK)),     sum(MVPsoundVel(indOK).*MVPpres(indOK)), sum(MVPpres(indOK).^2)];
b = [sum(CTDsoundVel_interp(indOK)); ...
     sum(CTDsoundVel_interp(indOK).*MVPsoundVel(indOK));...
     sum(CTDsoundVel_interp(indOK).*MVPpres(indOK))];
% coefs = A(1:2,1:2) \ b(1:2);
coefs = A \ b;

% Apply the coefficients
% MVPtempCorr = coefs(1) + coefs(2)*MVPtemp;
MVPsoundVelCorr = coefs(1) + coefs(2)*MVPsoundVel + 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(hSV);
subplot(1,2,1);
plot(MVPsoundVelCorr,MVPpres,'m');
legend('CTD','MVP','MVP post-cal');
subplot(1,2,2);
plot(MVPsoundVelCorr./CTDsoundVel_interp,MVPpres,'m');
legend('MVP/CTD','MVP post-cal/CTD');

% Compute the salinity
MVPpsalCorr = salinity_from_soundspeed(MVPsoundVelCorr,MVPtemp,MVPpres,MVPlat,1);

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

