eph2xyzn.m

这是国外关于卫星导航方面一书的源代码

function [xyzx,xyzy,xyzz,delta_ts,rela]=eph2xyzn(t,pr,rinexdat,wie)
%---------------------------------------------------------------------------------------
% GPSLab:   EPH2XYZN.M
%
%           WGS84 Coordinates of GPS Satellite Orbits from RINEX N-Files
%
% -----------------------------------------------------------------------------
% WHAT: Epochwise evaluation of XYZ(Sat) from Kepler-elements (Sat)
%
% HOW: [x,y,z,delta_ts,rela]=eph2xyzn(t,pr,[svprn,toc,a0,a1,a2,crs,dn,ma,cuc,ecc,cus, ...
%                            ... sqrta,toe,cic,omega0,cis,oinc,crc,aop,omegadot,idot, ...
%                            ... gweek,accu,health])
%
% INPUT: t,pr,[svprn,toc,a0,a1,a2,crs,dn,ma,cuc,ecc,cus,sqrta,toe,cic,omega0,cis, ...
%        ... oinc,crc,aop,omegadot,idot,gweek,accu,health]
%        
%        All angles have to be in "radians" (RINEX convention), NOT in "semi-circles"
%	      Time unit is "seconds" & linear unit is "meters" (RINEX convention)
%
%        t - time of epoch [seconds of the GPS week]
%
%        pr - pseudorange [meter] as given in RINEX observation file for epoch t
%        further input sequence in RINEX nav file compatible order (nav blocks after header),
%        but kill / replace the following values:
%                   year,month,day,h,m,s (2nd to 7th value of 1st row) 
%                      - replaced by toc in [sec of GPS week]
%                   IODE (1st value of 2nd row) - killed
%                   Codes,L2PFlag (2nd and 4th value of 6th row) - killed
%                   TGD,IODC (3rd and 4th value of 7th row) - killed
%                   MessageTransmissionTime (1st value of 8th row) & 3 spares - killed
%
%        wie - flag for decision between inertial system coordinates or ECEF coordinates:
%              wie==0 means inertial, wie==1 means ECEF
%
% OUTPUT: X(Sat), Y(Sat), Z(Sat), satellite clock correction for actual epoch in [sec],
%         relativistic correction for periodic effect due to orbit eccentricity [meter]
%
% Attention please: 1) Valid only for the assumption  toc = toe !!!  (simplification)
%                      toc = ref.time of clock   toe = ref.time of ephemeris
%                   2) Input/Output not vectorized
%                   3) No effective programming, for educational purposes only!!
%                      --> normally Kepler's equation and iteration should be
%                          extracted and put into own functions
%
% 1996-05-10, Zebhauser (besides literature the algorithm was partially adopted from or completed by 
%                        comparisons with code by Peter Dana from UniTex and DIPOP from Uni New Brunswick
% 1998-03-11, Zebhauser (Update to actual Eph format including toc,gweek,accu,health)
% 1999-02-19, Zebhauser (choosing possibility between inertial and ECEF added)
% 1999-09-13, Zebhauser (Check of assumption, that toc=toe)
% 1999-10-11, Zebhauser (Removing substraction of relativist. corr. related to TOC from a0,a1,a2)
%
% Literature: GPS SPS Signal Specification - 2nd edition (1995) S.38, ICD-GPS-200(IRN-200C-002) (1997) S.98-100,
%             Parkinson, Spilker, Axelrad, Enge (eds.) (1996/1+2), Leick (1995), W黚bena (1991), Seeber (1989), ...
%---------------------------------------------------------------------------------------

if nargin==3
   wie=1;
end

if wie~=1&wie~=0
   errordlg('so ned! Es gibt nur die Optionen 0 oder 1 f黵 die zweite Variable (wird =1 gesetzt).', ...
      'GPSLab: Fehler beim Aufruf von EPH2XYZN.M');
   wie=1;
end

svprn=rinexdat(1);
toc=rinexdat(2);        % noch einbauen! Aenderung bei Berechnung der Sat.uhrparameter
a0=rinexdat(3);
a1=rinexdat(4);
a2=rinexdat(5);
crs=rinexdat(6);
dn=rinexdat(7);
ma=rinexdat(8);
cuc=rinexdat(9);
ecc=rinexdat(10);
cus=rinexdat(11);
sqrta=rinexdat(12);
toe=rinexdat(13);
cic=rinexdat(14);
omega0=rinexdat(15);
cis=rinexdat(16);
oinc=rinexdat(17);
crc=rinexdat(18);
aop=rinexdat(19);
omegadot=rinexdat(20);
idot=rinexdat(21);
gweek=rinexdat(22);    % noch Abfrage in Verbindg. mit WEEK/HALFWEEK einbauen (Wochenuebergang)
accu=rinexdat(23);     % eher fuer Kontrollzwecke im Hauptprogramm (Warnung, wenn > 10 m)
health=rinexdat(24);   % eher fuer Kontrollzwecke im Hauptprogramm (Bei Abwahl von unhealthy Sats)
%---------------------------------------------------------------------------------------
if toc~=toe
   gpsjpg=imread('sat.jpg','JPEG');
   svstr=num2str(svprn);
   hmsgeph=msgbox(['Achtung - Kollision: Bedingung TOC=TOE nicht erf黮lt !!!' ...
' Referenzzeit f黵 Ephemeriden ist nicht identisch mit Referenzzeit f黵 Uhrparameter' ...
' bei dem Satelliten mit der PRN-Nummer ',svstr, ...
'. Satellitendaten per Hand aus der Ephemeridendatei l鰏chen oder andere Daten laden und neu starten.'], ...
'GPSLab Abbruch: Nicht abgefangener Sonderfall','custom',gpsjpg,summer(64));
clear gpsjpg;
return;
end

% 躡erpr黤ung der GPS-Woche in Ephemeriden und Header (s1head, s2head) auf Identit鋞 entf鋖lt.
% M黶ste in den aufrufenden Programmen geschehen.

%---------------------------------------------------------------------------------------
HALFWEEK=302400;			% half GPS week in sec
WEEK=604800;				% full -"-
GM=3.986004415E14;				% grav. const. (WGS84)
OMEGA_DOT_E=7.2921151467E-5;		% earth rot. rate (WGS84)
if wie==0
   OMEGA_DOT_E=0;
end
C=2.99792458E8;				% speed of light (WGS84)
% F=-2*GM^.5/C^2;                       % const. term of Kepler-part of the relativist. delay
F=-4.442807633E-10;
%---------------------------------------------------------------------------------------
n0=sqrt(GM/sqrta^6); 			% computed mean motion
n=n0+dn; 		                % corrected mean motion
rel=F*sqrta*ecc;                        % factor for relativist. correction, period. part cause of ecc.
%---------------------------------------------------------------------------------------

x=ma; 					% kepler's equation for eccentric anomaly ek ...
y=ma-(x-ecc*sin(x));                    % ... starting with  ek=ma (mean anomaly)
x1=x;
x=y;
for i=0:15				% determination of ek by iterating
x2=x1;
x1=x;
y1=y;
y=ma-(x-ecc*sin(x));
if(abs(y-y1)<1.0E-15)
break;
end
x=(x2*y-x*y1)/(y-y1);
end					% end of iterations
ek=x; 					% end of det. of ecc. anomaly
sinek=sin(ek);
cosek=cos(ek);
%---------------------------------------------------------------------------------------
u0=a0;     % substraction of relativist. corr. related to TOC from a0,a1,a2: cancelled !!!!!
u1=a1;     % this substraction was proposed by van Dierendonck (1978), but is not valid
u2=a2;     %  really. Could be that it was valid in the beginning of GPS era. Actual data
%             show that a0 is free of the relativistic effect due to eccentricity, and the
%             ICD200c notifies that (a little bit loosely, but does it ! )  Zeb, 1999-10-11
%---------------------------------------------------------------------------------------
% til here only variables related to toe,
% from now on computation of transmission time since toe
%---------------------------------------------------------------------------------------
t0s=t-pr/C;                             % time of transmission
tk=t0s-toe;                             % time of transmission since toe (time of ephemeris)
if (tk>HALFWEEK)
tk=tk-WEEK;
end
if(tk<-HALFWEEK)
tk=tk+WEEK;
end

mk=ma+n*tk; 				% mean anomaly

x=mk; 					% kepler's equation for eccentric anomaly ek
y=mk-(x-ecc*sin(x));
x1=x;
x=y;
for i=0:15				% determination of ek by iterating
x2=x1;
x1=x;
y1=y;
y=mk-(x-ecc*sin(x));
if(abs(y-y1)<1.0E-15)
break;
end
x=(x2*y-x*y1)/(y-y1);
end					% end of iterations
ek=x; 					% end of det. of ecc. anomaly
sinek=sin(ek);
cosek=cos(ek);

delta_ts=u0+u1*tk+u2*tk^2+rel*sinek;    % sat.clock corr.: offset,drift,ageing & relativist. corr.

tk=tk-delta_ts;                         % time of transmission since toe with clock and rel. corr.
if (tk>HALFWEEK)
tk=tk-WEEK;
end
if(tk<-HALFWEEK)
tk=tk+WEEK;
end

%---------------------------------------------------------------------------------------
% from now on using this transmission time since toe (corr.)
%---------------------------------------------------------------------------------------

mk=ma+n*tk; 				% mean anomaly
x=mk; 					% kepler's equation for eccentric anomaly ek
y=mk-(x-ecc*sin(x));
x1=x;
x=y;
for i=0:15				% determination of ek by iterating
x2=x1;
x1=x;
y1=y;
y=mk-(x-ecc*sin(x));
if(abs(y-y1)<1.0E-15)
break;
end
x=(x2*y-x*y1)/(y-y1);
end					% end of iterations
ek=x; 					% end of det. of ecc. anomaly

cosek=cos(ek); 				
sinek=sin(ek);
%---------------------------------------------------------------------------------------
cosvk=(cosek-ecc); 			% denominator for det. of true anomaly
sinvk=sqrt(1-ecc*ecc)*sinek; 		% numerator for det. of true anomaly
vk=atan2(sinvk,cosvk); 			% true anomaly
if (vk<0)
vk=vk+2*pi;
end
%---------------------------------------------------------------------------------------

pk=vk+aop; 				% argument of latitude

sin2pk=sin(2.0*pk);
cos2pk=cos(2.0*pk);
duk=cus*sin2pk+cuc*cos2pk;		% argument of latitude correction
drk=crc*cos2pk+crs*sin2pk; 		% radius correction
dik=cic*cos2pk+cis*sin2pk;		% correction to inclination

uk=pk+duk; 				% latitude
rk=sqrta*sqrta*(1-ecc*cosek)+drk;	% corrected radius
ik=oinc+dik+idot*tk;   			% corrected inclination
xkp=rk*cos(uk); 			% x in orbital plane
ykp=rk*sin(uk); 			% y in orbital plane
%---------------------------------------------------------------------------------------
ok=omega0+(omegadot-OMEGA_DOT_E)*tk-OMEGA_DOT_E*toe;  	% longitude of ascending node

ykpcosik=ykp*cos(ik);
sinok=sin(ok);
cosok=cos(ok);
%---------------------------------------------------------------------------------------
xyzx=xkp*cosok-ykpcosik*sinok;		% ecef sv coordinates x,y,z
xyzy=xkp*sinok+ykpcosik*cosok;
xyzz=ykp*sin(ik);
%---------------------------------------------------------------------------------------
% sv relativistic correction in meters, periodic effect because of the eccentricity
% of the elliptic orbit with max. range correction of 6.79 m
% (constant part/ secular drift corrected in sv clocks by adding -0.445 to the nominal frequency)
% integrated in a0,a1,a2 of clock parameters, related to toc

rela=C*F*sqrta*ecc*sinek;

%---------------------------------------------------------------------------------------
% ENDE
% 1996-05-10, Zebhauser
% 1998-03-11, Zebhauser (Update to actual Eph format including toc,gweek,accu,health)
% 1999-02-19, Zebhauser (choosing possibility between inertial and ECEF added)
%---------------------------------------------------------------------------------------