svpalm.m

GPS TOOLBOX包含以下内容： 1、GPS相关常量和转换因子； 2、角度变换； 3、坐标系转换： &#61656 点变换； &#61656 矩阵变换； &#61656 向量变换

%                                 svpalm.m
%  Scope:   This MATLAB macro computes ECEF satellite position based on 
%           satellite almanac data; WGS-84 constants are used.
%  Usage:   psa = svpalm(tsim,adata)
%  Description of parameters:
%           tsim      - input, tsim, GPS system time at time of transmission,
%                       i.e. GPS time corrected for transit time 
%                       (range/speed of light), in seconds
%           adata(1)  - input, toa, reference time almanac (time of 
%                       applicability), in seconds         
%           adata(2)  - input, smaxis (a), semi-major axis, in meters
%           adata(3)  - input, ecc (e), satellite eccentricity  
%           adata(4)  - input, izero (I_0), inclination angle at reference 
%                       time, in radians 
%           adata(5)  - input, razero (OMEGA_0), right ascension at reference 
%                       time, in radians (longitude of ascending node of orbit 
%                       plane at weekly epoch)
%           adata(6)  - input, argper (omega), argument of perigee, in radians
%           adata(7)  - input, mzero (M_0), mean anomaly at reference time, in 
%                       radians      
%           adata(8)  - input, radot (OMEGA_DOT), rate of right ascension, in 
%                       radians/second 
%           psa       - output, ECEF satellite position vector, the components 
%                       are in meters
%  External Matlab macros used:  wgs84con
%  Last update:  11/09/00
%  Copyright (C) 1996-00 by LL Consulting. All Rights Reserved.

function   psa = svpalm(tsim,adata) 

[nrow,ncol] = size(adata);
if  ncol ~= 8
   error('Error  -  SVPALM; check the dimension of the input data');
end

%  Constants and initialization

wgs84con;
% global constants used: gravpar, rot_rate

toa    = adata(1); 
smaxis = adata(2); 
ecc    = adata(3); 
izero  = adata(4); 
razero = adata(5); 
argper = adata(6); 
mzero  = adata(7); 
radot  = adata(8);

%  Compute GPS time of week at which SV position is computed

ttr = tsim;   
tk = ttr - toa;      %  time from almanac reference epoch, in seconds
if  tk > 302400.
  tk = tk - 604800;
elseif  tk < -302400.
  tk = tk + 604800;
end

%  Initialization 
          
mmot = sqrt(gravpar/smaxis/smaxis/smaxis);
r1mecc2 = sqrt(1. - ecc * ecc);
sinw = sin(argper);
cosw = cos(argper);
sinio = sin(izero);
cosio = cos(izero);

%  Compute mean anomaly

meanom = mzero + mmot * tk;

%  Compute eccentric anomaly  (Keppler's equation)

eccano = meanom + ecc * sin(meanom);                 %  initial guess 
for j = 1:6             % number of iterations can be changed if necessary
   temp1 = 1. - ecc * cos(eccano);
   temp2 = sin(eccano) - eccano * cos(eccano);
   eccano = (ecc * temp2 + meanom) / temp1;
end

%  Compute sine and cosine of eccentric anomaly, and rho

sine = sin(eccano);
cose = cos(eccano);
rho = 1. - ecc * cose;

%  Compute sine and cosine of true anomaly

cosv = (cose - ecc) / rho;
sinv = r1mecc2 * sine / rho;

%  Compute sine and cosine of argument of latitude

sinphi = cosv * sinw + cosw * sinv;
cosphi = cosv * cosw - sinv * sinw;

%  Correct argument of latitude  (no correction)

sinu = sinphi;
cosu = cosphi;

%  Compute satellite position in orbital plane

orbrad = smaxis * rho;
xprime = orbrad * cosu;
yprime = orbrad * sinu;

%  Correct inclination  (no correction)

sini = sinio;
cosi = cosio;

%  Compute longitude of ascending node, its sine and cosine

omega = razero + radot * tk - rot_rate * ttr;
sinome = sin(omega);
cosome = cos(omega);

%  Compute satellite position in ECEF

psa(1) = xprime * cosome - yprime * cosi * sinome;
psa(2) = xprime * sinome + yprime * cosi * cosome;
psa(3) = yprime * sini;