function coefs=tidefit(uu,tt,tstart,tstop,range)
%function coefs=tidefit(uu,tt,tstart,tstop,range)
%        U and T are the matrices containing the variable to be fit in U, and
%        the associated times in T.  Each row of U should be a time series.
%        TSTART and TSTOP indicate the time range to be plotted.
%        RANGE is an optional argument to specify the number of bins to plot.
%        It should be a vector, for example, RANGE=[1:10] will fit bins
%        through 10.  If no range is specified, all the bins are fit.
%
period=[12.42 6.21 4.14]        % can add 24 if interested
if nargin==4
  [mu,nu]=size(uu);
  range=[1:mu];
end
tinc=(tstop-tstart)/50;
tp=[tstart:tinc:tstop];
tp=tp(:);
ntp=length(tp);
freq=(2.*pi)*ones(period)./period;
nfreq=length(freq);
for nn=range
   u=uu(nn,:);
   if length(tt(:))==length(uu(:))
         t=tt(nn,:);
   else
         t=tt;
   end
   ind=find(u~=NaN);
   nind=length(ind);
   if nind>=9
      u=u(ind);
      t=t(ind);
      u=u(:);
      t=t(:);
%------ set up A -------
       A=zeros(nind,nfreq*2+1);
       A(:,1)=ones(nind,1);
       for j=[1:nfreq]
           A(:,2*j)=cos(freq(j)*t);
           A(:,2*j+1)=sin(freq(j)*t);
       end
%-------solve [A coeff = u] -----------------
       coef=A\u;
%-------generate solution components-----
       up=zeros(ntp,nfreq*2+1);
       up(:,1)=coef(1)*ones(ntp,1);
       for j=[1:nfreq]
          up(:,j+1)=coef(j*2)*cos(freq(j)*tp)+coef(j*2+1)*sin(freq(j)*tp);
       end
       up(:,nfreq+2)=sum(up(:,[1:nfreq+1])')';  % sum of all comps
       
       error=std(A*coef-u)/sqrt(nind-length(period)*2+1);
   else
       up=ones(ntp,nfreq+2)*NaN;
       coef=ones(nfreq*2+1,1)*NaN;
   end

   plot(t,u,'o',tp,up,tp,zeros(tp),'w')
   title(sprintf('Least Squares Fit at Bin %g   error= %5.2g cm/s',nn,error));
   xlabel('Time, hours');
   ylabel('cm/second')
   pause
   coefs(:,nn)=coef;
end

