function [Ip,Is] = dfft2(yo,v1,v2,zpt)
% Disentangling via FFT.
% This version for PC, for any number of spectra.
% Modifications for filter function, light ratio correction.
% To use:  [Ip,Is] = dfft2(y,v1,v2,zpt);
%
% by D. Holmgren, Apr 28 95 holmgren@phobos.astro.uwo.ca
% bugs corrected May 6 95.
[nl,N] = size(yo); norm = 1/nl;
% set systemic velocity to zero...
v1 = v1 - zpt;
v2 = v2 - zpt;
% FFT each spectrum...
Y = zeros(nl,N);
for k = 1:N
 Y(:,k) = fft(yo(:,k)) .* norm;
end
% set up component spectrum arrays and set I1(0) = I2(0) = 0.
I1 = zeros(nl,1); I2 = I1;
% set up frequency array...
% y = linspace(0.,0.5,nl)';
% y = linspace(0.,1.,nl)';
df = 1./nl; y = df .* (0:1:nl-1)';
y = (2*pi) .* y;
% set up filter function...
f = hanning(nl); f = fftshift(f);
% f = ones(nl,1);
% solve for Fourier coefficients I1(k), I2(k) of component spectra...
% Only go as far as nl-1 because PC can crash at y(nl).
for k = 2:nl
 a12 = 0; a21 = 0; b1 = 0; b2 = 0;
% construct matrix elements from RVs at each phase...
 for l = 1:N
  arg = y(k)*(v1(l)-v2(l));
  a12 = a12 + f(k)*(cos(arg) + i*sin(arg));
  arg = y(k)*(v2(l)-v1(l));
  a21 = a21 + f(k)*(cos(arg) + i*sin(arg));
  arg = y(k)*v1(l);
  b1 = b1 + f(k)*Y(k,l)*(cos(arg) + i*sin(arg));
  arg = y(k)*v2(l);
  b2 = b2 + f(k)*Y(k,l)*(cos(arg) + i*sin(arg));
 end
% set up and solve system...
 a = [N*f(k) a12; a21 N*f(k)];
 b = [b1; b2];
if real(det(a)) == 0,
 I1(k) = 0;
 I2(k) = 0;
else
 x = a\b; I1(k) = x(1); I2(k) = x(2);
end
end
% reconstruct component spectra...
I1 = ifft(I1); I2 = ifft(I2);
% correct normalization...
Ip = real(I1) .* nl; Is = real(I2) .* nl;
% plot results...
w = (0:1:nl-1)';
plot(w,Ip,w,Is)
xlabel('Wavelength')
ylabel('Residual Flux')
title('Reconstructed component spectra')
shg