EAS-591T – Space Geodetic Measurements of Active Crustal
Motions
LAB 4
From GPS ephemerides to ECEF satellite positions
GPS satellite ephemerides are broadcast to users as part of the GPS signal in the navigation message. They are calculated by predicting orbits forward on the basis of a curve-fit to 4 to 6 hours of actual orbit data. GPS broadcast ephemerides are accurate to ~10 m. They include the following parameters:
- toe: epoch of reference of the ephemerides in seconds of the current GPS week. It is the center of the interval over which the ephemeris is valid.
- m0, Dn, e, sqrt(a), W0, i0, w0: orbital elements (= Keplerian parameters) at toe.
- Wdot, idot: rate of right ascension and inclination angle.
- Cwc Cws Crc Crs Cic Cis: corrections coefficients to the argument of perigee, orbit radius and inclination for the perturbations caused by variations in the Earth’s gravity field, solar radiation pressure, and attraction from Sun and Moon.
Assignment:
1.
Write a program
to convert the ephemerides given in a RINEX navigation file into an ECEF
coordinate system. Input file = epgga2.010:
·
Write a
function to read the broadcast ephemeris file. This function should return a
matrix with 17 rows (the parameters needed later on in the calculations) and as
many columns as there are ephemerides. You will need the 17 following
parameters:
|
svprn |
satellite PRN
number |
|
Mo |
mean anomaly |
|
roota |
sqrt(semi-major
axis) |
|
deltan |
variation of mean
angular velocity |
|
ecc |
eccentricity |
|
omega0 |
argument of
perigee |
|
cuc, cus, crc, crs, cic,
cis |
correction
coefficients |
|
i0 |
inclination |
|
idot |
rate of inclination |
|
Omega0 |
right ascension |
|
Omegadot |
rate of right
ascension |
|
toe |
time of ephemeris |
The call to
that function could look like:
eph = read_rinexn('epgga2.010');
·
Write a
function to convert the ephemerides matrix (output of previous function) for a
given satellite at a given time (i.e. 3 input arguments). This function should
return the corresponding X,Y,Z coordinates in ECEF frame, for instance:
[X,Y,Z] = get_satpos(t,sv,eph);
A possible algorithm is given below, together
with the formulas that give the values for W, w, i, E, r, and v
needed in the equations above:
1. Extract needed
parameters from the GPS broadcast ephemerides file. Note that the RINEX navigation files provide information every 2 hours
only. Therefore, you first need to find the correct time of ephemeris (toe)
in the navigation file given the time for which you need to extrapolate the satellite
position (tdata).
2. Compute
basic parameters at request time tdata:
Time elapsed
since toe: 
Mean anomaly
at t: 
Iterative
solution for E: 
(use iterative solution)
True anomaly v: 
3. Correct for orbital perturbations:
Argument of
perigee: 
Radial distance: 
Inclination: 
4. Compute the right
ascension, accounting for
Earth's rotation (we = mean angular velocity of the Earth) and
variations of W (= Wdot):
5. Convert satellite position from orbital frame to ECEF frame:
Write the satellite position in the orbital
frame: 
Build the
rotation matrix to ECEF frame:
Apply rotation:
6. Check your program.
Use values from the sp3 file igs11484.sp3 as “ground truth”:
t=346500 and sv = 31 Ž X=11660.379642, Y=11313.211213, Z=-21326.822815
2.
Compute ECEF
coordinates for satellite 31 every 15 minutes and compare with the coordinates
given in the corresponding sp3 file (precise IGS orbits).
·
Write a
function that extract the XYZ position of a satellite ‘sv’ from an sp3 file and
returns that position:
[Xs,Ys,Zs,Ts] =
read_sp3(sp3file,sv);
·
Compute the 3-D
residual position = vector difference between sp3 and broadcast XYZ position.
Plot the norm of this vector as a function of time for satellite 31.
3.
Comment on the
differences between the 2 orbits.
Short description of the RINEX navigation file format. The format of fully described in http://www.ngs.noaa.gov/CORS/Rinex2.html.
+----------------------------------------------------------------------------+
| TABLE
A4 |
| GPS NAVIGATION MESSAGE FILE - DATA RECORD DESCRIPTION |
+--------------------+------------------------------------------+------------+
| OBS. RECORD | DESCRIPTION |
FORMAT |
+--------------------+------------------------------------------+------------+
|PRN / EPOCH
/ SV CLK| - Satellite PRN number |
I2, |
| | - Epoch: Toc - Time of
Clock | |
| | year (2 digits, padded with 0 | |
| | if necessary) | 1X,I2.2, |
| | month |
1X,I2, |
| | day | 1X,I2, |
| | hour |
1X,I2, |
| | minute |
1X,I2, |
| | second |
F5.1, |
| | - SV clock bias (seconds) | 3D19.12 |
| | - SV clock drift (sec/sec) | |
| | - SV clock drift rate
(sec/sec2) | |
+--------------------+------------------------------------------+------------+
| BROADCAST
ORBIT - 1| - IODE Issue of Data, Ephemeris | 3X,4D19.12 |
| | - Crs (meters) | |
| | - Delta n (radians/sec) | |
| | - M0 (radians) | |
+--------------------+------------------------------------------+------------+
| BROADCAST
ORBIT - 2| - Cuc (radians) | 3X,4D19.12 |
| | - e Eccentricity | |
| | - Cus (radians) | |
| | - sqrt(A) (sqrt(m)) | |
+--------------------+------------------------------------------+------------+
| BROADCAST
ORBIT - 3| - Toe Time of Ephemeris | 3X,4D19.12 |
| | (sec of GPS week) |
|
| | - Cic (radians) | |
| | - OMEGA (radians) | |
| | - CIS (radians) | |
+--------------------+------------------------------------------+------------+
| BROADCAST
ORBIT - 4| - i0
(radians) | 3X,4D19.12 |
| | - Crc (meters) | |
| | - omega (radians) | |
| | - OMEGA DOT (radians/sec) | |
+--------------------+------------------------------------------+------------+
| BROADCAST
ORBIT - 5| - IDOT (radians/sec) | 3X,4D19.12 |
| | - Codes on L2
channel | |
| | - GPS Week # (to go with
TOE) | |
| | Continuous number, not mod(1024)! | |
| | - L2 P data flag | |
+--------------------+------------------------------------------+------------+
| BROADCAST
ORBIT - 6| - SV accuracy
(meters) | 3X,4D19.12 |
| | - SV health (bits 17-22 w 3 sf 1) | |
| | - TGD (seconds) | |
| | - IODC Issue of Data,
Clock | |
+--------------------+------------------------------------------+------------+
| BROADCAST
ORBIT - 7| - Transmission time of message *) | 3X,4D19.12 |
| | (sec of GPS week, derived e.g. | |
| | from Z-count in Hand Over Word (HOW) |
|
| | - Fit interval (hours) | |
| | (see ICD-GPS-200, 20.3.4.4) | |
| | Zero if not known |
|
| | - spare | |
| | - spare | |
+--------------------+------------------------------------------+------------+