Thank you to Charlie Kluepfel <ChasKlu@aol.com> for providing the following information:
Here is the commented program for solving the spherical triangle with the north pole at the top, the northernmost point of the satellite's orbit (coincidental with NYC at the time the satellite passed) at the bottom (south end) of the left (west) side, and the current subsatellite point making the third point of the triangle. The triangle is a right triangle (the orbit is perpendicular to the meridian at its northernmost point), thereby simplifying the spherical law of cosines. The latitude of NYC is taken as 41, which is correct to the nearest degree. Aside from the trigonometry, conceptually, the satellite travels at a constant 15 degrees per hour along its orbit. At the northernmost and southernmost points it is traversing lines of longitude that are closer together than they are at the equator and therefor it covers more than 15 degrees of longitude per hour, while the earth itself goes through only 15 degrees per hour (at the north pole it doesn't move at all, at the equator the surface goes at 1000 mph, at NYC, 755 mph, considering velocities). The opposite effect occurs on each of the two passages over the equator, where part of the satellite's 15-degree-per-hour motion is a north-south component, and only part is an east-west component, while the meridians are fully separated at 15 degrees of longitude per 15 degrees of great circle, so the satellite's east-west motion lags the earth's surface's. The result is a figure-8.
The program is in QuickBasic:
'set up trig functions in degrees:
pi = 4 * ATN(1)
dr = pi / 180
DEF fnsn (x) = SIN(dr * x)
DEF fncs (x) = COS(dr * x)
DEF fnas (x) = ATN(x / SQR(1 - x * x)) / dr
DEF fnac (x) = 90 - fnas(x)
OPEN "analemma" FOR OUTPUT AS #1
sin41 = fnsn(41): cos41 = fncs(41) 'NYC taken as at 41 degrees north lat
FOR i = 0 TO 360 STEP 15 ' hourly results: 15 degrees per hour
' solve a spherical triangle connecting the points
' north pole
' New York's position when satellite was overhead
' position directly underneath satellite now
' side from pole to NYC's old position is complement of NYC's latitude
' (90 - 41 = 49)
' side from old NYC position to current satellite position is
' 15 degrees per elapsed hour
' use law of cosines to find third side, but use arcsin instead of arccos
' as we seek the latitude of the satellite, which is the complement of
' the remaining side of the triangle:
lat = fnas(sin41 * fncs(i))
' (cosine of included angle, opposite side being found is zero as this
' is a right triangle-- meridian of NYC is perpendicular to orbit at its
' maximum latitude)
' solve for the angle at the north pole, which is the longitudinal motion:
lon = fnac((fncs(i) - sin41 * fnsn(lat)) / (cos41 * fncs(lat)))
IF i > 180 THEN lon = 360 - lon 'correct quadrant
' NYC has moved also, so subtract out NYC's rotational motion, which is
' a linear change in longitudinal position as the earth rotates:
lon = lon - i
'lon is relative motion eastward from NYC, so last item printed is
' the 74 degree *west* longitude of NYC *minus* the rel longitude
' of the satellite.
PRINT #1, USING "##### ###.# ###.# ###.#"; i; lat; lon; 74 - lon
NEXT
CLOSE 1
The resulting output is:
0 41.0 0.0 74.0
15 39.3 4.5 69.5
30 34.6 7.4 66.6
45 27.6 8.0 66.0
60 19.1 6.5 67.5
75 9.8 3.6 70.4
90 0.0 0.0 74.0
105 -9.8 -3.6 77.6
120 -19.1 -6.5 80.5
135 -27.6 -8.0 82.0
150 -34.6 -7.4 81.4
165 -39.3 -4.5 78.5
180 -41.0 -0.0 74.0
195 -39.3 4.5 69.5
210 -34.6 7.4 66.6
225 -27.6 8.0 66.0
240 -19.1 6.5 67.5
255 -9.8 3.6 70.4
270 -0.0 0.0 74.0
285 9.8 -3.6 77.6
300 19.1 -6.5 80.5
315 27.6 -8.0 82.0
330 34.6 -7.4 81.4
345 39.3 -4.5 78.5
360 41.0 -0.0 74.0
The left column represents degrees along the satellite's orbit, second is the latitude the satellite is above, third is the longitude difference, east of NYC; finally the actual west longitude of the satellite as NYC is 74 degrees west.
When you plot it out it looks like a figure-8.
Back to Marilyn Misses the Satellite