Orbital Coordinate Systems, Part II

By Dr. T.S. Kelso

November/December 1995

In our last column, we were working our way through the process of
calculating the position and velocity of an observer on the Earth's surface in
the Earth-Centered Inertial (ECI) reference frame. The goal, of course, was to
be able to determine the position of the observer relative to an orbiting
satellite to aid the process of visually acquiring or otherwise tracking the
satellite from the ground.

At the end of the last column, we had seen that it would be impossible to
determine an observer's position in the ECI reference frame without knowing that
observer's local sidereal time, that is, the angle between the observer's
meridian (longitude) and the vernal equinox. But because this time is measured
relative to the stars and not the Sun, the equations to do the calculations are
a bit complicated. Let's start of this column with a brief review of the
equation for the local sidereal time and then explore, bit by bit, how to
develop a computer routine to calculate this value for us. We will then continue
the process by calculating the remaining information needed to determine an
observer's position in the ECI reference frame.

As we saw last time, the local sidereal time, θ(τ), can be
calculated by adding the observer's east longitude, λE,
to the Greenwich sidereal time (GST), θg(τ). Using an
equation from Page 50 of the
Explanatory
Supplement to the Astronomical Almanac, we can first determine θg(0h),
the Greenwich sidereal time at 0h (midnight) UTC, from which

θg(τ) =
θg(0h) + ωe·Δτ

where Δτ is the UTC time of interest and ωe
= 7.29211510 × 10-5 radians/second is the
Earth's rotation rate. Recalling that θg(0h)
is given as

where Tu = du/36525 and
du is the number of days of Universal Time elapsed since JD
2451545.0 (2000 January 1, 12h UT1), we can now set about determing
what we know and what we need to calculate.

Let's develop a top-down algorithm to see what we need to do. Our goal in
this portion of the algorithm is to determine the value of θg
at a particular time, τ. To do this, we need values of
θg(0h), ωe, and Δτ. The
value of ωe is fixed and Δτ is the fraction of the
day since 0h UTC or Δτ = fraction(τ) in days.
Therefore, we still need to know Tu, which depends upon
du. But du depends upon knowing the number
of days of Universal Time elapsed since JD 2451545.0. So, if we know the Julian
Date (JD) of τ, we have everything we need. Our top-down algorithm then
looks something like this:

Therefore, our most basic calculation is the determination of the Julian Date. The
Julian Date can be calculated from the equation in Section 12.95 (Page 606)
[actually, this should be Section 12.92 (Page 604)] of the
Explanatory
Supplement to the Astronomical Almanac, or using the approach on Page 61 of
Astronomical
Algorithms by Jean Meeus. This latter text is an excellent reference source for
many relevant calculations in orbital mechanics and is highly recommended.

Using Meeus' approach, the Pascal code for calculating the Julian Date of January
0.0 of any year would be:

As an example, the Julian Date of 0h UTC on 1995 October 01 would be written
as:

JD := Julian_Date_of_Year(1995) + DOY(1995,10,1);

with a result of 2449991.5 = 2449717.5 + 274. Therefore, du
equals -1553.5 days and Tu equals -1553.5/36525. From this
value of Tu, θg(0h) can now be
calculated to be -343378.2154 seconds. Since there are 86,400 seconds in a day,
this is the same time (angle) as 2221.7846 seconds, so the equivalent GMST is
0h 37m 01s.7846 or an angle of 9.257
degrees.

Now, if the time of interest on 1995 October 01 was 9h UTC, we
would have to add ωe·Δτ to the value of
θg(0h), where Δτ = 32,400 seconds (9
hours). Being careful to use the proper units, our new GMST is 9h
38m 30s.4928 or 144.627 degrees.

We can consolidate our calculation of the Greenwich Mean Sidereal Time into
the following simple Pascal function where the input is the Julian Date of the
time of interest (in our example of 9h UTC on 1995 October 01, JD is
2449991.875) and the output is the GMST in radians. For our test case, this
should be 2.524218 radians.

Now, we can complete our calculation of the ECI position of our observer.
We'll start by assuming a spherical Earth and then go back and rework our
solution, in our next column, using an oblate Earth. If our observer is located
at 40° North latitude and 75° West longitude (near Philadelphia), we can
easily calculate the z coordinate according to

z = Re sin φ

where Re = 6378.135 km and φ = 40°. To
calculate x and y, we use

x = R cos θ

y = R sin θ

where

θ = θg + λE

R = Re cos φ.

Using our example, λE = -75° and θg =
144°.627, so θ = 69°.627 (remember, this is the local sidereal
time—the angle between the observer's meridian and the vernal equinox) and
R = 4885.935 km. Therefore, the ECI position of our observer at the time
of interest is:

x = 1700.938 km, y = 4580.302 km, z = 4099.786 km.

After all that work calculating the local sidereal time, the rest of the
calculation was relatively easy! We can encapsulate this calculation using the
following simple Pascal procedure:

In the inputs for this routine, lat and lon are in
radians, alt is in kilometers, and time is the Julian
Date of interest. The outputs (x, y, and
z) are in kilometers.

So, how do we now calculate the look angle to a satellite? If the satellite's
position in the ECI coordinate system is defined as [xs,
ys, zs] and the observer is
[xo, yo, zo],
then the range vector is simply

[rx, ry,
rz] = [xs - xo, ys
- yo, zs - zo].

This vector, however, is in the ECI system and to generate look angles, we
need it to be in the topocentric-horizon system shown in Figure 1. That system
has its z axis pointing toward the zenith, the x axis pointing
South, and the y axis pointing East.

Figure 1. Topocentric-Horizon Coordinate System

To transform to the topocentric-horizon system, we must first rotate through
an angle θ (the local sidereal time) about the z axis (Earth
rotation axis) and then through an angle φ (the observer's latitude) about
the y axis. The coordinates (rS, rE,
rZ) become:

rS = sin φ cos θ
rx + sin φ sin θ ry - cos φ
rz

rE = -sin θ rx
+ cos θ ry

rZ = cos φ cos θ
rx + cos φ sin θ ry + sin φ
rz

The range to the satellite is

r = √ [rS2
+ rE2 + rZ2],

the elevation is

El = sin-1(rZ / r),

and the azimuth is

Az = tan-1(-rE / rS).

The minus sign is necessary because azimuth is measured clockwise from North
instead of counter-clockwise from South (which would be standard for a
right-handed orthogonal coordinate system). Care must be taken with the azimuth
to ensure the proper quadrant is selected for the arctangent.

We can calculate the look angles using the procedure described with the
Pascal procedure Calculate_Look, shown below. The inputs are the
satellite's ECI coordinates (xs, ys, and
zs) in kilometers, the observer's latitude and longitude
(lat and lon) in radians and altitude
(alt) in kilometers, along with the time of interest
(time) as a Julian Date. The outputs are the azimuth and elevation
(az and el) in radians and the range (rg)
in kilometers.

To sum things up, in order to calculate the look angles for a satellite
relative to an observer on the ground, we must first calculate the satellite's
position in the ECI coordinate system, then calculate the observer's ECI
position, take the difference of the two vectors, and then transform (rotate)
the vector from the ECI coordinate frame to the topocentric-horizon frame. The
most difficult part of this process is in calculating the Earth's rotation angle
when determining the observer's position.

As always, if you have questions or comments on this column, feel free to
send me e-mail at
TS.Kelso@celestrak.com or write care of
Satellite Times. Until next time, keep looking up!