t_tide.m

用Matlab编写的一款计算调和常数的程序包

function [nameu,fu,tidecon,xout]=t_tide(xin,varargin);
% T_TIDE Harmonic analysis of a time series
% [NAME,FREQ,TIDECON,XOUT]=T_TIDE(XIN) computes the tidal analysis 
% of the (possibly complex) time series XIN.
%
% [TIDESTRUC,XOUT]=T_TIDE(XIN) returns the analysis information in
% a structure formed of NAME, FREQ, and TIDECON.
%
% XIN can be scalar (e.g. for elevations), or complex ( =U+sqrt(-1)*V
% for eastward velocity U and northward velocity V.
%
% Further inputs are optional, and are specified as property/value pairs
% [...]=T_TIDE(XIN,property,value,property,value,...,etc.)
%      
% These properties are:
%
%       'interval'       Sampling interval (hours), default = 1. 
%          
%   The next two are required if nodal corrections are to be computed,
%   otherwise not necessary. If they are not included then the reported
%   phases are raw constituent phases at the central time. 
%
%   If your time series is longer than 18.6 years then nodal corrections
%   are not made -instead we fit directly to all satellites (start time
%   is then just used to generate Greenwich phases).
%
%       'start time'     [year,month,day,hour,min,sec]
%                        - min,sec are optional OR 
%                        decimal day (matlab DATENUM scalar)
%       'latitude'       decimal degrees (+north) (default: none).
%
%   Where to send the output.
%       'output'         where to send printed output:
%                        'none'    (no printed output)
%                        'screen'  (to screen) - default
%                        FILENAME   (to a file)
%
%   Correction factor for prefiltering.
%       'prefilt'        FS,CORR
%                        If the time series has been passed through
%                        a pre-filter of some kind (say, to reduce the
%                        low-frequency variability), then the analyzed
%                        constituents will have to be corrected for 
%                        this. The correction transfer function 
%                        (1/filter transfer function) has (possibly 
%                        complex) magnitude CORR at frequency FS (cph). 
%                        Corrections of more than a factor of 100 are 
%                        not applied; it is assumed these refer to tidal
%                        constituents that were intentionally filtered 
%                        out, e.g., the fortnightly components.
%
%   Adjustment for long-term behavior ("secular" behavior).
%       'secular'        'mean'   - assume constant offset (default).
%                        'linear' - get linear trend.
%                     
%   Inference of constituents.
%       'inference'      NAME,REFERENCE,AMPRAT,PHASE_OFFSET
%                        where NAME is an array of the names of 
%                        constituents to be inferred, REFERENCE is an 
%                        array of the names of references, and AMPRAT 
%                        and PHASE_OFFSET are the amplitude factor and
%                        phase offset (in degrees)from the references. 
%                        NAME and REFERENCE are Nx4 (max 4 characters
%                        in name), and AMPRAT and PHASE_OFFSET are Nx1
%                        (for scalar time series) and Nx2 for vector 
%                        time series (column 1 is for + frequencies and
%                        column 2 for - frequencies).
%
%   Shallow water constituents
%       'shallow'        NAME
%                        A matrix whose rows contain the names of 
%                        shallow-water constituents to analyze.
%
%   Resolution criterions for least-squares fit.        
%       'rayleigh'       scalar - Rayleigh criteria, default = 1.
%                        Matrix of strings - names of constituents to
%                                   use (useful for testing purposes).
%  
%   Calculation of confidence limits.
%       'error'          'wboot'  - Boostrapped confidence intervals 
%                                   based on a correlated bivariate 
%                                   white-noise model.
%                        'cboot'  - Boostrapped confidence intervals 
%                                   based on an uncorrelated bivariate 
%                                   coloured-noise model (default).
%                        'linear' - Linearized error analysis that 
%                                   assumes an uncorrelated bivariate 
%                                   coloured noise model. 
%                                   
%   Computation of "predicted" tide (passed to t_predic, but note that
%                                    the default value is different).
%       'synthesis'      0 - use all selected constituents
%                        scalar>0 - use only those constituents with a 
%                                   SNR greater than that given (1 or 2 
%                                   are good choices, 2 is the default).
%                              <0 - return result of least-squares fit 
%                                   (should be the same as using '0', 
%                                   except that NaN-holes in original 
%                                   time series will remain and mean/trend
%                                   are included).
%
%   Least squares soln computational efficiency parameter
%	'lsq'		'direct'  - use A\x fit
%			'normal'  - use (A'A)\(A'x) (may be necessary
%				    for very large input vectors since
%                                   A'A is much smaller than A)
%			'best'	  - automatically choose based on
%				    length of series (default).
%
%       It is possible to call t_tide without using property names,
%       in which case the assumed calling sequence is
%
%          T_TIDE(XIN,INTERVAL,START_TIME,LATITUDE,RAYLEIGH)
%
%
%  OUTPUT: 
%
%    nameu=list of constituents used
%    fu=frequency of tidal constituents (cycles/hr)
%    tidecon=[fmaj,emaj,fmin,emin,finc,einc,pha,epha] for vector xin
%           =[fmaj,emaj,pha,epha] for scalar (real) xin
%       fmaj,fmin - constituent major and minor axes (same units as xin)       
%       emaj,emin - 95% confidence intervals for fmaj,fmin
%       finc - ellipse orientations (degrees)
%       einc - 95% confidence intervals for finc
%       pha - constituent phases (degrees relative to Greenwich)
%       epha - 95% confidence intervals for pha
%    xout=tidal prediction
%
% Note: Although missing data can be handled with NaN, it is wise not
%       to have too many of them. If your time series has a lot of 
%       missing data at the beginning and/or end, then truncate the 
%       input time series.  The Rayleigh criterion is applied to 
%       frequency intervals calculated as the inverse of the input 
%       series length.
%
% A description of the theoretical basis of the analysis and some
% implementation details can be found in:
%
% Pawlowicz, R., B. Beardsley, and S. Lentz, "Classical Tidal 
%   "Harmonic Analysis Including Error Estimates in MATLAB 
%    using T_TIDE", Computers and Geosciences, 28, 929-937 (2002).
%
% (citation of this article would be appreciated if you find the
%  toolbox useful).


% R. Pawlowicz 11/8/99 - Completely rewritten from the transliterated-
%                        to-matlab IOS/Foreman fortran code by S. Lentz
%                        and B. Beardsley.
%              3/3/00  - Redid errors to take into account covariances 
%                        between u and v errors.
%              7/21/00 - Found that annoying bug in error calc! 
%              11/1/00 - Added linear error analysis.
%              8/29/01 - Made synth=1 default, also changed behavior 
%                        when no lat/time given so that phases are raw
%                        at central time. 
%              9/1/01  - Moved some SNR code to t_predic.
%              9/28/01 - made sure you can't choose Z0 as constituent.
%              6/12/01 - better explanation for variance calcs, fixed
%                        bug in typed output (thanks Mike Cook).
%              8/2/03 - Added block processing for long time series (thanks
%                       to Derek Goring).
%              9/2/03 - Beta version of 18.6 year series handling
%              12/2/03 - Bug (x should be xin) fixed thanks to Mike Cook (again!)

%
% Version 1.2



% ----------------------Parse inputs-----------------------------------

ray=1;
dt=1;
fid=1;
stime=[];
lat=[];
corr_fs=[0 1e6];
corr_fac=[1  1];
secular='mean';
inf.iname=[];
inf.irefname=[];
shallownames=[];
constitnames=[];
errcalc='cboot';
synth=2;
lsq='best';

k=1;
while length(varargin)>0,
  if ischar(varargin{1}),
    switch lower(varargin{1}(1:3)),
      case 'int',
        dt=varargin{2};
      case 'sta',
        stime=varargin{2};
	if length(stime)>1, 
	  stime=[stime(:)' zeros(1,6-length(stime))]; 
	  stime=datenum(stime(1),stime(2),stime(3),stime(4),stime(5),stime(6));
	end;
      case 'lat',
         lat=varargin{2};
      case 'out',
         filen=varargin{2};
	 switch filen,
	   case 'none',
	     fid=-1;
	   case 'screen',
	     fid=1;
	   otherwise
	     [fid,mesg]=fopen(filen,'w');
	     if fid==-1, error(msg); end;
	  end;
      case 'ray',
         if isnumeric(varargin{2}),
           ray=varargin{2};
	 else
	   constitnames=varargin{2};
	   if iscellstr(constitnames), constitnames=char(constitnames); end;
	 end;
       case 'pre',
         corr_fs=varargin{2};
	 corr_fac=varargin{3};
         varargin(1)=[];
      case 'sec',
         secular=varargin{2};
      case 'inf',
         inf.iname=varargin{2};
	 inf.irefname=varargin{3};
	 inf.amprat=varargin{4};
	 inf.ph=varargin{5};
	 varargin(1:3)=[];
      case 'sha',
         shallownames=varargin{2};
      case 'err',
         errcalc=varargin{2};
      case 'syn',
         synth=varargin{2};
      case 'lsq',
         lsq=varargin{2};	 
      otherwise,
         error(['Can''t understand property:' varargin{1}]);
    end;
    varargin([1 2])=[]; 
  else  
    switch k,
      case 1,
        dt=varargin{1};
      case 2,
        stime=varargin{1};
      case 3,
        lat=varargin{1};
      case 4,
        ray=varargin{1};
      otherwise
        error('Too many input parameters');
     end;
     varargin(1)=[];
  end;
  k=k+1;
end;
 
[inn,inm]=size(xin);
if ~(inn==1 | inm==1), error('Input time series is not a vector'); end;

xin=xin(:); % makes xin a column vector
nobs=length(xin);

if strcmp(lsq(1:3),'bes'),  % Set matrix method if auto-choice.
 if nobs>10000,
    lsq='normal';
 else
    lsq='direct';
 end;
end;
 
if nobs*dt> 18.6*365.25*24,  % Long time series
  longseries=1; ltype='full';
else
  longseries=0; ltype='nodal';
end;
        		
nobsu=nobs-rem(nobs-1,2);% makes series odd to give a center point

t=dt*([1:nobs]'-ceil(nobsu/2));  % Time vector for entire time series,
                                 % centered at series midpoint. 

if ~isempty(stime),
  centraltime=stime+floor(nobsu./2)./24.0*dt;
else
  centraltime=[];
end;

% -------Get the frequencies to use in the harmonic analysis-----------

[nameu,fu,ju,namei,fi,jinf,jref]=constituents(ray/(dt*nobsu),constitnames,...
                                           shallownames,inf.iname,inf.irefname,centraltime);

mu=length(fu); % # base frequencies
mi=length(fi); % # inferred

% Find the good data points (here I assume that in a complex time 
% series, if u is bad, so is v).

gd=find(finite(xin(1:nobsu)));
ngood=length(gd);
fprintf('   Points used: %d of %d\n',ngood,nobs)



%----------------------------------------------------------------------
% Now solve for the secular trend plus the analysis. Instead of solving
% for + and - frequencies using exp(i*f*t), I use sines and cosines to 
% keep tc real.  If the input series is real, than this will 
% automatically use real-only computation (faster). However, for the analysis, 
% it's handy to get the + and - frequencies ('ap' and 'am'), and so 
% that's what we do afterwards.

% The basic code solves the matrix problem Ac=x+errors where the functions to
% use in the fit fill up the A matrix, which is of size (number points)x(number
% constituents). This can get very, very large for long time series, and
% for this the more complex block processing algorithm was added. It should
% give identical results (up to roundoff error)

if strcmp(lsq(1:3),'dir'),

  if secular(1:3)=='lin',
    tc=[ones(length(t),1) cos((2*pi)*t*fu') sin((2*pi)*t*fu') t*(2/dt/nobsu)];
  else
    tc=[ones(length(t),1) cos((2*pi)*t*fu') sin((2*pi)*t*fu') ];
  end;
  
  coef=tc(gd,:)\xin(gd);

  z0=coef(1);
  ap=(coef(2:(1+mu))-i*coef((2+mu):(1+2*mu)))/2;  % a+ amplitudes
  am=(coef(2:(1+mu))+i*coef((2+mu):(1+2*mu)))/2;  % a- amplitudes
  if secular(1:3)=='lin',
    dz0=coef(end);
  else
    dz0=0;
  end;    
  xout=tc*coef;  % This is the time series synthesized from the analysis

else  % More complicated code required for long time series when memory may be
      % a problem. Modified from code submitted by Derek Goring (NIWA Chrischurch)
      
      % Basically the normal equations are formed (rather than using Matlab's \
      % algorithm for least squares); this can be done by adding up subblocks
      % of data. Notice how the code is messier, and we have to recalculate everything
      % to get the original fit.

  nsub=5000;  % Block length - doesn't matter really but should be small enough to
              % get allocated quickly	      
  if secular(1:3)=='lin',
    lhs=zeros(mu*2+2,mu*2+2); rhs=zeros(mu*2+2,1);
    for j1=1:nsub:ngood
      j2=min(j1 + nsub - 1,ngood);
      E=[ones(j2-j1+1,1) cos((2*pi)*t(gd(j1:j2))*fu') sin((2*pi)*t(gd(j1:j2))*fu') t(gd(j1:j2))*(2/dt/nobsu)];
      rhs=rhs + E'*xin(gd(j1:j2));
      lhs=lhs + E'*E;
    end;
  else  
    lhs=zeros(mu*2+1,mu*2+1); rhs=zeros(mu*2+1,1);
    for j1=1:nsub:ngood
      j2=min(j1 + nsub - 1,ngood);
      E=[ones(j2-j1+1,1) cos((2*pi)*t(gd(j1:j2))*fu') sin((2*pi)*t(gd(j1:j2))*fu')];
      rhs=rhs + E'*xin(gd(j1:j2));