function [y2,x,ex]=fitsine(t,y,f,n,t_int,x_ap);

keybPause=0;    % plot results and pause after fitting

% function [y2,x,ex]=fitsine(t,y,f,n,t_int,x_ap)

% Decompose signal y(t) into a mean and sinusoids 
% of known frequency, unknown amplitude and phase.
% f is a 1 by n vector of frequencies,
% in units of cycles per XXX (units of vector t).
%  n: amplitude of "noise" (anything besides mean + sinusoid)
%   in units of y.  If not specified, this routine assumes
%   n is std(y).
% t_int is the integral time scale, used to determine
%   the spacing of "independent" observations.
% x_ap, if provided, set the prior values of the
%   amplitudes in vector x.
%
% Output: y2=[cos(2*pi*t*f),sin(2*pi*t*f),ones(size(t))]*x
%   x is the amplitudes, ex contains the estimated 
%   errors for amplitudes x.

%if exist('x_ap') keypause=1; end;

if ~exist('n')
    n=std(y);
end;
if ~exist('t_int')
    t_int=0;
end;
if ~exist('x_ap')
    x_ap=[];
end;
if isnan(sum(x_ap))
    x_ap=[];
end;

if (size(t,2)>size(t,1))
   t=t'; y=y';
end;

yorig=y;torig=t;
t=t+1979;torig=torig+1979;

if (length(t)>1 & t_int>0)  % subsampling code
    dt2=1/365.24;
    t2=min(t):dt2:max(t);
    y2=NaN*ones(size(t2));
    num2=y2;
    for i=1:length(t2);
        q=find(t>=t2(i)-dt2/2 & t<=t2(i)+dt2/2);
        if (length(q)) 
            y2(i)=mean(y(q)); 
            num2(i)=length(q);
        end;
    end;
    q=find(isfinite(y2));
    t=t2(q)';y=y2(q)';num=num2(q)';
end;

% first guess of mean:  least-squares adjust from here
if length(x_ap)
    y_ap=[cos(2*pi*t*f),sin(2*pi*t*f),ones(size(t))]*x_ap;
else
    ybar=mean(y);
    y_ap=ybar*ones(size(t));
end;
y=y-y_ap;

A=([cos(2*pi*t*f),sin(2*pi*t*f),ones(size(t))]);
 
% tapered, weighted least squares.
% Prior amplitudes:
%   sinusoidal components: range(y)/2
%       (evenly-distributed between the frequencies)
%   adjustment to mean: +-range(y)/2
Rxx=diag( [range(y)/2/length(f)*ones(1,2*size(f,2)) ...
        range(y)/2].^2);

% calculate an Rnn with Lagrangian integral 
% time scale of t_int (units of t)
if (t_int>0)
  dT=abs(ones(length(y),1)*t' -t*ones(1,length(y)));
  td=sqrt(pi/2)*t_int/exp(-.5);
  te=2*sqrt(2)/pi*td;
  Rnn=(n^2)*cos(pi*dT/(2*td)).*exp(-(dT/te).^2);
  whitenoise_to_eddy=.1;
  Rnn=(1-whitenoise_to_eddy)*Rnn+ ...
    whitenoise_to_eddy*diag(n^2*ones(length(y),1));

  C=Rxx*A'/(A*Rxx*A'+Rnn);    % use full structure of Rnn
else
  % quick code, appropriate for diagonal Rnn=n^2*eye
  C=inv(inv(Rxx)+A'*(1/n^2)*A)*A'*(1/n^2);
end;

x=C*y;
ex=sqrt(diag(Rxx-C*A*Rxx));

y2=A*x + y_ap;

xrat=abs(x)./sqrt(diag(Rxx));
if (length(x_ap))  x=x+x_ap;
else
    x(length(x))=x(length(x))+ybar;
end;

ybar=x(length(x));

if (keybPause) % plot results and keyboard
    subplot(211)
    tt=(min(t):(max(t)-min(t))/100:max(t))';
    yy=[cos(2*pi*tt*f),sin(2*pi*tt*f),ones(size(tt))]*x;
    plot(torig,yorig,'.',tt,yy);
    tt2=(0:.01:1)';
    yy2=[cos(2*pi*tt2*f),sin(2*pi*tt2*f),ones(size(tt2))]*x;
    subplot(212)
    plot(torig-floor(torig),yorig,'.',tt2,yy2)
    fprintf('mean: %5.2f +- %5.2f\n',x(length(x)),ex(length(ex)));
    fprintf('sine amp: %5.2f +- %5.2f\n', ...
        sqrt(sum(x(1:2:length(x)-1).^2)), ...
        sqrt(sum(ex(1:2:length(ex)-1).^2)) );
    keyboard;
end;
if (max(xrat)>2 & ~length(x_ap))    
    y2=NaN*ones(size(y));
    x(1:length(x)-1)=NaN;
    ex(1:length(ex)-1)=NaN;
    x(length(x))=mean(yorig);
end;
