function [ucoef,unew]=tidefit(jd,uu,goplot,period)
%function [ucoef,unew]=tidefit(jd,uu,[goplot],[period])
%  with goplot=1, it plots out the fit
%	if goplot=0, or no argument, it doesn't plot
%  period is optional argument
%   note:  shape of uu(nt,m) where nt is the time dimension and m is the number
%     of separate timeseries
%   jd can be a vector or a matrix, with the same shape as uu
%   ucoef are the coefficients, starting with mean, then cos(f1), then sin(f1), 
%             up to the number of periods.  Hence there are 2f+1 coefficients.

  
if nargin<3; goplot=0; end;
if nargin<4; period=[12.42 6.21 4.14]; end        % can add 24 if interested

[n,m]=size(uu);
[njd,mjd]=size(jd);
if (n~=njd)&(m~=mjd); jd=jd'; end
if n<5; uu=uu';jd=jd'; [n,m]=size(uu); end;
if n<m; disp('note: time should be first dimension of uu'); end

range=m;

tt=jd*24;
freq=(2.*pi)*ones(period/period(1))./period;
nfreq=length(freq);

for nn=1:range
fprintf('%4.0f',nn);
u=uu(:,nn);
   if length(tt(:))==length(uu(:))
         t=tt(:,nn);
   else
         t=tt;
   end
	ind=find(u~=NaN);
   nind=length(ind);

   if nind>length(period)*2+1
      u=u(ind);
      t=t(ind);
      u=u(:);
      t=t(:);
      nt=length(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(nt,nfreq*2+1);
       up(:,1)=coef(1)*ones(nt,1);
       for j=[1:nfreq]
          up(:,j+1)=coef(j*2)*cos(freq(j)*t)+coef(j*2+1)*sin(freq(j)*t);
       end
       up(:,nfreq+2)=sum(up(:,[1:nfreq+1])')';  % sum of all comps
       unew(:,nn)=up(:,nfreq+2);
		 ucoef(:,nn)=coef;
		 error=std(A*coef-u)/sqrt(nind-length(period)*2+1);
       if goplot
          jdp=t/24;
          plot(jdp,u,'oc5',jdp,up,jdp,unew(:,nn),'x',jdp,zeros(jdp),'w')
        	 title(sprintf('Least Squares Fit, series %g   error= %5.2g cm/s',nn,error));
	       xlabel('julian day');
	       ylabel('cm/second')
          pause
       end
else
       unew(:,nn)=ones(t)*nan;
		 ucoef(:,nn)=ones(2*nfreq+1,1)*nan;
end   % if nind> ....
end	%  loop
fprintf('\n')

