Celestial Calculations: A Gentle Introduction to Computational Astronomy

136 Chapter 6 Note that equation 6.2.7 calculates the number of Julian centuries the standard epoch is from 12h UT (noon) on January 0, 1900. Using the equation of the center to find the true anomaly, let’s now work through a complete example to find the Sun’s location on February 5, 2015, at 12h LCT. Assume that an observer is in the Eastern Standard Time zone at 78◦ W longitude, 38◦ N latitude and is not on daylight saving time. 1. Use the equations from chapter 3 to perform all needed time conversions. That is, compute UT, GST, and LST from the given LCT. Adjust the date if necessary. (Ans: UT = 17.0h, GST = 2.035013h, LST = 20.835013h, Date = 2/5/2015.) 2. Compute the Julian day number JDe for the standard epoch. Be sure to include the fractional part of the day. (Ans: JDe = 2,451,545.0 for J2000.) 3. Compute the Julian day number JD for the desired date. Use the Greenwich date and UT from step 1, not the LCT time and date, and be sure to include the fractional part of the day. From step 1, we need the Julian day number for 17.0h UT on 2/5/2015. (Ans: JD = 2,457,059.20833.) 4. Compute De, the total number of elapsed days since the standard epoch, by subtracting JDe from JD. (Ans: De = 5514.20833 days.) 5. Use equation 6.2.3 to compute M . to the range of 0◦ to 360◦. (Ans: M = 5432.592589.◦) 6. Add or subtract multiples of 360◦ to adjust M (Hint: use M mod 360◦.) (Ans: M = 32.592589◦.) 7. Use equation 6.2.4 to approximate the equation of the center. (Ans: Ec = 1.031320◦.) 8. Add Ec to M to get the true anomaly. to the range of 0◦ to 360◦. (Ans: υ = 33.623909◦.) 9. Add or subtract multiples of 360◦ to adjust υ (Ans: no adjustment is necessary.) 10. Add υ and g to get the ecliptic longitude. (Ans: λ = 316.562255◦.) 11. If λ > 360◦, subtract 360◦. (Ans: no adjustment is necessary.)

The Sun 137 At this point, the Sun’s ecliptic coordinates are λ ecliptic longitude, 0◦ ecliptic latitude. 12. Convert the ecliptic coordinates from the prior step to equatorial coordi- nates. (Ans: α = 21.267801h, δ = −15.872529◦.) 13. Finally, convert the equatorial coordinates to horizon coordinates. (Remember that converting equatorial to horizon coordinates requires that the right ascension α from the previous step be converted to an hour angle.) (Ans: h = 35◦47 01 , A = 172◦17 46 .) As we pointed out earlier, we can solve either the equation of the cen- ter or Kepler’s equation to obtain the true anomaly from the mean anomaly. Let’s repeat the previous example to determine the Sun’s ecliptic coordinates by using the simple iterative scheme from subsection 4.5.5 to solve Kepler’s equation. Steps 1–6 of the algorithm above are identical whether solving the equation of the center or solving Kepler’s equation. Thus, our objective is to solve Kepler’s equation to determine the true anomaly from the mean anomaly found in step 6 (M = 32.592589◦). Recall from subsection 4.5.5 that we will iteratively compute the equation Ei = Mr + e sin Ei−1, where both the eccentric anomaly and mean anomaly are expressed in radians, not degrees. Thus, we must first convert M from degrees to radians, giving us Mr = Mπ = 32.592589◦π ≈ 0.568848 radians. 180◦ 180◦ This gives us a first estimate of the eccentric anomaly as E0 = Mr = 0.568848 radians. The next few iterations follow, where is the difference between successive approximations to the eccentric anomaly. E1 ≈ 0.577848, = 0.009000, E2 ≈ 0.577974, = 0.000126, E3 ≈ 0.577976, = 0.000002. Only 3 iterations are necessary to obtain a sufficiently accurate approxima- tion for the eccentric anomaly (0.577976 radians, which is ≈ 33.115588◦). The rapid convergence to a solution is achieved because the Earth-Sun orbital eccentricity is so small.

138 Chapter 6 Given the eccentric anomaly, we can obtain the true anomaly by applying equation 4.5.8, which is υ = 1+e E tan 1 − e tan 2 . 2 Doing so gives υ ≈ 33.642307◦. This differs from the result in step 8 by about 0.018◦, or by a little over 1 arcminute. Solving Kepler’s equation replaces steps 7 and 8, after which we merely continue with steps 9–13 to determine the Sun’s horizon coordinates. Using the true anomaly obtained by solving Kepler’s equation instead of the value from the equation of the center, we obtain λ = 316.580653◦, α = 21.269018h, δ = −15.867003◦, h = 35◦47 13 , A = 172◦16 25 . The method presented here for determining the position of the Sun, whether by solving the equation of the center or Kepler’s equation, is usually accurate to within about 1 , which is sufficient for the purposes of this book. For the best accuracy, however, one should use a standard epoch as close as possible to the date for which the Sun’s position is desired, or use equations 6.2.7 through 6.2.10 to adjust the Sun’s orbital elements to the desired date. More accurate algorithms can be found from other sources, such as The Astronomical Almanac or Meeus’s Astronomical Algorithms. Although more accurate algorithms exist, the method presented here is sufficient to demon- strate how orbital elements, the anomalies, and the various coordinate systems are interrelated. 6.3 Sunrise and Sunset As with locating a star in horizon coordinates, there is no guarantee that the Sun will be visible in the sky at the given date, time, and location. For exam- ple, let’s use the algorithm from the previous section (and solve the equation of the center) to compute the Sun’s location at 20h LCT on February 5, 2015, for an observer in the Eastern Standard Time zone at 78◦ W longitude, 38◦ N latitude. Assume the observer was not on daylight saving time. For that observer, date, and time, the Sun’s horizon coordinates were h = −28◦01 04 , A = 271◦26 04 , which was well below the observer’s horizon. This should not be a surprise because it is easy to know when the Sun will or will not