www.pudn.com > GPS_satellite_positions.rar > 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übbena (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ür 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üllt !!!' ...
' Referenzzeit für Ephemeriden ist nicht identisch mit Referenzzeit für Uhrparameter' ...
' bei dem Satelliten mit der PRN-Nummer ',svstr, ...
'. Satellitendaten per Hand aus der Ephemeridendatei löschen oder andere Daten laden und neu starten.'], ...
'GPSLab Abbruch: Nicht abgefangener Sonderfall','custom',gpsjpg,summer(64));
clear gpsjpg;
return;
end
% Überprüfung der GPS-Woche in Ephemeriden und Header (s1head, s2head) auf Identität entfällt.
% Müsste in den aufrufenden Programmen geschehen.
%---------------------------------------------------------------------------------------
HALFWEEK=302400; % half GPS week in sec
WEEK=604800; % full -"-
GM=3.986004418E14; % 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)
%---------------------------------------------------------------------------------------