function [smooth] = tsg_average(hTsgGUI, PARA, iTsg)
% Perform the average of a time series at the position of a WS sample
%
% The average is made for a period equal to 'tsg.cst.TSG_DT_SMOOTH'
% Data exceeding the average over that period by 'tsg.cst.TSG_STDMAX'
% are not taken into account.
%
% Input
% hTsgGUI ............ Handle to the main user interface
%
% No computation : NaN
%

% Get the tsg structure from the application
% ------------------------------------------
tsg = getappdata( hTsgGUI, 'tsg_data');

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

% Memory allocation - nval only used for debug
% --------------------------------------------
% smooth = zeros( size(tsg.(PARA)) );
% nval   = zeros( size(tsg.(PARA)) );

% Loop over the tsg.SSPS time series
% -----------------------------------
%h = waitbar(0,'Please wait. Compute a smooth time series ....');
%iEnd = length(tsg.(PARA));
%for i = 1:iEnd
  
  % Display a wait bar
  % ------------------
%  waitbar(i/iEnd);
  
  % Select the param data over 'tsg.cst.TSG_DT_SMOOTH' time interval
  % ind1 : indices of tsg.PARA in the 'tsg.cst.TSG_DT_SMOOTH' time interval
  % ind2 : indices of tsg.PARA not rejected by the S.D. test
  % -----------------------------------------------------------------------
  ind1 = find( tsg.DAYD >= tsg.DAYD(iTsg) - tsg.cst.TSG_DT_SMOOTH/2 & ...
               tsg.DAYD <= tsg.DAYD(iTsg) + tsg.cst.TSG_DT_SMOOTH/2 & ...
               tsg.([PARA '_QC'])(iTsg) <= PROBABLY_GOOD);

	ind2 = ind1;
  if ~isempty(ind2)
    
    currentStd  = Inf;
    previousStd = 0;
    
    % Compare Standard Deviation to the MAX acceptable STD
    % ----------------------------------------------------
    while currentStd > tsg.([PARA '_STDMAX']) && currentStd ~= previousStd
      
      previousStd = currentStd;
      
      % Standard deviation and average over timeInterval
      % ------------------------------------------------
			currentStd = nanstd(  tsg.(PARA)(ind2) );
      meanParam  = nanmean( tsg.(PARA)(ind2) );
            
			% Indices of 'good' values of Param
      % ---------------------------------
      ind2 = ind1( tsg.(PARA)(ind1) >= meanParam - currentStd & ...
                   tsg.(PARA)(ind1) <= meanParam + currentStd );
    end
    
    smooth = nanmean( tsg.(PARA)(ind2) );
    
  else
    smooth = NaN;
  end
  
  % nval only used for debug
  % ------------------------
  % nval( i ) = length( ind2 );
  
%end 
%close(h)

% Transfer the smooth timeseries to the TSG structure
% nval only used for debug
% ---------------------------------------------------
%tsg.ssps_smooth      = smooth;
%tsg.ssps_smooth_nval = nval;

% Update the tsg structure in the application
% --------------------------------------------
%setappdata( hTsgGUI, 'tsg_data', tsg);