function [demeaned,pattern] = vdemean(data, varargin)
% [demeaned, pattern] = vdemean(data, varargin);
% data is Nbins X Ntimes 
%
% varargin       value       default              what
% 'meandir'      'v','h'       'v'          'v':  remove a vertical reference 
%                                                     layer from each profile
%                                           'h': remove a time-varying signal
%                                                     from each bin
% 'meantype'   'refavg','mean'  'refavg'      use this method to calculate
%                                                  value to be removed
% 'pattern'     vector     (no action)        take out this specific pattern
%                                                   (dimensions must match
%                                                     'meandir')
% 'refbins'      [indices]  (all)               use these indices for pattern
%
% eg:
%   demeaned = vdemean(data, 'meandir', 'h') % refavg timeseries removed 
%   demeaned = vdemean(data, 'refbins', [3:10], 'meantype', 'mean');
%                                                  %old fashioned ref layer
%
 

% surely there's one of these somewhere already?
% JH 2001/06/19


cfg.pattern = [];
cfg.refbins = [];
cfg.meantype = 'refavg';
cfg.meandir = 'v';


[nrows, ncols] = size(data);

if (nargin == 0)
   help vdemean
   return
end

cfg = fillstruct(cfg, varargin,0);

if (~ismember(cfg.meandir, {'h','v'}))
   help vdemean
   return
end



if strcmp(cfg.meandir, 'h') %take out a timeseries
   if isempty(cfg.pattern)
      if isempty(cfg.refbins)
         %use the whole thing
         cfg.refbins = 1:length(data(:,1));
      end
      % now we have ref bins
      if strcmp(cfg.meantype', 'mean')
         pattern = mstdgap(data(cfg.refbins,:));
      else
         if exist('refavg') ~= 2
            addpath(fullfile(fileparts(which('adcppath')), 'codas3'))
         end
         [vp, pattern] = refavg(data, cfg.refbins, 2);
      end
   else
      pattern = cfg.pattern;
   end
   if length(pattern) ~= length(data(1,:))
      fprintf('pattern has wrong dimensions. should match ''meandir''\n')
      cfg
      demeaned = [];
      return
   end
   demeaned = data - repmat(pattern, nrows, 1);
else % take out a profile
   if isempty(cfg.pattern)
      if isempty(cfg.refbins)
         %use the whole thing
         cfg.refbins = 1:length(data(:,1));
      end
      % now we have ref bins
      if strcmp(cfg.meantype', 'mean')
         pattern = mstdgap(data(cfg.refbins,:));
      else
         if exist('refavg') ~= 2
            addpath(fullfile(fileparts(which('adcppath')), 'codas3'))
         end
         [pattern] = refavg(data, cfg.refbins, 2);
      end
   else
      pattern = cfg.pattern(:);
   end
   if length(pattern) ~= length(data(:,1))
      fprintf('pattern has wrong dimensions. should match ''meandir''\n')
      cfg
      demeaned = [];
      return
   end
   demeaned = data - repmat(pattern, 1, ncols);
end


