matlab t tide,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).% % (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 vectornobs=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 pointt=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 frequenciesmi=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 analysiselse % 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));