corTsgGradient.m

function [error] = corTsgGradient(hMainFig, PARA, dateMin, dateMax)
% Correct the TSG salinity time series with the Water sample.
% The correction is a linear interpolation of the sample/tsg difference 
% between 2 dates as a function of salinity or temperature (NOT time)
% Only adviced when the sample/tsg difference shifts suddenly when
% crossing a strong salinity gradient, to avoid a correction discontinuity 
% 
% Input
% hMainFig ..... Handle to the main GUI
% PARA ..........Cell array
%                   PARA{1} contains the characters (SSP, SSJT, SSTP)
%                   PARA{2} contains either the cahracters (SSPS, SSJT, SSTP)
%                           or (SSPS_CAL, SSJT_CAL, SSTP_CAL)
% 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
%
% 2009/12/01 G. Alory
% This method is adapted from G. Reverdin, and used when Nuka Arctica 
% crosses the salinity front at the southern tip of Groenland
%
% $Id$

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

% -------------------------------------------------------------------------
% 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 = keys(tsg.qc.hash);

% TODO: define size of keptCode
% -----------------------------
%keptCode = zeros(numel(qc_list), 1);

% iterate (loop) on each key store inside hastable
% ------------------------------------------------
keptCode = [];

nKeptCode = 0;
for key = qc_list

  % get handle of checkbox
  % ----------------------
  hCb = findobj(hMainFig, 'tag', ['TAG_CHECK_CORRECTION_' char(key)]);
    
  if get( hCb, 'value' )
    nKeptCode = nKeptCode + 1;
    keptCode(nKeptCode) = tsg.qc.hash.(key).code;
  end
end

% Get PROBABLY_GOOD, PROBABLY_BAD and VALUE_CHANGED codes
% -------------------------------------------------------
PROBABLY_GOOD = tsg.qc.hash.PROBABLY_GOOD.code;
PROBABLY_BAD  = tsg.qc.hash.PROBABLY_BAD.code;
VALUE_CHANGED = tsg.qc.hash.VALUE_CHANGED.code;

% Intialisation
% 01/09/2009 : intialisation to NaN for real and 0 for byte (QC)
% BE CAREFUL:
% netcdf toolbox failed with assertion when we write NaN to ncbyte variable
% -------------------------------------------------------------------------
if isempty( tsg.([PARA{1} '_ADJUSTED_ERROR']) )  
  tsg.([PARA{1} '_ADJUSTED'])       = NaN*ones(size(tsg.(PARA{1})));
  tsg.([PARA{1} '_ADJUSTED_QC'])    = zeros(size(tsg.([PARA{1} '_QC'])));
  tsg.([PARA{1} '_ADJUSTED_ERROR']) = NaN*ones(size(tsg.(PARA{1}))); 
end

if dateMax > dateMin

  % Find samples within TIME_WINDOWS with Good, probably Good, QC
  % -------------------------------------------------------------
  ind = find( tsg.DAYD_EXT    >= dateMin &  tsg.DAYD_EXT    <= dateMax &...
              tsg.([SAMPLE '_EXT_QC']) <= PROBABLY_GOOD);

  if ~isempty(ind)

    % detect NaN in sample.SSPS_DIF due to bad QC code for tsg.SSPS
    % -------------------------------------------------------------
    ind2 = find(~isnan(tsg.EXT_DIF(ind)));

    if ~isempty(ind2) && length(ind2) >= 2

      % Locate front with big shift in TSG/sample difference  
      % -------------------------------------------------------------
%       tsg_diff=diff(tsg.EXT_SMOOTH(ind(ind2)));
%       corr_diff=diff(tsg.EXT_DIF(ind(ind2)));
    
      % The correction is applied to the TSG between dateMin and dateMax using
      % a linear interpolation only on measurements with keptCode Quality Codes
      % ------------------------------------------------------------------------
      dtTsgQCkept=[];
      for icode = 1 : length( keptCode )
        dtTsg = find( tsg.DAYD    >= tsg.DAYD_EXT(ind(ind2(1)))  & tsg.DAYD <= tsg.DAYD_EXT(ind(ind2(end))) &...
                      tsg.([PARA{1} '_QC']) == keptCode( icode ));

        if ~isempty( dtTsg )

        dtTsgQCkept=[dtTsgQCkept; dtTsg];
                
        % Compute the correction + the error
        % ----------------------------------
        tsg_diff=diff(tsg.EXT_SMOOTH(ind(ind2)));

        for iint=1:(length(ind2)-1)
        % in each interval between 2 samples, the correction is a linear interpolation
        % of the sample/tsg difference as a function of Tsg SSPS or SSTP,
        % the error an interpolation of 0.5*abs(sample/tsg diff.)
        
            dtTsg2 = find( tsg.DAYD(dtTsg) >= tsg.DAYD_EXT(ind(ind2(iint)))  &...
                tsg.DAYD(dtTsg) <= tsg.DAYD_EXT(ind(ind2(iint+1))) &...
                tsg.([PARA{1} '_QC'])(dtTsg) == keptCode( icode ));
            
            if ~isempty(dtTsg2)
                
                tsg_rel=(tsg.(PARA{2})(dtTsg(dtTsg2))-tsg.EXT_SMOOTH(ind(ind2(iint))))/tsg_diff(iint);
                tsg_rel(tsg_rel<0)=0;
                tsg_rel(tsg_rel>1)=1;
                tsg.([PARA{1} '_ADJUSTED'])(dtTsg(dtTsg2)) =...
                    (1-tsg_rel)*tsg.EXT_DIF(ind(ind2(iint))) + tsg_rel*tsg.EXT_DIF(ind(ind2(iint+1)));
                tsg.([PARA{1} '_ADJUSTED_ERROR'])(dtTsg(dtTsg2)) =...
                    (1-tsg_rel)*abs(tsg.EXT_DIF(ind(ind2(iint)))/2) + tsg_rel*abs(tsg.EXT_DIF(ind(ind2(iint+1)))/2);

                % Compute the corrected value : original value + correction
                % ---------------------------------------------------------
                tsg.([PARA{1} '_ADJUSTED'])(dtTsg(dtTsg2)) =...
                    tsg.(PARA{2})(dtTsg(dtTsg2)) + tsg.([PARA{1} '_ADJUSTED'])(dtTsg(dtTsg2));
            end
        end
        % Transfer the QC
        % ---------------
        tsg.([PARA{1} '_ADJUSTED_QC'])(dtTsg) = tsg.([PARA{1} '_QC'])(dtTsg);
        end
      end

      % The error minimum cannot be lower than the sensor accuracy
      % ----------------------------------------------------------
      accuracy= tsg_accuracy(hMainFig, PARA{1},dtTsgQCkept );
      error_too_small=find(abs(tsg.([PARA{1} '_ADJUSTED_ERROR'])(dtTsgQCkept)) < accuracy);

      if ~isempty(error_too_small)
          if length(ind2) > 2
              tsg.([PARA{1} '_ADJUSTED_ERROR'])(dtTsgQCkept(error_too_small))=accuracy(error_too_small);
          else
              tsg.([PARA{1} '_ADJUSTED_ERROR'])(dtTsgQCkept(error_too_small))=-accuracy(error_too_small);
          end
      end

      % Case with 1 sample only: we can't apply this correction
      % --------------------------------------------------------------------------
    else

      msgbox( {'Function corTsgGradient'; ' '; ...
               'At least 2 samples are needed';...
               'The soft cannot make the interpolation'},...
              'Correction Method', 'warn');
      error = -1;
      
      return

    end

  end

  % Update tsg application data
  % ---------------------------
  setappdata( hMainFig, 'tsg_data', tsg);

  % everything OK
  % -------------
  error = 1;

else

  % DateMax <= DateMin
  % ------------------
  error = -1;

end