function mask=mask2(vector,delta,width,range_limits);
% Continuum mask creation
%
% Syntax:
% mask=mask2(vector,delta,width)
% mask=mask2(vector,delta,width,range_limits)  if you wish to work
%                      only within a specific range.
%
% This function creates a mask for a spectrum to isolate regions of
% continuum.  The routine looks for regions with local standard deviations
% less than DELTA of the std. dev. of the whole spectrum.  This leads to a
% true-false mask (a vector of 1's and 0's).  Then it box-car
% averages this true-false mask with a box of width, WIDTH, which is turned
% into a mask.
%
% VECTOR is the input spectrum for which the mask is to be created.
% DELTA is the factor less than the global standard deviation that the
%   local standard deviations must be less than to be considered continuum.
% WIDTH is the size of the box to be used to smooth out the distibution.
%
% Note, the routine isn't particularily robust and requires a lot of
% guess-work as to the parameters above.  DELTA should be around .2
% and WIDTH should be around 5.  WIDTH must be a positive integer.
%
% Eric Tittley (95 01 30) etittley@phobos.astro.uwo.ca

 range_flag=1;
 if (nargin<3) error('Too few arguments in call to MASK'); end
 if (nargin==3) range_flag=0; end
 if (nargin>4) error('Too many arguments in call to MASK'); end
 
 n=length(vector);
 vector=vector(:);

 if(range_flag)
  range=[range_limits(1):range_limits(2)];
  large_deviation=std(vector(range));
 else large_deviation=std(vector);
 end

 local_deviation=large_deviation*(1&vector)';
 array=[1:21]'*vector(11:n-10)'; % predefined for speed
 for k=-10:10
  array(k+11,:)=vector(11-k:n-10-k)';
 end
 local_deviation(11:n-10)=std(array);  %row vector

% The above does this, just much faster
% for k=11:n-10
%  local_deviation(k)=std(vector(k-10:k+10));
% end
 
 mask=local_deviation<(delta*large_deviation); %row vector
 if(width>0)
  array=[1:2*width+1]'*mask(width+1:n-width);
  for k=-width:width
   array(k+width+1,:)=mask(width+1-k:n-width-k);
  end
  mask(width+1:n-width)=round(mean(array)); %row vector
 end

% The above does this, just much faster
% average=0*mask;
% for k=width+1:n-width average(k)=round(mean(mask(k-width:k+width))); end
% mask=average;

% Remove any obvious discrepancies.  First, flatten the spectrum. Then
% find all regions more than 3 sigmas away.
 good=0;
 x=find(mask)';
 y=vector(x);

while(~good);
 good=1;
 pout=polyfit(x,y,1);
 residuals=polyval(pout,x)-y;
 crap=residuals>(2*std(residuals));
 if(~isempty(find(crap)))
  mask(x(find(crap)))=0*[1:length(find(crap))];
  good=0;
 end
 x=find(mask)';
 y=vector(x);
% plot(x,y);
% good=input('Good enough? (0 for no, 1 for yes, -1 for exit)');
% if(good<0) error('User break from MASK2') end
end

% Set regions outside the range to zero.
 if(range_flag)
  if(range_limits(1)<=1) range_limits(1)=2; end
  if(range_limits(2)>=length(mask)) range_limits(2)=length(mask)-1; end
  mask(1:range_limits(1)-1)=0*[1:range_limits(1)-1];
  mask(range_limits(2)+1:length(mask))=0*[range_limits(2)+1:length(mask)];
 end
 mask=mask';