clc
clear
close all

% Get TSG and bucket data
% -----------------------
% tsg    = getappdata( hMainFig, 'tsg_data');
% bucket = getappdata( hMainFig, 'bucket_data');

% Lecture fichier TSG
% -------------------
filename = 'F:\work\M_TsgQc\tsg_data\past0601.txt';
[tsg, error] = tsg_read( filename);

% Lecture fichier bucket
% -----------------------
filename = 'F:\work\M_TsgQc\tsg_data\past0601.btl';
[bucket, error] = tsg_readBucket( filename);

% Salinity in one 1 hour interval, should not depart the average for more
% than SAL_STD_LIM standard deviation
% -----------------------------------------------------------------------
SAL_STD_LIM = 0.1;

% 1 hour interval expressed in MATLAB serial Number
% -------------------------------------------------
INTERVAL_SMOOTHING = datenum(0, 0, 0, 1, 0 , 0);

% dt between 2 tsg measurements
% -----------------------------
TSG_SAMPLING_TIME = datenum(0, 0, 0, 0, 5 , 0);
% ou plus general
TSG_SAMPLING_TIME = tsg.TIME(2) - tsg.TIME(1);

% Correction Windows in days
% --------------------------
TIME_WINDOWS = 10;


% *************************************************************************
%                             First Step
%
% Running average over INTERVAL_SMOOTHING hour.
% Salinity in one INTERVAL_SMOOTHING hour interval, should not depart the
% average for more than SAL_STD_LIM
%
% *************************************************************************

[psal_smooth, nval] = ...
      tsg_moveaverage(tsg.TIME, tsg.PSAL, INTERVAL_SMOOTHING, SAL_STD_LIM);

  
% In the future we could use additional samples (i.e. CTD data)
% These data should be merged with the bucket measurement.
% All the samples would be stored in a structure 'sample'
% Function to be built
% ---------------------------------------------------------------------
sample.TIME        = bucket.TIME_WS;
sample.LATITUDE    = bucket.LATITUDE_WS;
sample.LONGITUDE   = bucket.LONGITUDE_WS;
sample.PSAL        = bucket.PSAL_WS;
sample.PSAL_QC     = bucket.PSAL_QC_WS;
sample.PSAL_DIF    = zeros(size(sample.PSAL));
sample.PSAL_SMOOTH = zeros(size(sample.PSAL));
  
% *************************************************************************
%                             Second Step
% 
% Co-location of samples and TSG measurements
% Compute the sample-TSG differences

sample = tsg_diffTsgSample(tsg, psal_smooth, sample, TSG_SAMPLING_TIME);


% Exclusion des points aberrants (ecart trop grand)
% -------------------------------------------------
%A = find( sal_diff  < (median(sal_diff) + 3 * std(sal_diff)));
%B = find( sal_diff  > (median(sal_diff) - 3 * std(sal_diff)));
%C = find( sal_diff == (median(sal_diff) + 3 * std(sal_diff)));


% *************************************************************************

% Calcul des valeurs medianes de correction: achaque point de
% comparaison, on attribue la mediane des corrections dans une fenetre de
% 10 jours. Puis on interpole entre chaque mediane pour avoir la correction
% a appliquer en chaque point de donnee TSG

% Fenetre de temps = 10 jours pour le calcul des valeurs medianes de
% correction
% ------------------------------------------------------------------

COR_mediane    = [0 0];
erreur_mediane = [0 0];

for i = 1:length(sample)
    
    is = find( sample.TIME > sample.TIME(i) - TIME_WINDOWS/2 & ...
               sample.TIME < sample.TIME(i) + TIME_WINDOWS/2);
           
    if ~isempty(is)
        COR_mediane    = [COR_mediane;...
                         [ sample.TIME(i) median(sample.PSAL_DIF(is))]];
        erreur_mediane = [erreur_mediane;...
                         [std(sample.PSAL_DIF(is))/sqrt(length(is)) length(is)]];
    end
end


COR_mediane    = COR_mediane(2:length(COR_mediane),:);
erreur_mediane = erreur_mediane(2:length(erreur_mediane),:);

erreur_mediane( find(erreur_mediane(:,2) < 4), 1) = 1;

% Complete les series de correction PAS TRES BON
if COR_mediane(1,1) ~= tsg.TIME(1)
    COR_mediane = [[tsg.TIME(1) COR_mediane(1,2)];COR_mediane];
end

% Pas tr�s bon : La correction ne doit pas �tre appliqu�e + de 10 jours
% avant ou apr�s le dernier echantillon
% --------------------------------------------------------------------
if COR_mediane(length(COR_mediane),1) ~= tsg.TIME(end)
    COR_mediane = [COR_mediane; ...
                  [tsg.TIME(end) ...
                  COR_mediane(length(COR_mediane),2)]];
end


COR = interp1(COR_mediane(:,1), COR_mediane(:,2), tsg.TIME);

tsg.PSAL_ADJ = tsg.PSAL + COR;   %salinite initiale corrige



% ******************************************************************
%                                   Trac�
% ******************************************************************    
figure
plot( tsg.TIME, tsg.PSAL, 'k' );
hold on
plot( tsg.TIME, psal_smooth, 'r' );

figure
hist(nval, 1:13)

figure
plot( tsg.TIME, psal_smooth, 'k' );
hold on
ind = find( sample.PSAL_QC == 1 ); 
plot( sample.TIME(ind), sample.PSAL(ind), '.b' );

ind = find( sample.PSAL_QC == 0 ); 
if ~isempty(ind)
    plot( sample.TIME(ind), sample.PSAL(ind), '.r' );
end

figure
plot( sample.TIME, sample.PSAL_DIF, '.k' );

figure
plot( tsg.TIME, tsg.PSAL, 'k' );
hold on
plot( tsg.TIME, tsg.PSAL_ADJ, 'r' );
plot( sample.TIME, sample.PSAL, '.b', 'Markersize',  6 );