t_tide.m

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

% Compute frequencies from astronomical considerations.

if minres>1/(18.6*365.25*24),                       % Choose only resolveable pairs for short
  [const,sat,cshallow]=t_getconsts(centraltime);    % Time series   
  ju=find(const.df>=minres);
else                                                % Choose them all if > 18.6 years.
  [const,sat,cshallow]=t_get18consts(centraltime);
  ju=[2:length(const.freq)]';  % Skip Z0
  for ff=1:2,                  % loop twice to make sure of neightbouring pairs
    jck=find(diff(const.freq(ju))<minres);
    if (length(jck)>0)
       jrm=jck;
       jrm=jrm+(abs(const.doodsonamp(ju(jck+1)))<abs(const.doodsonamp(ju(jck))));
       disp(['  Warning! Following constituent pairs violate Rayleigh criterion']);
       for ick=1:length(jck);
	 disp(['     ',const.name(ju(jck(ick)),:),' vs ',const.name(ju(jck(ick)+1),:) ' - not using ',const.name(ju(jrm(ick)),:)]);
       end;
       ju(jrm)=[];
    end
  end;
end;
  
if ~isempty(constit),     % Selected if constituents are specified in input.
  ju=[];
  for k=1:size(constit,1),
   j1=strmatch(constit(k,:),const.name);
   if isempty(j1),
     disp(['Can''t recognize name ' constit(k,:) ' for forced search']);
   elseif j1==1,
     disp(['*************************************************************************']);
     disp(['Z0 specification ignored - for non-tidal offsets see ''secular'' property']);
     disp(['*************************************************************************']);
   else  
     ju=[ju;j1];
   end;
  end;
  [dum,II]=sort(const.freq(ju)); % sort in ascending order of frequency.
  ju=ju(II);
end;


disp(['   number of standard constituents used: ',int2str(length(ju))])

if ~isempty(shallow),          % Add explictly selected shallow water constituents.
 for k=1:size(shallow,1),
   j1=strmatch(shallow(k,:),const.name);
   if isempty(j1),
     disp(['Can''t recognize name ' shallow(k,:) ' for forced search']);
   else
     if isnan(const.ishallow(j1)),
       disp([shallow(k,:) ' Not a shallow-water constituent']);
     end;
     disp(['   Forced fit to ' shallow(k,:)]);
     ju=[ju;j1];
   end;
 end;
 
end;
      
nameu=const.name(ju,:);
fu=const.freq(ju);


% Check if neighboring chosen constituents violate Rayleigh criteria.
jck=find(diff(fu)<minres);
if (length(jck)>0)
   disp(['  Warning! Following constituent pairs violate Rayleigh criterion']);
   for ick=1:length(jck);
   disp(['     ',nameu(jck(ick),:),'  ',nameu(jck(ick)+1,:)]);
   end;
end

% For inference, add in list of components to be inferred.

fi=[];namei=[];jinf=[];jref=[];
if ~isempty(infname),
  fi=zeros(size(infname,1),1);
  namei=zeros(size(infname,1),4);
  jinf=zeros(size(infname,1),1)+NaN;
  jref=zeros(size(infname,1),1)+NaN;

  for k=1:size(infname,1),
   j1=strmatch(infname(k,:),const.name);
   if isempty(j1),
     disp(['Can''t recognize name' infname(k,:) ' for inference']);
   else
    jinf(k)=j1;
    fi(k)=const.freq(j1);
    namei(k,:)=const.name(j1,:);
    j1=strmatch(infref(k,:),nameu);
    if isempty(j1),
      disp(['Can''t recognize name ' infref(k,:) ' for as a reference for inference']);
    else
      jref(k)=j1;
      fprintf(['   Inference of ' namei(k,:) ' using ' nameu(j1,:) '\n']);
    end;
   end;
  end;    
  jinf(isnan(jref))=NaN;
end;

%----------------------------------------------------------------------
function y=fixgaps(x);
% FIXGAPS: Linearly interpolates gaps in a time series
% YOUT=FIXGAPS(YIN) linearly interpolates over NaN in the input time 
% series (may be complex), but ignores trailing and leading NaNs.

% R. Pawlowicz 11/6/99
% Version 1.0

y=x;

bd=isnan(x);
gd=find(~bd);

bd([1:(min(gd)-1) (max(gd)+1):end])=0;


y(bd)=interp1(gd,x(gd),find(bd)); 


%----------------------------------------------------------------------
function ain=cluster(ain,clusang);
% CLUSTER: Clusters angles in rows around the angles in the first 
% column. CLUSANG is the allowable ambiguity (usually 360 degrees but
% sometimes 180).

ii=(ain-ain(:,ones(1,size(ain,2))))>clusang/2;
ain(ii)=ain(ii)-clusang;
ii=(ain-ain(:,ones(1,size(ain,2))))<-clusang/2;
ain(ii)=ain(ii)+clusang;


%----------------------------------------------------------------------
function [NP,NM]=noise_realizations(xres,fu,dt,nreal,errcalc);
% NOISE_REALIZATIONS: Generates matrices of noise (with correct
% cross-correlation structure) for bootstrap analysis.
%

% R. Pawlowicz 11/10/00
% Version 1.0

if strmatch(errcalc,'cboot'),
  [fband,Pxrave,Pxiave,Pxcave]=residual_spectrum(xres,fu,dt);
  
  Pxcave=zeros(size(Pxcave));  %% For comparison with other technique!
  %fprintf('**** Assuming no covariance between u and v errors!*******\n');

elseif strmatch(errcalc,'wboot'),
  fband=[0 .5];
  nx=length(xres);
  A=cov(real(xres),imag(xres))/nx;
  Pxrave=A(1,1);Pxiave=A(2,2);Pxcave=A(1,2);
else
  error(['Unrecognized type of bootstap analysis specified: ''' errcalc '''']);
end;
  
nfband=size(fband,1);

Mat=zeros(4,4,nfband);
for k=1:nfband,

  % The B matrix represents the covariance matrix for the vector
  % [Re{ap} Im{ap} Re{am} Im{am}]' where Re{} and Im{} are real and
  % imaginary parts, and ap/m represent the complex constituent 
  % amplitudes for positive and negative frequencies when the input
  % is bivariate white noise. For a flat residual spectrum this works 
  % fine.
 
  % This is adapted here for "locally white" conditions, but I'm still
  % not sure how to handle a complex sxy, so this is set to zero
  % right now.
  
  p=(Pxrave(k)+Pxiave(k))/2;
  d=(Pxrave(k)-Pxiave(k))/2;
  sxy=Pxcave(k);
  
  B=[p    0   d   sxy;
     0    p  sxy  -d;
     d   sxy  p    0
     sxy -d   0    p];

  % Compute the transformation matrix that takes uncorrelated white 
  % noise and makes noise with the same statistical structure as the 
  % Fourier transformed noise.
  [V,D]=eig(B);
  Mat(:,:,k)=V*diag(sqrt(diag(D)));
end;

% Generate realizations for the different analyzed constituents.

N=zeros(4,nreal);
NM=zeros(length(fu),nreal);
NP=NM;
for k=1:length(fu);
  l=find(fu(k)>fband(:,1) & fu(k)<fband(:,2));
  N=[zeros(4,1),Mat(:,:,l)*randn(4,nreal-1)];
  NP(k,:)=N(1,:)+i*N(2,:);
  NM(k,:)=N(3,:)+i*N(4,:);
end;

%----------------------------------------------------------------------
function [ercx,eicx]=noise_stats(xres,fu,dt);
% NOISE_STATS: Computes statistics of residual energy for all 
% constituents (ignoring any cross-correlations between real and
% imaginary parts).

% S. Lentz  10/28/99
% R. Pawlowicz 11/1/00
% Version 1.0

[fband,Pxrave,Pxiave,Pxcave]=residual_spectrum(xres,fu,dt);
nfband=size(fband,1);
mu=length(fu);

% Get the statistics for each component.
ercx=zeros(mu,1);
eicx=zeros(mu,1);
for k1=1:nfband;
   k=find(fu>=fband(k1,1) & fu<=fband(k1,2));
   ercx(k)=sqrt(Pxrave(k1));
   eicx(k)=sqrt(Pxiave(k1));
end

%----------------------------------------------------------------------
function [fband,Pxrave,Pxiave,Pxcave]=residual_spectrum(xres,fu,dt)
% RESIDUAL_SPECTRUM: Computes statistics from an input spectrum over
% a number of bands, returning the band limits and the estimates for
% power spectra for real and imaginary parts and the cross-spectrum.          
%
% Mean values of the noise spectrum are computed for the following 
% 8 frequency bands defined by their center frequency and band width:
% M0 +.1 cpd; M1 +-.2 cpd; M2 +-.2 cpd; M3 +-.2 cpd; M4 +-.2 cpd; 
% M5 +-.2 cpd; M6 +-.21 cpd; M7 (.26-.29 cpd); and M8 (.30-.50 cpd). 

% S. Lentz  10/28/99
% R. Pawlowicz 11/1/00
% Version 1.0

% Define frequency bands for spectral averaging.
fband =[.00010 .00417;
        .03192 .04859;
        .07218 .08884;
        .11243 .12910;
        .15269 .16936;
        .19295 .20961;
        .23320 .25100;
        .26000 .29000;
        .30000 .50000];

% If we have a sampling interval> 1 hour, we might have to get
% rid of some bins.
%fband(fband(:,1)>1/(2*dt),:)=[];

nfband=size(fband,1);
nx=length(xres);

% Spectral estimate (takes real time series only).

[Pxr,fx]=psd(real(xres),nx,1/dt); % Call to SIGNAL PROCESSING TOOLBOX - see note in t_readme. If you have an error here you are probably missing this toolbox
[Pxi,fx]=psd(imag(xres),nx,1/dt); % Call to SIGNAL PROCESSING TOOLBOX - see note in t_readme.
[Pxc,fx]=csd(real(xres),imag(xres),nx,1/dt); % Call to SIGNAL PROCESSING TOOLBOX - see note in t_readme.


df=fx(3)-fx(2);
Pxr(round(fu./df)+1)=NaN ; % Sets Px=NaN in bins close to analyzed frequencies
Pxi(round(fu./df)+1)=NaN ; % (to prevent leakage problems?).
Pxc(round(fu./df)+1)=NaN ; 

Pxrave=zeros(nfband,1);
Pxiave=zeros(nfband,1);
Pxcave=zeros(nfband,1);
% Loop downwards in frequency through bands (cures short time series
% problem with no data in lowest band).
%
% Divide by nx to get power per frequency bin, and multiply by 2
% to account for positive and negative frequencies.
%
for k=nfband:-1:1,
   jband=find(fx>=fband(k,1) & fx<=fband(k,2) & finite(Pxr));
   if any(jband),
     Pxrave(k)=mean(Pxr(jband))*2/nx;
     Pxiave(k)=mean(Pxi(jband))*2/nx;
     Pxcave(k)=mean(Pxc(jband))*2/nx;
   elseif k<nfband,
     Pxrave(k)=Pxrave(k+1);   % Low frequency bin might not have any points...
     Pxiave(k)=Pxiave(k+1);   
     Pxcave(k)=Pxcave(k+1);   
   end;
end



%----------------------------------------------------------------------
function [emaj,emin,einc,epha]=errell(cxi,sxi,ercx,ersx,ercy,ersy)
% [emaj,emin,einc,epha]=errell(cx,sx,cy,sy,ercx,ersx,ercy,ersy) computes
% the uncertainities in the ellipse parameters based on the 
% uncertainities in the least square fit cos,sin coefficients.
%
%  INPUT:  cx,sx=cos,sin coefficients for x 
%          cy,sy=cos,sin coefficients for y
%          ercx,ersx=errors in x cos,sin coefficients
%          ercy,ersy=errors in y cos,sin coefficients
%          
%  OUTPUT: emaj=major axis error
%          emin=minor axis error
%          einc=inclination error (deg)
%          epha=pha error (deg)

% based on linear error propagation, with errors in the coefficients 
% cx,sx,cy,sy uncorrelated. 

% B. Beardsley  1/15/99; 1/20/99
% Version 1.0

r2d=180./pi;
cx=real(cxi(:));sx=real(sxi(:));cy=imag(cxi(:));sy=imag(sxi(:));
ercx=ercx(:);ersx=ersx(:);ercy=ercy(:);ersy=ersy(:);

rp=.5.*sqrt((cx+sy).^2+(cy-sx).^2);
rm=.5.*sqrt((cx-sy).^2+(cy+sx).^2);
ercx2=ercx.^2;ersx2=ersx.^2;
ercy2=ercy.^2;ersy2=ersy.^2;

% major axis error
ex=(cx+sy)./rp;
fx=(cx-sy)./rm;
gx=(sx-cy)./rp;
hx=(sx+cy)./rm;
dcx2=(.25.*(ex+fx)).^2;
dsx2=(.25.*(gx+hx)).^2;
dcy2=(.25.*(hx-gx)).^2;
dsy2=(.25.*(ex-fx)).^2;
emaj=sqrt(dcx2.*ercx2+dsx2.*ersx2+dcy2.*ercy2+dsy2.*ersy2);

% minor axis error
dcx2=(.25.*(ex-fx)).^2;
dsx2=(.25.*(gx-hx)).^2;
dcy2=(.25.*(hx+gx)).^2;
dsy2=(.25.*(ex+fx)).^2;
emin=sqrt(dcx2.*ercx2+dsx2.*ersx2+dcy2.*ercy2+dsy2.*ersy2);

% inclination error
rn=2.*(cx.*cy+sx.*sy);
rd=cx.^2+sx.^2-(cy.^2+sy.^2);
den=rn.^2+rd.^2;
dcx2=((rd.*cy-rn.*cx)./den).^2;
dsx2=((rd.*sy-rn.*sx)./den).^2;
dcy2=((rd.*cx+rn.*cy)./den).^2;
dsy2=((rd.*sx+rn.*sy)./den).^2;
einc=r2d.*sqrt(dcx2.*ercx2+dsx2.*ersx2+dcy2.*ercy2+dsy2.*ersy2);

% phase error
rn=2.*(cx.*sx+cy.*sy);
rd=cx.^2-sx.^2+cy.^2-sy.^2;
den=rn.^2+rd.^2;
dcx2=((rd.*sx-rn.*cx)./den).^2;
dsx2=((rd.*cx+rn.*sx)./den).^2;
dcy2=((rd.*sy-rn.*cy)./den).^2;
dsy2=((rd.*cy+rn.*sy)./den).^2;
epha=r2d.*sqrt(dcx2.*ercx2+dsx2.*ersx2+dcy2.*ercy2+dsy2.*ersy2);