Moved to Wordpress

I decided to move to Wordpress.

jasonlivermore.wordpress.com

More on Orbits, Part 6

All this effort should produce a list of satellites with reasonably good accuracy. There are some much more sophisticated software models out there that take into account other perturbations like atmospheric drag, the pull of the sun and moon, and so on. I have ignored these for simplicity. I used Python to code this up, because it's such a cool scripting language and has so many modules. Then I thought about plotting the satellite positions on Google Maps, and started looking into the API. There were some handy examples that I borrowed from, to generate an HTML file that I could then load. Of course, my thoughts turned to hosting it, and it turns out that Google will host Python and Java apps with its App Engine. Some more digging up of examples and documentation led to this: http://gpslookangle.appspot.com I'll put more info into the "About" link on that page.

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.

More on Orbits, Part 4

With our polar coordinates (r, ν) in hand, there are a few more transformations.

The argument of latitude u is obtained by adding the argument of perigee ω to ν. Argument of perigee (or more generally, periapsis) is taken from the GPS almanac.

u = ν + ω

Next we convert to cartesian coordinates to get position in the orbital plane:

x' = r * cos(u)
y' = r * sin(u)

Now we need to account for rotation of the earth. The right ascension of the ascending node Ω₀ is given in the almanac, and it is valid at the beginning of the referenced GPS week, not at the time of applicability. So we add the rotation speed of the earth times the number of seconds since the start of the week.

Ω = Ω₀ - Ω_e(dot) * t

This will rotate our coordinates so we can get a longitude of the satellite relative to earth longitude. The quantity Ω_e(dot) is the rotation rate of the earth, 7.292115 x 10-5 radians per second.

With all these new quantities in hand, we'll use a rotation matrix to move the orbit into earth coordinates: latitude and longitude.

More on Orbits, Part 3

So far we haven't accounted for the eccentricity of the orbit. We can get this number from the appropriate line in the almanac file. A value of zero means the orbit is circular. As the value moves towards 1, the orbit is more and more enlongated. GPS satellites have very circular orbits, with values less than 0.02. Comets typically have values much closer to 1.

From mean anomaly, we find the eccentric anomaly E, from mean anomaly M that we found in the previous post, and eccentricity ε.

M = E - ε sin E

This equation needs to be solved iteratively. Newton's method can be used to get a very accurate value after a few loops, by finding the root of the equation:

f(E) = E - ε sin E - M = 0

The derivative with respect to E is f '(E) , or: 1 - ε cos E

We can start with a first guess E₀ = M. Take E₁ = E₀ - f(E₀)/f '(E₀)

So E₁ = E₀ - (E₀ - ε sin E₀ - M) / (1 - ε cos E₀)

Then find E₂ with E₁ in the above function, and so on. You can check E_n in the original equation and keep going until E_n causes M to converge to the correct value.

E gives us the angle from the center of the ellipse to the satellite. Note that the center is the point between the two foci of the ellipse. The earth will be at one focus, not at the center.

At this point we can get the radius r which represents the distance of the satellite from the center of the earth.

r = a * (1 - ε * cos E)

Finally, to get the angle of the satellite position relative to earth, we convert to true anomaly ν.

cos ν = (cos E - ε) / (1 - ε * cos E)

Now we have the position in polar coordinates (r, ν) of the satellite in the orbital plane. There is still a little more work to get a current latitude and longitude position.

More on Orbits, Part 2

I should add a note here about how to identify the satellites. All the satellites transmit on the same frequency, and each satellite has a unique pseudo-random noise (PRN) code that spreads its message data over a wide bandwidth. This is a direct sequence spread spectrum method known as CDMA. There are 32 unique codes so there will be up to 32 active satellites, and they are identified in the almanac by their PRN codes.

The GPS almanac references its time relative to the start of a GPS week. This week starts at midnight UTC time on the first day of the calendar week (Saturday night turning into Sunday morning.) The orbital parameters given in the almanac are valid at the reference time, or epoch, which is given as "Time of Applicability". This time is the number of seconds after the start of the referenced GPS week.

The satellite position is given as mean anomaly at the epoch, on the line in the almanac file that starts with "Mean Anom". Mean anomaly is an angle measurement that increases linearly over time. We get the current mean anomaly using the difference between our time of interest t and the epoch t_oe:

t_k = t - t_oe

Then the current mean anomaly is

M = M₀ + n * t_k

There is our mean motion n from the previous post. Now we have an angle measurement that will help us find the position of the satellite in its orbit at time t.

More on Orbits

Awhile back, I wrote about my orbit simulator. This has been on my mind again, as I'm dealing with some issues at work that have to do with GPS receivers. GPS is useful to synchronize widely dispersed systems such as cellular base stations. It not only provides location, but also very accurate timing. In fact, it is the timing that allows you to triangulate your exact position.

My concern is in predicting which satellites are in view at any given time. There are some websites that track satellite orbits and predict passes over your location, but I'm looking for something more specific to what I want. How to go about this? It turns out that a lot of the same equations I used in my orbit simulator come into play for predicting locations of real world satellites.

The first good reference is the ICD-GPS-200C document which is pretty easy to find on the web. The equations starting in Table 20-IV are what we need.

Since the satellite orbits are perturbed by various influences, the GPS almanac is updated daily. The orbits are generally stable, but will have some small variance from day to day. Almanacs are available here.

The first equation I'll deal with is the mean motion, n.

n = sqrt(GMe / a³)

GMe is 3.986 x 10¹⁴ and relates to the gravitational constant and the mass of the earth. The semi-major axis a has to do with the size of the elliptical orbit. The almanac gives the square root of a in the line that starts with SQRT(A).

The mean motion relates to the inverse of the orbital period and gives the number of radians per second that a given satellite is moving around its orbit. This will be important for finding the position at a given time, because the almanac just gives a snapshot of where it was, or will be, at some reference time.

That is all for this post. The next step is finding the position in orbit at a particular time.