t_xtide.m

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

  switch format(1:2),    
    case 'ti',
  
    fprintf('High/Low Predictions for %s\n',xharm.station(ista,:));
    fprintf('Time offset %.1f from UTC\n\n',xharm.timezone(ista));
    
    outstr=repmat(' ',length(flat),41);
    outstr(:,1:20)=datestr(hi.mtime);
    outstr(:,22:27)=reshape(sprintf('%6.2f',hi.value),6,length(flat))';
    if currents,
      ll=hi.type==1;
      outstr(ll,31:41)=repmat(' Flood Tide',sum(ll),1);
      ll=hi.type==0;
      outstr(ll,31:41)=repmat(' Ebb Tide  ',sum(ll),1);
    else
      ll=hi.type==1;
      outstr(ll,31:41)=repmat(' High Tide ',sum(ll),1);
      ll=hi.type==0;
      outstr(ll,31:41)=repmat(' Low Tide  ',sum(ll),1);
    end;
    disp(outstr)   
           
    case 'ra'
     plot(tim,pred)
     datetick;
     title(['Tidal prediction for ',deblank(xharm.station(ista,:)) ' beginning ' datestr(tim(1))]); 
     ylabel(deblank(xharm.units(ista,:)));

    case 'fu'
     plot(tim,pred.yout)
     datetick;
     title(['Tidal prediction for ',deblank(xharm.station(ista,:)) ' beginning ' datestr(tim(1))]); 
     ylabel(deblank(xharm.units(ista,:)));
     
    case 'in',
        
    fprintf('Station: %s\n',pred.station);
    if pred.longitude<0, lon='W'; else lon='E'; end;
    if pred.latitude<0,  lat='S'; else lat='N'; end;
    fprintf('Location: %d %.1f'' %c, %d %.1f'' %c\n',fix(abs(pred.latitude)),rem(abs(pred.latitude),1)*60,...
         lat,fix(abs(pred.longitude)),rem(abs(pred.longitude),1)*60,lon);
    fprintf('Time offset %.1f from UTC\n\n',pred.timezone);
     	
   end;
   clear pred
end;  
  
%%%%%%%%%%%%%%%%%%%%
function [xtide,xharm]=read_xtidefile(fid);
% Reads the xtide harmonics file and creates a data structure
% with all that info for faster access


l=fgetl_nocom(fid);

ncon=sscanf(l,'%d');

xtide=struct('name',repmat(' ',ncon,8),'speed',zeros(ncon,1),...
	     'startyear',0,'equilibarg',zeros(ncon,68),'nodefactor',zeros(ncon,68));

for k=1:ncon,
 l=fgetl_nocom(fid);
 xtide.name(k,:)=l(1:8);
 xtide.speed(k)=sscanf(l(9:end),'%f');
end;

xtide.startyear=sscanf(fgetl_nocom(fid),'%d');

nyear=sscanf(fgetl_nocom(fid),'%d');

for k=1:ncon,
  l=fgetl(fid);
  xtide.equilibarg(k,:)=fscanf(fid,'%f',nyear);
  l=fgetl(fid);
end;
l=fgetl(fid); % Skip *END*

nyear=sscanf(fgetl_nocom(fid),'%d');

for k=1:ncon,
  l=fgetl(fid);
  xtide.nodefactor(k,:)=fscanf(fid,'%f',nyear);
  l=fgetl(fid);
end;
l=fgetl(fid); % Skip *END*

% Now read in all harmonic data


%nsta=1754; % This is number of stations in harmonics (1998-07-18)
%nsta=3351; % This is number of stations in v1.42 or harmonics file
nsta=3316; % This is number in v1.51

xharm=struct('station',repmat(' ',nsta,79),'units',repmat(' ',nsta,8),...
	     'longitude',zeros(nsta,1),'latitude',zeros(nsta,1),...
	     'timezone',zeros(nsta,1),'datum',zeros(nsta,1),...
	     'A',zeros(nsta,ncon),'kappa',zeros(nsta,ncon));

nh=0;
while length(l)>0 & l(1)~=-1,
 
  l=[l '   '];
  nh=nh+1;
  while ~strcmp(l(1:3),'# !'),
    l=[fgetl(fid) '   '];
  end;
  while strcmp(l(1:3),'# !'),
   switch l(4:7),
    case 'unit',
     tmp=deblank(l(findstr(l,':')+2:end));
     xharm.units(nh,1:length(tmp))=tmp;
    case 'long',
      xharm.longitude(nh)=sscanf(l(findstr(l,':')+1:end),'%f');
    case 'lati'  
      xharm.latitude(nh)=sscanf(l(findstr(l,':')+1:end),'%f');
   end;
   l=fgetl(fid);
  end; 
  tmp=deblank(l);
  if tmp(1)~='#', % Not commented out
    xharm.station(nh,1:length(tmp))=tmp;

    tmp=fgetl(fid);
    k=min(findstr(tmp,':'));
    tim=sscanf(tmp(1:k-1),'%d')+sscanf(tmp(k+[1:2]),'%d')/60;
    xharm.timezone(nh)=tim;

    tmp=fgetl(fid);
    xharm.datum(nh)=sscanf(tmp,'%f');

    for k=1:ncon,
      l=fgetl(fid);
      if l(1)~='x',
	ll=min([findstr(' ',l) find(abs(l)==9)]); % space or tab
	tmp=sscanf(l(ll+1:end),'%f',2);
	xharm.A(nh,k)=tmp(1);
	xharm.kappa(nh,k)=tmp(2);
      end;
    end;
    l=fgetl(fid);
  else
    nh=nh-1;  
  end;
  
  if rem(nh,50)==0, fprintf('.'); end;
end;
fprintf('\n');

% Convert internally to sparse matrix storage (much smaller).
xharm.A=sparse(xharm.A);
xharm.kappa=sparse(xharm.kappa);

return;
  
%%%%%%%%%%%%%%%%%%%%  
function l=fgetl_nocom(fid);
% Gets a line that isn't a comment line
%
l=fgetl(fid);
while length(l)>0 & l(1)=='#',
  l=fgetl(fid);
end;
  
%%%%%%%%%%%%%%%%%%%%%
function [units,convf]=convert_units(unt,origunits);
% Conversion factors from origianl units if requested and possible
% (no conversions from knots to feet).
%
  if strcmp(unt(1:3),origunits(1:3)) | strcmp(unt(1:3),'ori'),
    units=origunits;
    convf=1;
  else
   switch unt(1:3),
    case 'fee',
       if strcmp(origunits(1:3), 'met'),
  	  units='feet';
  	  convf=3.2808399;
       else
  	  units=origunits;
  	  convf=1;
       end;
    case 'met',
       if strcmp(origunits(1:3), 'fee'),
  	  units='meters';
  	  convf=0.3048;
       else
  	  units=origunits;
  	  convf=1;
       end;
    case 'm/s',
       if strcmp(origunits(1:3), 'kno'),
  	  units='meters/sec';
  	  convf=0.51444444;
       else
  	  units=origunits;
  	  convf=1;
       end;
    case 'kno',
       if strcmp(origunits(1:3), 'm/s'),
  	  units='knots';
  	  convf=1.9438445;
       else
  	  units=origunits;
  	  convf=1;
       end;
    otherwise
      error('Unknown units')
    end;
  end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [d,hdg]=t_gcdist(lat1,lon1,lat2,lon2)
% function [d,hdg]=t_gcdist(lat1,lon1,lat2,lon2)
% Function to calculate distance in kilometers and heading between two
% positions in latitude and longitude.
% Assumes -90 > lat > 90  and  -180 > long > 180
%    north and east are positive
% Uses law of cosines in spherical coordinates to calculate distance
% calculate conversion constants
%
%  Code from Richard Dewey.

raddeg=180/pi;
degrad=1/raddeg;
% convert latitude and longitude to radians
lat1=lat1.*degrad;
lat2=lat2.*degrad;
in1=find(lon1>180);lon1(in1)=lon1(in1)-360;
in2=find(lon2>180);lon2(in2)=lon2(in2)-360;
lon1=-lon1.*degrad;
lon2=-lon2.*degrad;
% calculate some basic functions
coslat1=cos(lat1);
sinlat1=sin(lat1);
coslat2=cos(lat2);
sinlat2=sin(lat2);
%calculate distance on unit sphere
dtmp=cos(lon1-lon2);
dtmp=sinlat1.*sinlat2 + coslat1.*coslat2.*dtmp;

% check for invalid values due to roundoff errors
in1=find(dtmp>1.0);dtmp(in1)=1.0;
in2=find(dtmp<-1.0);dtmp(in2)=-1.0;

% convert to meters for earth distance
ad = acos(dtmp);
d=(111.112) .* raddeg .* ad;

% now find heading
hdgcos = (sinlat2-sinlat1.*cos(ad))./(sin(ad).*coslat1);

% check value to be legal range
in1=find(hdgcos>1.0);hdgcos(in1)=1.0;
in2=find(hdgcos<-1.0);hdgcos(in2)=-1.0;
hdg = acos(hdgcos).*raddeg;

% if longitude is decreasing then heading is between 180 and 360
test = sin(lon2-lon1);
in1=find(test>0.0);
hdg(in1)=360-hdg(in1);