Orbital Coordinate Systems, Part IIBy Dr. T.S. Kelso 



Orbital Coordinate Systems, Part IIBy Dr. T.S. Kelso 
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 EarthCentered 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}(0^{h}), the Greenwich sidereal time at 0^{h} (midnight) UTC, from which
θ_{g}(τ) = θ_{g}(0^{h}) + ω_{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}(0^{h}) is given as
θ_{g}(0^{h}) = 24110^{s}.54841 + 8640184^{s}.812866 T_{u} + 0^{s}.093104 T_{u}^{2}  6.2 × 10^{6} T_{u}^{3}
where T_{u} = d_{u}/36525 and d_{u} is the number of days of Universal Time elapsed since JD 2451545.0 (2000 January 1, 12^{h} UT1), we can now set about determing what we know and what we need to calculate.
Let's develop a topdown 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}(0^{h}), ω_{e}, and Δτ. The value of ω_{e} is fixed and Δτ is the fraction of the day since 0^{h} UTC or Δτ = fraction(τ) in days. Therefore, we still need to know T_{u}, which depends upon d_{u}. But d_{u} 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 topdown algorithm then looks something like this:
Calculate: θ_{g} which depends on (θ_{g}(0^{h}), ω_{e}, Δτ); ω_{e} constant; Δτ = fraction(τ)
Calculate: θ_{g}(0^{h}) which depends on (T_{u}); T_{u} = d_{u}/36525
Calculate: d_{u} which depends on (JD(τ))
Calculate: JD(τ)
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:
Function Julian_Date_of_Year(year : double) : double; { Calculate Julian Date of 0.0 Jan year } var A,B : longint; begin year := year  1; A := Trunc(year/100); B := 2  A + Trunc(A/4); Julian_Date_of_Year := Trunc(365.25 * year) + Trunc(30.6001 * 14) + 1720994.5 + B; end; {Function Julian_Date_of_Year}
To calculate the Julian Date of any calendar date, we simply combine the Julian Date of that year with the day of the year, where the day of the year can be calculated as:
Function DOY(yr,mo,dy : word) : word; const days : array [1..12] of word = (31,28,31,30,31,30,31,31,30,31,30,31); var i,day : word; begin day := 0; for i := 1 to mo1 do day := day + days[i]; day := day + dy; if ((yr mod 4) = 0) and (((yr mod 100) <> 0) or ((yr mod 400) = 0)) and (mo > 2) then day := day + 1; DOY := day; end; {Function DOY}
As an example, the Julian Date of 0^{h} 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, d_{u} equals 1553.5 days and T_{u} equals 1553.5/36525. From this value of T_{u}, θ_{g}(0^{h}) 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 0^{h} 37^{m} 01^{s}.7846 or an angle of 9.257 degrees.
Now, if the time of interest on 1995 October 01 was 9^{h} UTC, we would have to add ω_{e}·Δτ to the value of θ_{g}(0^{h}), where Δτ = 32,400 seconds (9 hours). Being careful to use the proper units, our new GMST is 9^{h} 38^{m} 30^{s}.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 9^{h} 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.
Function ThetaG_JD(jd : double) : double; { Reference: The 1992 Astronomical Almanac, page B6. } var UT,TU,GMST : double; begin UT := Frac(jd + 0.5); jd := jd  UT; TU := (jd  2451545.0)/36525; GMST := 24110.54841 + TU * (8640184.812866 + TU * (0.093104  TU * 6.2E6)); GMST := Modulus(GMST + 86400.0*1.00273790934*UT,86400.0); ThetaG_JD := twopi * GMST/86400.0; end; {Function ThetaG_JD}
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 = R_{e} sin φ
where R_{e} = 6378.135 km and φ = 40°. To calculate x and y, we use
x = R cos θ
y = R sin θ
where
θ = θ_{g} + λ_{E}
R = R_{e} 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:
Procedure Calculate_User_Pos(lat,lon,alt,time : double; var x,y,z : double); { Reference: The 1992 Astronomical Almanac, page K11. } const re = 6378.135; var theta : double; begin theta := Modulus(ThetaG_JD(time) + lon,twopi); r := (re + alt)*Cos(lat); x := r*Cos(theta); y := r*Sin(theta); z := (re + alt)*Sin(lat); end; {Procedure Calculate_User_Pos}
In the inputs for this routine,
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 [x_{s}, y_{s}, z_{s}] and the observer is [x_{o}, y_{o}, z_{o}], then the range vector is simply
[r_{x}, r_{y}, r_{z}] = [x_{s}  x_{o}, y_{s}  y_{o}, z_{s}  z_{o}].
This vector, however, is in the ECI system and to generate look angles, we need it to be in the topocentrichorizon 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. TopocentricHorizon Coordinate System
To transform to the topocentrichorizon 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 (r_{S}, r_{E}, r_{Z}) become:
r_{S} = sin φ cos θ r_{x} + sin φ sin θ r_{y}  cos φ r_{z}
r_{E} = sin θ r_{x} + cos θ r_{y}
r_{Z} = cos φ cos θ r_{x} + cos φ sin θ r_{y} + sin φ r_{z}
The range to the satellite is
r = √ [r_{S}^{2} + r_{E}^{2} + r_{Z}^{2}],
the elevation is
El = sin^{1}(r_{Z} / r),
and the azimuth is
Az = tan^{1}(r_{E} / r_{S}).
The minus sign is necessary because azimuth is measured clockwise from North instead of counterclockwise from South (which would be standard for a righthanded 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
Procedure Calculate_Look(xs,ys,zs,lat,lon,alt,time : double; var az,el,rg : double); var xo,yo,zo, rx,ry,rz, theta, top_s,top_e,top_z : double; begin Calculate_User_Pos(lat,lon,alt,time,xo,yo,zo); theta := Modulus(ThetaG_JD(time) + lon,twopi); rx := xs  xo; ry := ys  yo; rz := zs  zo; top_s := Sin(lat)* Cos(theta)*rx + Sin(lat)* Sin(theta)*ry  Cos(lat)*rz top_e :=  Sin(theta)*rx + Cos(theta)*ry; top_z := Cos(lat)* Cos(theta)*rx + Cos(lat)* Sin(theta)*ry + Sin(lat)*rz; az := ArcTan(top_e/top_s); if top_s > 0 then az := az + pi; if az < 0 then az := az + twopi; rg := Sqrt(rx*rx + ry*ry + rz*rz); el := ArcSin(top_z/rg); end; {Procedure Calculate_Look}
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 topocentrichorizon 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 email at TS.Kelso@celestrak.com or write care of Satellite Times. Until next time, keep looking up!