corTsgMedian.m

function [error] = corTsgMedian(hMainFig, PARA, dateMin, dateMax)
% 
% Correct the TSG salinity time series with the Water sample.
% Use the median value of TIME_WINDOWS water sample to compute the
% correction. see the documentation
% 
% Input
% hMainFig ..... Handle to the main GUI
% PARA ......... Parameter (SSPS, SSJT, SSTP)
% dateMin ...... the correction is applied between dateMin and date Max
% dateMax ...... the correction is applied between dateMin and date Max
%
% Output
% Error ........  1 everything OK
%       ........ -1 dateMax <= date Min
%
% The error is maximum, equal to 1, if the median is computed with less
% than 4 samples
% The error minimum cannot be lower than 0.01

% Get application data
% --------------------
tsg  = getappdata( hMainFig, 'tsg_data');
SAMPLE = tsg.preference.sample;

% Shorten the variable name
% -------------------------
TIME_WINDOWS = tsg.cst.COR_TIME_WINDOWS;

% -------------------------------------------------------------------------
% Get from the checkbox the QC code on which the correction will be applied
% -------------------------------------------------------------------------

% get list of keys from hashtable tsg.qc.hash, defined inside
% tsg_initialisation.m
% -----------------------------------------------------------
qc_list = get(tsg.qc.hash);

% iterate (loop) on each key store inside hastable
% ------------------------------------------------
keptCode = [];
nKeptCode = 0;
for i=1:numel(qc_list)

  % get key and some values in hashtable
  % ------------------------------------
  key   = qc_list{i};
  
  % get handle of checkbox
  % ----------------------
  hCb = findobj(hMainFig, 'tag', ['TAG_CHECK_CORRECTION_' key]);
    
  if get( hCb, 'value' )
    nKeptCode = nKeptCode + 1;
    keptCode(nKeptCode) = get(tsg.qc.hash, key, 'code');
  end
end

% Get PROBABLY_GOOD, PROBABLY_BAD and VALUE_CHANGED codes
% -------------------------------------------------------
PROBABLY_GOOD = get(tsg.qc.hash, 'PROBABLY_GOOD', 'code');
PROBABLY_BAD  = get(tsg.qc.hash, 'PROBABLY_BAD',  'code');
BAD           = get(tsg.qc.hash, 'BAD',           'code');
VALUE_CHANGED = get(tsg.qc.hash, 'VALUE_CHANGED', 'code');

% intialisation
% -------------
if isempty( tsg.([PARA '_ADJUSTED']) )  
  tsg.([PARA '_ADJUSTED'])       = tsg.(PARA);
  tsg.([PARA '_ADJUSTED_QC'])    = tsg.([PARA '_QC']);
end
% If only calibration have been applied ERROR is empty
% ---------------------------------------------------
if isempty( tsg.([PARA '_ADJUSTED_ERROR']) )  
  tsg.([PARA '_ADJUSTED_ERROR']) = NaN* ones(size(tsg.(PARA))); 
end

% Create a structure with an NaN
% No other solution, as I can't add a structure to an empty one
% -------------------------------------------------------------
cor = struct('DAYD', NaN, 'DIFF', NaN, 'ERROR', NaN, 'NVALUE', NaN);

if dateMax > dateMin
  
  % Find the indices of samples within the time limits.
  % --------------------------------------------------
  indSample = find(tsg.DAYD_EXT >= dateMin & tsg.DAYD_EXT <= dateMax);
  
  indCor = 0;
  for i = 1:length(indSample)

    % Find samples within TIME_WINDOWS with Good, probably Good, probably Bad QC
    % ---------------------------------------------------------------
    ind = find( tsg.DAYD_EXT >= (tsg.DAYD_EXT(indSample(i)) - TIME_WINDOWS/2) &...
                tsg.DAYD_EXT <= (tsg.DAYD_EXT(indSample(i)) + TIME_WINDOWS/2) &...
                tsg.([SAMPLE '_EXT_QC']) <= PROBABLY_GOOD);
                
    if ~isempty(ind)
      
      % detect NaN in sample.SSPS_DIF due to bad QC code in tsg.SSPS
      % ------------------------------------------------------------
      ind2 = find(~isnan(tsg.EXT_DIF(ind)));

      % Compute the median difference and error within TIME_WINDOWS
      % -----------------------------------------------------------
      if ~isempty(ind2)
        if ~isempty(tsg.EXT_DIF(ind(ind2)))

          A = tsg.EXT_DIF(ind(ind2));
          meanA = mean(A);
          stdA  = std(A);

          % Standard deviation test: keep these values
          % ------------------------------------------
          ind3 = find( A >= meanA-3*stdA & A <= meanA+3*stdA);

          B = tsg.EXT_DIF(ind(ind2(ind3)));
          if ~isempty( B )
            indCor = indCor + 1;
            cor.DAYD(indCor)   = tsg.DAYD_EXT((indSample(i)));
            cor.DIFF(indCor)   = median(B);
            cor.ERROR(indCor)  = nanstd(B)/sqrt(length(B));
            cor.NVALUE(indCor) = length(B);
          end

          % Standard deviation test: don't keep these values
          % ------------------------------------------------------
          ind4 = find( A < meanA-3*stdA | A > meanA+3*stdA);
          if ~isempty( ind4 )
            tsg.([SAMPLE '_EXT_QC'])(ind(ind2(ind4))) = BAD;
          end
        end
      end
    end
  end

  % Eliminate the first element if NaN
  % ----------------------------------
  if isnan(cor.DAYD(1))
    cor.DAYD(1)   = [];
    cor.DIFF(1)   = [];
    cor.ERROR(1)  = [];
    cor.NVALUE(1) = [];
  end

  if ~isempty( cor.DAYD )

    % The error is maximum, equal to 1, if the median is computed with less
    % than 4 samples
    % -----------------------------------------------------------------------
    cor.ERROR( cor.NVALUE < 4 ) = 1;
    
    % The error minimum cannot be lower than 0.01
    % -------------------------------------------
    cor.ERROR( cor.ERROR < 0.01 ) = 0.01;

    % The correction is applied between dateMin and dateMax
    % We attribute to dateMin the first correction computed
    % and to dateMax the last one
    %
    % Find the tsg date in the interval dateMin-dateMax
    % -------------------------------------------------
    dtTsg = find(tsg.DAYD >= dateMin  & tsg.DAYD <= dateMax);

    if cor.DAYD(1) ~= dateMin
      cor.DAYD   = [tsg.DAYD(dtTsg(1)) cor.DAYD];
      cor.DIFF   = [cor.DIFF(1)        cor.DIFF];
      cor.ERROR  = [cor.ERROR(1)       cor.ERROR];
      cor.NVALUE = [cor.NVALUE(1)      cor.NVALUE];
    end
    if cor.DAYD(end) ~= dateMax
      cor.DAYD   = [cor.DAYD   tsg.DAYD(dtTsg(end))];
      cor.DIFF   = [cor.DIFF   cor.DIFF(end)];
      cor.ERROR  = [cor.ERROR  cor.ERROR(end)];
      cor.NVALUE = [cor.NVALUE cor.NVALUE(end)];
    end
    
    % Test if the date are strictly increasing otherwise interp1
    % cannot be used
    % ----------------------------------------------------------
    if ~isempty( find(diff(cor.DAYD) == 0) )
      msgbox( {'Function corTsgMedian'; ' '; ...
               'At least 2 samples have the same date';...
               'The soft cannot make the interpolation'},...
              'Correction Method', 'warn');
      error = -1;
      return
    end

    % The correction is applied to the TSG between dateMin and dateMax using
    % a linear interpolation only on measurements with keptCode
    % ----------------------------------------------------------------------
    for icode = 1 : length( keptCode )
      dtTsg = find( tsg.DAYD    >= dateMin  & tsg.DAYD <= dateMax &...
                    tsg.([PARA '_QC']) == keptCode( icode ));

      tsg.([PARA '_ADJUSTED'])(dtTsg)       = tsg.(PARA)(dtTsg) + ...
                                  interp1(cor.DAYD, cor.DIFF, tsg.DAYD(dtTsg));
      tsg.([PARA '_ADJUSTED_ERROR'])(dtTsg) = ...
                                  interp1(cor.DAYD, cor.ERROR, tsg.DAYD(dtTsg));
    end
        
  end

  % Update tsg application data
  % ---------------------------
  setappdata( hMainFig, 'tsg_data', tsg);
  
  % everything OK
  % -------------
  error = 1;

else

  % DateMax <= DateMin
  % ------------------
  msgbox( {'Function corTsgMedian'; ' '; ...
           'Date limits are not correct' },...
           'Correction Method', 'warn', 'modal');

  error = -1;
  
end