More on Orbits, Part 5

One more quantity we need is the orbital inclination i which is found in the almanac.

Now taking the x' and y' coordinates, and angles Ω and i, we get earth-centered cartesian coordinates:

x = x' * cos Ω - y' * cos i * sin Ω
y = x' * sin Ω + y' * cos i * cos Ω
z = y' * sin i

Finally, the satellite longitude is obtained by the arctangent of y/x. We have to be careful of which quadrant we are in. Many programming languages have a handy function, atan2(x, y), that returns the angle in the proper quadrant. Note that some languages have the arguments in the other order, (y, x).

lon = atan2(x, y)

Latitude can be found with arctangent of z / sqrt(x² + y²):

lat = atan2( sqrt(x² + y²), z)

Now with a (lat, lon) ordered pair for the satellite coordinates, we can finally get the elevation angle above the horizon from our location (base_lat, base_lon). This involves some trigonometry. We need the angle gamma, which is between a line from the center of the earth to our base location, and a line from the center of the earth to the satellite. It is given by:

cos(gamma) = cos(base_lat) * cos(lat) * cos(lon - base_lon) + sin(base_lat) * sin(lat)

With this we can get the satellite range, the distance d from our base location:

d = r * sqrt(1 + Re² / r² - 2 * Re / r * cos(gamma) )

where r is the satellite radius that was calculated earlier; and Re is radius of the earth, 6370 km.

Then cosine of the elevation cos(El) = r * sin(gamma) / d

One of the complications is that the given equations work well if the satellite is on the same side of the earth. There has to be a way to check if cos(El) is actually valid. A negative elevation will still give a positive cosine, and the arccos will return a positive elevation. Using the law of cosines, we can get another angle, psi, that is between a line from the center of the earth to our station, and a line from our station to the satellite. Using the law of cosines:

cos(psi) = (d² + Re² - r²) / (2 * d * Re)

The elevation will be psi minus 90 degrees. So if psi is < 90 degrees then the elevation is negative. If the cosine is positive, then psi < 90 degrees and we want to disregard satellites with such an angle. They are below the horizon.