function reconstructed=lucy(observation,profile,niterations,low_filter_flag,damping,sigma,params)
%Spectral deconvolution using a damped Lucy-Richardson method with suppressed
%high-frequency modification.
%
%RECONSTRUCTED=LUCY(OBSERVATION,PROFILE,NITERATIONS)                        w/o damping nor filtering
%RECONSTRUCTED=LUCY(OBSERVATION,PROFILE,NITERATIONS,0)                      same as above
%RECONSTRUCTED=LUCY(OBSERVATION,PROFILE,NITERATIONS,1)                      filtering w/o damping
%RECONSTRUCTED=LUCY(OBSERVATION,PROFILE,NITERATIONS,0,DAMPING,SIGMA)        damping w/o filtering
%RECONSTRUCTED=LUCY(OBSERVATION,PROFILE,NITERATIONS,1,DAMPING,SIGMA)        both filtering and damping, interactive
%RECONSTRUCTED=LUCY(OBSERVATION,PROFILE,NITERATIONS,1,DAMPING,SIGMA,PARAMS) both filtering and damping, non-interactive
%RECONSTRUCTED=LUCY(OBSERVATION,PROFILE,NITERATIONS,1,0,0,PARAMS)           filtering w/o damping, non-interactive
%
% Returns the deconvolution of OBSERVATION using a RICHARDSON-LUCY method with either
% damping and/or low-pass filtering.
%
% OBSERVATION is a vector with the observed distribution
% PROFILE is a probability distribution of the same length
%         as observation.  The profile must be normalized.
% NITERATIONS is the number of iterations to go through.
% DAMPING is the damping factor.  Typically, it should be an EVEN integer 4 or
%         greater.
% SIGMA is a vector of the same length as OBSERVATION giving the noise levels
%       of the observed values at the various pixels. NOTE: this is not the
%       signal to noise.
% PARAMS is a vector of length 3 which contains the parameters to pass to
%        OPT_FILTER if non-interactive mode is requested.  See OPT_FILTER
%        for details about these parameters. 
% RECONSTRUCTED is the output and represents the deconvolved
%               or extracted original distribution.
%
% The routine, generally, solves the integral (sum, in discrete
% mathematics):
% observation(x)=sum_over_y(reconstructed(y)*profile(x-y))
% for the distribution, reconstructed.
%
% Requires OPT_FILT.M
%
% See L.B.Lucy, Astronomical Journal, V 79, p. 745
%  or http://phobos.astro.uwo.ca/~etittley/projects/lucy.ps

if (nargin<3 | nargin>7) error('Incorrect number of parameters passed to LUCY'), end
if nargin==3 low_filter_flag=0; end
dampflag=0;
if nargin>4
 if damping~=0 
  dampflag=1;
  sigma=sigma/2; % Empirically, this helps a great deal
 end
end


param_flag=0;
if nargin==7
 if length(params)==3
  param_flag=1;
  amp=params(1);
  width=params(2);
  slope=params(3);
 else error('Parameter PARAMS not a vector of length 3 in call to LUCY')
 end
end

profile=profile(:);

% Check to see that observation and profile are the same length
n=length(observation);
m=length(profile);
if m~=n error('Vectors observation and profile passed to lucy must be the same length'), end

% Find the maximum of the profile.  Note, this assumes that the
% maximum actually falls in the center of a pixel!  If it doesn't,
% a shift equal to the deviation will occur.  However, the
% maximum can occur at _any_ pixel.
profile=[profile;profile;profile];
mid=find(max(profile)==profile);
midm=mid(2);
profile_shift=profile(midm:midm+n-1);
profile_fft=fft(profile_shift);
Q=profile(midm-[0:n-1]);
Q_fft=fft(Q);

% A vector the same length as the input vectors, to be used in
% summations.
% template=[1:n];

if low_filter_flag
% Get the filter to be used to suppress high-frequency modifications
 if param_flag [clean,filter]=opt_filt(observation,amp,width,slope);
 else [clean,filter]=opt_filt(observation);
 end
end

reconstructed=observation;	% Initial guess

for iteration=1:niterations	% START OF MAIN LOOP

% Calculate what would be observed if the calculated distribution
% were convolved with the profile. (eq. 14 in the ref.)
% for k=1:n
%  implied_observation(k)=sum(reconstructed.*profile(midm+k-template));
% end
 implied_observation=real(ifft(fft(reconstructed).*profile_fft)); %Fourier Convolution for the above

% Calculate the next guess to the distribution
% (eq. 15 in the ref.)  Note: if implied_observation==observation,
% then the sum is only of the profile, which is unity, so no
% correction is made.
 if ~dampflag factor=observation./implied_observation;  %Undamped
 else factor=1+(observation-implied_observation)./implied_observation.*exp(-(sigma./(observation-implied_observation)).^damping); %Damped
 end

% for k=1:n
%  correction(k)=sum(profile(midm+template-k).*factor);
% end
 correction=real(ifft(fft(factor).*Q_fft)); %Fourier Convolution for above

 if low_filter_flag
% Filter the correction factor
  correction_tmp=real(ifft(fft(correction).*filter'));
  correction=correction_tmp+mean(correction-correction_tmp);
 end

 reconstructed=reconstructed.*correction;

end				%END OF MAIN LOOP