Home > bobstuff > tide_fit.m

tide_fit

PURPOSE ^

TIDE_FIT Simple least-squares tidal analysis.

SYNOPSIS ^

function [ucoef,unew]=tide_fit(jd,uu,periods,goplot)

DESCRIPTION ^

 TIDE_FIT Simple least-squares tidal analysis. 

  [ucoef,unew]=tide_fit(jd,uu,periods,goplot)

  Input:  
      uu = matrix of time series, one time series in each column
     jd  = julian day vector  (as in julian.m and gregorian.m) 
 periods = vector of periods to fit (hours) (e.g. [12.42 6.21])  for M2, M4
  goplot = optional argument, which, if supplied will cause a plot of 
           the tidal fit to be generated.  (the actual value you supply
           is arbitrary)
  Output: 
   ucoef = the coefficients, starting with mean, then cos(f1), then sin(f1), 
           up to the number of periods.  Hence there are 2f+1 coefficients.
           To get amplitude and phase, use TIDE_ELL.
   unew =  time series of predicted tide 

 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 Version 1.0 (11/24/1997)   Rich Signell (rsignell@usgs.gov) 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [ucoef,unew]=tide_fit(jd,uu,periods,goplot)
0002 % TIDE_FIT Simple least-squares tidal analysis.
0003 %
0004 %  [ucoef,unew]=tide_fit(jd,uu,periods,goplot)
0005 %
0006 %  Input:
0007 %      uu = matrix of time series, one time series in each column
0008 %     jd  = julian day vector  (as in julian.m and gregorian.m)
0009 % periods = vector of periods to fit (hours) (e.g. [12.42 6.21])  for M2, M4
0010 %  goplot = optional argument, which, if supplied will cause a plot of
0011 %           the tidal fit to be generated.  (the actual value you supply
0012 %           is arbitrary)
0013 %  Output:
0014 %   ucoef = the coefficients, starting with mean, then cos(f1), then sin(f1),
0015 %           up to the number of periods.  Hence there are 2f+1 coefficients.
0016 %           To get amplitude and phase, use TIDE_ELL.
0017 %   unew =  time series of predicted tide
0018 %
0019 %
0020 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0021 % Version 1.0 (11/24/1997)   Rich Signell (rsignell@usgs.gov)
0022 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0023 
0024 % if no periods are supplied, just do m2, m4 and m6
0025 if nargin<3; periods=[12.42 6.21 4.14]; end    %M2, M4, M6
0026 
0027 
0028 [m,n]=size(uu);
0029 
0030 freq=(2.*pi)*ones(size(periods))./periods;
0031 nfreq=length(freq);
0032 
0033 tt=jd*24;
0034 nt=length(tt);
0035 nfreq=length(freq);
0036 
0037 for nn=1:n
0038 fprintf('Series %4.0f. Hit <Return> to continue\n',nn);
0039    u=uu(:,nn);
0040    t=tt;
0041    ind=find(~isnan(u));
0042    nind=length(ind);
0043 
0044    if nind>nfreq*2+1
0045       u=u(ind);
0046       t=t(ind);
0047       u=u(:);
0048       t=t(:);
0049 %------ set up A -------
0050        A=zeros(nind,nfreq*2+1); 
0051        A(:,1)=ones(nind,1);
0052        for j=[1:nfreq]
0053            A(:,2*j)=cos(freq(j)*t);
0054            A(:,2*j+1)=sin(freq(j)*t);
0055        end
0056 %-------solve [A coeff = u] -----------------
0057        coef=A\u;
0058 %-------generate solution components-----
0059        up=zeros(nt,nfreq*2+1);
0060        up(:,1)=coef(1)*ones(nt,1);
0061        for j=[1:nfreq]
0062           up(:,j+1)=coef(j*2)*cos(freq(j)*t)+coef(j*2+1)*sin(freq(j)*t);
0063        end
0064        up(:,nfreq+2)=sum(up(:,[1:nfreq+1])')';  % sum of all comps
0065        unew(:,nn)=up(:,nfreq+2);
0066          ucoef(:,nn)=coef;
0067          error=std(A*coef-u)/sqrt(nind-length(periods)*2+1);
0068        if exist('goplot')
0069           jdp=(t-t(1))/24;
0070         h=plot(jdp,[u unew(:,nn)],jdp,zeros(size(jdp)),'k'); 
0071         hold on
0072         plot(jdp,u,'x',jdp,unew(:,nn),'o');
0073         hold off
0074         legend('data','tide fit')
0075           title(sprintf...
0076             ('Least Squares Fit, series %g   error= %5.2g',nn,error));
0077           xlabel('Days from start of record');
0078           pause
0079        end
0080 else
0081        unew(:,nn)=ones(t)*nan;
0082          ucoef(:,nn)=ones(2*nfreq+1,1)*nan;
0083 end   % if nind> ....
0084 end    %  loop

Generated on Thu 01-Dec-2005 07:46:03 by m2html © 2003