Clojure Astronomy

Mapping orbital positions

In the previous post I showed how to map an elliptical orbit around the sun, but I didn't show where any planets were in their orbits. Due to the inconstant speed of objects in elliptical orbits, finding their exact positions requires some numerical approximation.

Orbital parameters for a body will usually give three parameters that we can use to determine its location on any date: orbital period, mean anomaly at epoch, and the epoch itself. We'll use the common epoch of J2000, since that's what Wikipedia gives in its sidebars on the planets.

(require '[cljc.java-time.offset-date-time :as odt]
         '[cljc.java-time.temporal.chrono-unit :as cu]
         '[cljc.java-time.zone-offset :as zo])

(def J2000 (odt/of 2000 1 1 12 0 0 0 zo/utc))

The positional parameters for J2000:

Body	Period (days)	Mean anomaly at epoch
Mercury	87.97	174.796°
Venus	224.70	50.115°
Earth	365.26	358.617°
Mars	686.98	19.412°
Jupiter	4332.59	20.02°
Saturn	10759.22	317.02°
Uranus	30688.50	142.2386°
Neptune	60195.00	256.228°
Pluto	90560.00	14.53°

The first step in calculating an orbiting body's position is to convert the mean anomaly at the epoch to the mean anomaly on some other date. The mean anomaly indicates how far around the orbit the body would be on a circular orbit, so to find the mean anomaly for another date, we calculate how many orbits the body could have completed since the epoch and add it to the mean anomaly at the epoch:

(defn mean-anomaly [{:keys [mean-anomaly-at-epoch period epoch]} time]
  (let [days-difference (cu/between cu/days epoch time)]
    (mod (Math/toDegrees
           (+ (Math/toRadians mean-anomaly-at-epoch)
              (/ (* 2 Math/PI days-difference)
                 period)))
         360)))

Testing it out:

(def today (odt/of 2022 11 17 0 0 0 0 zo/utc))

(map (juxt :name
           #(mean-anomaly % today))
     bodies)

(["Mercury" 166.34735155412636]
 ["Venus" 115.90322524154362]
 ["Earth" 313.38114007621516]
 ["Mars" 77.70540009898468]
 ["Jupiter" 354.24677890130374]
 ["Saturn" 236.5755811666645]
 ["Uranus" 240.24925545725597]
 ["Neptune" 306.19560528283085]
 ["Pluto" 47.74333922261484])

Before we can get the actual angle of a body around the sun, we need to convert mean anomaly to eccentric anomaly, which is the angle of the body on an ellipse instead of a circle, but with the sun at the center instead of one focus. While the eccentric anomaly can be used to calculate the mean anomaly with

it's not so easy to solve for E, so we'll have to find it numerically. We can use Newton's method, providing the mean anomaly as an initial guess:

(defn newtons-method [f f' x0 tolerance]
  (loop [x x0]
    (let [x' (- x (/ (f x) (f' x)))]
      (if (< tolerance (Math/abs (- x x')))
        (recur x')
        x'))))

(defn eccentric-anomaly [mean-anomaly eccentricity]
  (let [mean-anomaly (Math/toRadians mean-anomaly)]
    (Math/toDegrees
      (newtons-method
        #(- % (* eccentricity (Math/sin %)) mean-anomaly)
        #(- 1 (* eccentricity (Math/cos %)))
        mean-anomaly
        1e-7))))

The values are only slightly different, even for more eccentric Mercury and Pluto:

(map (juxt :name
           (fn [{:keys [eccentricity] :as body}]
             (-> body
                 (mean-anomaly today)
                 (eccentric-anomaly eccentricity))))
     bodies)

(["Mercury" 168.6633301234182]
 ["Venus" 116.25121452249438]
 ["Earth" 312.6773257557783]
 ["Mars" 83.0171319577618]
 ["Jupiter" 353.9515591810506]
 ["Saturn" 233.9580169682439]
 ["Uranus" 237.95832738958006]
 ["Neptune" 305.7922948814424]
 ["Pluto" 60.10127519036517])

With the eccentric anomaly we can finally calculate the true anomaly, which is the actual angle of a body in an elliptical orbit with the sun at one focus. The calculation is

I've chosen this form not just for the symmetry of the expression but also because this fractional form allows us to use Java's Math.atan2 function to avoid quadrant calculations:

(defn true-anomaly [eccentric-anomaly eccentricity]
  (let [eccentric-anomaly (Math/toRadians eccentric-anomaly)]
    (Math/toDegrees
      (* 2 (Math/atan2 (* (Math/sqrt (+ 1 eccentricity))
                          (Math/sin (/ eccentric-anomaly 2)))
                       (* (Math/sqrt (- 1 eccentricity))
                          (Math/cos (/ eccentric-anomaly 2))))))))

(map (juxt :name
           (fn [{:keys [eccentricity] :as body}]
             (-> body
                 (mean-anomaly today)
                 (eccentric-anomaly eccentricity)
                 (true-anomaly eccentricity))))
     bodies)

(["Mercury" 170.78759627285143]
 ["Venus" 116.59868632821262]
 ["Earth" 311.96946226397176]
 ["Mars" 88.36707110285339]
 ["Jupiter" 353.6487977391962]
 ["Saturn" 231.38165947045587]
 ["Uranus" 235.69476079075423]
 ["Neptune" 305.38795245840817]
 ["Pluto" 73.44020657430707])

Let's finally draw bodies in their correct positions. Here are a couple functions similar to ellipse and orbit from the previous post:

(defn elliptical-position [{:keys [semi-major-axis eccentricity] :as body} time]
  (let [true-anom (-> body
                      (mean-anomaly time)
                      (eccentric-anomaly eccentricity)
                      (true-anomaly eccentricity))
        dist (distance semi-major-axis eccentricity)]
    (polar->cartesian (dist true-anom) true-anom)))

(defn orbit-position
  [{:keys [argument-of-perihelion
           inclination
           longitude-of-ascending-node]
    :as body}
   time]
  (->> (elliptical-position body time)
       (transform (rotate-z argument-of-perihelion))
       (transform (rotate-x inclination))
       (transform (rotate-z longitude-of-ascending-node))))

(require '[blog.math :refer [linear]]
         '[blog.svg :refer [view-box translate points->d]])

To check our work, compare the above view with WolframAlpha or a more sophisticated map for November 17, 2022.

November 2022