function d1 = decpro(wi,ip,prof)
% Function to do Fourier deconvolution.
% Applied to stellar line profile given an instrumental
% profile.  To use:
%  d1 = decpro(wi,ip,prof);
%
% This routine works. 20/2/95
% Sun version: filtering menu is a GUI!
%
% by D. Holmgren, holmgren@phobos.astro.uwo.ca
ni = length(ip); ns = length(prof);
% w = hanning(ns);
n2 = ns/2;
%  Apply Hanning window to residual flux profile...
% prof = w .* (1 - prof);
% Modified to use 10% cosine bell...
prof = 1 - prof;
for k = 1:10
 indx = ns + 1 - k;
 w = 0.5 * ( 1 - cos(pi*k/10));
 prof(k) = prof(k)*w;
 prof(indx) = prof(indx)*w;
end
% Set up FFT arrays...
nfft = input('FFT length: ');
y = zeros(nfft,1); z = y;
y(1:n2) = prof(n2+1:ns);
y(nfft-n2+1:nfft) = prof(1:n2);
% Get wavelength intervals...
dlami = input('     Dlam for IP: ');
dlams = input('Dlam for Stellar: ');
% Put instrumental profile onto stellar wavelength scale...
npts = 0;
for k = 1:ns
 wl = (k-1)*dlams;
 if wl >= wi(ni), break, end
 z(k) = spline(wi,ip,wl);
 npts = npts + 1;
end
z(nfft-npts+2:nfft) = flipud(z(2:npts));
% plot this as a check...
plot(fftshift(z))
xlabel('Delta lambda (A)')
ylabel('Profile')
title('Instrumental Profile on Stellar Scale')
pause
% normalise interpolated instrumental profile...
si = sum(z); z = z ./ si;
% Do the FFTs...
C = fft(z,nfft);
P = fft(y,nfft);
% ...and the deconvolution...
D1 = P ./ C;
% Allow for Fourier filtering...
% Note that this is standard optimal filtering.
% See DFG's book about this.
ny = menu('Fourier filtering','Apply filter','No filtering');
if ny == 1,
 bta = input('Enter S/N: '); bta = 1/bta;
 fwhm = input('Enter FWHM of IP: ');
 alpha = fwhm/nfft;
 disp(alpha); pause(1)
 phi = filt(bta,alpha,nfft/2);
 f = [phi; flipud(phi)];
 D1 = D1 .* f;
end
% d1 = deconvolved line profile (converted to real and
%      unwrapped).
d1 = ifft(D1);
d1 = fftshift(real(d1));
% Plot original profile and deconvolved...
x = (0:1:nfft-1)';
plot(x,fftshift(y),'o',x,d1)
xlabel('Wavelength')
ylabel('Residual Flux')
title('Deconvolved Profile')
% Don't need shg for Sun Matlab...
% shg