t_xtide.m

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

function pred=t_xtide(varargin);
% T_XTIDE Tidal prediction
% YOUT=T_XTIDE(STATION) makes a tidal prediction
% for the current day using the harmonics file from XTIDE. 
% if STATION is a string then the first match found in the database is
% used, you can request matches to other stations by appending '(2)' to
% the string. If you don't know the station name but want to find the
% nearest to a given LONG,LAT, try T_XTIDE(LONG,LAT).
%
% The times of predicted tides are given by the next numerical argument
% (if any), e.g. [...]=T_XTIDE(STATION,TIM). 
% TIM can be: a vector of matlab-format decimal days (from DATENUM).
%           : a scalar <1000, taken as the number of days from present
%           : a scalar >1000, taken as the starting time in matlab-format 
%             for a 2 day time series. 
%           : not given, in which case the current time is used as a start 
%             time.
%
% Times are usually taken to be in 
% standard time for the given location (no daylight savings adjustments); 
% if in doubt use the 'info' or 'full' options where offset from UTC is given.
%
% 
% Other optional arguments can be specified following this using 
% property/value pairs: 
%
%     'format'     'raw' (default)
%                    YOUT is just a time series to match the time in TIM
%
%                  'times'
%                    YOUT is a structure of high/low tide information
%                    between times min(TIM) and max(TIM).
%
%                  'info'
%                    YOUT is a structure giving station information
%                    (location, time zone, units, harmonic constituents)
%
%                  'full'
%                    Combination of 'raw' and 'info' in a structure YOUT.
%
%     'units'     {'meters' | 'feet' | 'm/s' | 'knots' | 'original' }
%                    Units of result (default is original units)
%
% If no output argument is specified data is plotted and/or displayed.
%
% If the chosen name matches several stations then the first in the list is
% chosen. Specific choices can be made by appending a '(2)' or '(3)' (etc.)
% to the name, e.g.  T_XTIDE('tofino (2)',...).
%
%  Requires the xtide harmonics file  - get this from 
%            http://bel-marduk.unh.edu/xtide/files.html

% R. Pawlowicz 1/Dec/2001
% Version 1.0
%          16/May/02 - added lat/long options (thanks to Richard Dewey).


% Get the harmonics data from a) a mat-file if it exists, or b) from a harmonics
% file.

if ~exist('t_xtide.mat','file'), % Read the harmonics file and make a mat file

  filnam='/usr/share/xtide/harmonics.txt';
  
  fprintf('\n********Can''t find mat-file t_xtide.mat ********\n\n');
  fprintf('Attempting to generate one from an xtide harmonics file....\n\n');
  fprintf('Latest version available from http://bel-marduk.unh.edu/xtide/files.html\n\n');
  
  % Input name
  fid=-1;
  while fid==-1,
    rep=filnam;
    while (lower(rep(1))~='y'),
     filnam=rep;
     rep='n';
     rep=input(['Harmonics filename: ' filnam '? (y/Y/new file name):'],'s');
     if isempty(rep), rep='y'; end;
    end; 
    
    fid=fopen(filnam);
    if fid==-1,
      fprintf(['\n****** Can''t open filename ->' filnam '<-\n\n']);
    end;
  end;
    
  fprintf('Reading harmonics file (this will take a while)\n');
  [xtide,xharm]=read_xtidefile(fid);
  
  fprintf('Saving harmonic information to t_xtide.mat\n');
  save t_xtide xtide xharm
   
else
  load t_xtide
end;

if nargin>0,

  if isstr(varargin{1}),  % Station  name given
    % Identify station - look for exact match first
    ista=strmatch(lower(varargin{1}),lower(xharm.station),'exact');
    % otherwise go for partial matches
    if isempty(ista),
      % First check to see if a number was selected:
      inum=-10;
      while inum<-1,
        inum=inum+1;
        ll=findstr(lower(varargin{1}),sprintf('(%d)',-inum));
        if ~isempty(ll),
          inum=abs(inum);
	      varargin{1}=deblank(varargin{1}(1:ll-1));
        end;
      end;  
      ista=strmatch(lower(varargin{1}),lower(xharm.station));
      if length(ista)>1,
        if inum>0 & inum<=length(ista),
          ista=ista(inum);
        else	
          fprintf('Ambiguous Station Choice - Taking first of:\n');
          for kk=1:length(ista),
	        fprintf('%5d: %s\n',ista(kk),deblank(xharm.station(ista(kk),:)));
	        fprintf('      Long: %.4f  Lat: %.4f \n',xharm.longitude(ista(kk)),xharm.latitude(ista(kk)));
          end;
          fprintf('\n');
          ista=ista(1);
        end 	
      elseif length(ista)==1 & inum>1,
        fprintf('***Can''t find variant (%d) of station - Taking only choice\n',inum);
      elseif length(ista)==0,  
        error('Could not match station');
      end;    
      varargin(1)=[];
     end;
  
   else   % Lat/long?
      [dist,hdg]=t_gcdist(xharm.latitude,xharm.longitude,varargin{2},varargin{1});
      [mind,ista]=min(dist);
      if length(ista)>1,
        fprintf('Ambiguous Station Choice - Taking first of:\n');
        for kk=1:length(ista),
	      fprintf('%5d: %s\n',ista(kk),deblank(xharm.station(ista(kk),:)));
	      fprintf('      Long: %.4f  Lat: %.4f \n',xharm.longitude(ista(kk)),xharm.latitude(ista(kk)));
        end;
        fprintf('\n');
        ista=ista(1);
      else
 	    fprintf('%5d: %s\n',ista,deblank(xharm.station(ista,:)));
	    fprintf('      Long: %.4f  Lat: %.4f \n',xharm.longitude(ista),xharm.latitude(ista)); 
      end;
      varargin(1:2)=[];
   end;
  
  % Time vector (if available) otherwise take current time.

  if length(varargin)>0 & ~isstr(varargin{1}),
    tim=varargin{1};
    tim=tim(:)';
    varargin(1)=[];
    if length(tim)==1,
      if tim<1000,
        dat=clock;
        tim=datenum(dat(1),dat(2),dat(3))+[0:1/48:tim];
      else
        tim=tim+[0:1/48:2]; % 2 days worth.
      end;	 	
    end;
  else 
    dat=clock;
    tim=datenum(dat(1),dat(2),dat(3))+[0:.25:48]/24;
  end;
   
  % Parse properties

  format='raw';
  unt='original';
  
  k=1;
  while length(varargin)>0,
      switch lower(varargin{1}(1:3)),
	case 'for',
	 format=lower(varargin{2});
	case 'uni',
	 unt=lower(varargin{2}); 
	otherwise,
           error(['Can''t understand property:' varargin{1}]);
      end;
      varargin([1 2])=[]; 
  end;
 
  % if we want a time series
  pred=[];
  % Convert units if requested.
  [units,convf]=convert_units(unt,xharm.units(ista,:));
  if strcmp(format(1:2),'ra') | strcmp(format(1:2),'fu') | strcmp(format(1:2),'ti')
    
    % Data every minute for hi/lo forecasting.
    if strcmp(format(1:2),'ti'),
      tim=tim(1):(1/1440):tim(end); 
    end;

    % Convert into time since the beginning of year
    mid=datevec(mean(tim));
    iyr=mid(1)-xtide.startyear+1;
    lt=length(tim);
    xtim=(tim-datenum(mid(1),1,1))*24; % Hours since beginning of year

    %-----------------------------------------------------
    % Sum up everything for the prediction!

    pred=xharm.datum(ista)+sum( ...
      repmat(xtide.nodefactor(:,iyr).*xharm.A(ista,:)',1,lt).* ...
      cos( ( xtide.speed*xtim + repmat(xtide.equilibarg(:,iyr)-xharm.kappa(ista,:)',1,lt) )*(pi/180) ),1);
    %-----------------------------------------------------

    pred=pred*convf;
    
    % Compute times of hi/lo from every-minute data
    if strcmp(format(1:2),'ti'),
     % Check if this is a current station
      if ~isempty(findstr('Current',xharm.station(ista,:))), currents=1; else currents=0; end;
      dpred=diff(pred);
      ddpred=diff(dpred>0);

      flat=find(ddpred~=0)+1;
      slk=find(sign(pred(1:end-1))~=sign(pred(2:end)));
      
      hi.mtime=tim(flat);
      hi.value=pred(flat);

      hi.type=zeros(size(flat));
      hi.type(find(ddpred(flat-1)<0))=1;  % 0=lo, 1=hi
      hi.units=deblank(units);
      
      pred=hi;
    end;
  end;
  
  % Create information structure
  
  if strcmp(format(1:2),'in') | strcmp(format(1:2),'fu'),
    if ~isempty(pred), 
      pred.yout=pred; 
      pred.mtime=tim; 
    else
      kk=find(xharm.A(ista,:)~=0);
      pred.freq=xtide.name(kk,:);
      pred.A=full(xharm.A(ista,kk)')*convf;
      pred.kappa=full(xharm.kappa(ista,kk)'); 
    end;
    pred.station=deblank(xharm.station(ista,:));
    pred.longitude=xharm.longitude(ista);
    pred.latitude=xharm.latitude(ista);
    pred.timezone=xharm.timezone(ista);
    pred.units=deblank(units);
    pred.datum=xharm.datum(ista)*convf;
  end;
 
end;

% If no output parameters then we plot or display things

if nargout==0,