function flattened=contfind(spec,mask_delta,mask_width,order)

% Fit a continnum to the spectrum, SPEC
%
% Syntax: flat=contfind(spec);
%              Returns the contiuum-flattened spectrum determined
%              by the method for MASK_FLAG=0.
% Syntax: flat=contfind(spec,mask_delta,mask_width,order)
%              Returns the contiuum-flattened spectrum determined
%              by the method for MASK_FLAG=1.
%
% The choice of methods is given by the value of MASK_FLAG.
%  MASK_FLAG=0: the continuum is first assumed to be represented by
%               the peaks in the spectrum, which are found by
%               performing a LINDFIND() on the inverted spectrum.
%               A linear fit is made to these points, by which the 
%               spectrum is divided.
%               After this, the continuum is approximated by the top
%               30 points in the spectrum, to which a linear fit is
%               made leading to the final division to output a
%               continuum-flattened spectrum.
%  MASK_FLAG=1  A mask is created to determine the location of
%               continuum, which is used to estimate the regions of 
%               continuum to which a polynomial of order ORDER is fit.
%               The output is the spectrum divided by this fit.
%
% Note: Requires the functions, lindfind.m, mask2.m, continuum.m, 
% and polyrecon.m to be on hand
%
% Eric Tittley (95 07)

 if nargin<1 error('Too few arguments in call to contfind'); end
 mask_flag=0;
 if nargin>1 mask_flag=1; end
 if (nargin>1) &  (nargin<4) error('Mask parameters set incorrectly in call to contfind'), end
 if nargin>4 error('Too many arguments in call to contfind'); end

% The range over which to do the fit.  This should be an input 
% parameter, but there are too many, as it is.
if length(spec)<1700 range=[1,length(spec)];
else range=[200,1600];
end

% The range of sorted amplitudes to use for the final continuum fit.
% This should be an input parameter, but there are too many, as it is.
 sort_range=[3:300];
 
 wave=[1:length(spec)]';
 spec=spec(:);

if mask_flag
 mask=mask2(spec,mask_delta,mask_width,range);
 a=find(mask);
 pout=polyfit(wave(a),spec(a),order);
else
 level=mean(-spec);
 a=linefind(wave,-spec,level);
 pout=polyfit(wave(a),spec(a),1);
end
 continuum=polyval(pout,wave);
 spec2=spec./continuum;
 [dummy,a]=sort(spec2(range(1):range(2)));
 a=sort(a( length(a)-sort_range(length(sort_range)) : length(a)-sort_range(1) ))+range(1)-1;
 pout=polyfit(wave(a),spec2(a),1);
 continuum=polyval(pout,wave);
 flattened=spec2./continuum;