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)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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