Hazard

Drawing the Solar System January 2, 2017

Several years ago I wanted to calculate the distance between Earth and Mars given an arbitrary future date. I wasn’t satisfied by the free solar system simulators I found, so I looked into calculating it myself. The write-ups I found usually offered long, non-intuitive expressions of four or five trig functions. I eventually discovered PyEphem would do what I wanted, but I never quite let go of being able to do the calculations myself; there had to be a more intuitive way to express the math. Turns out, there is.

A planet’s elliptical orbit around the sun can be described by a handful of parameters. Wikipedia is pretty good about listing them in a sidebar. Here are Earth’s orbital parameters:

We need eight values: the epoch, semi-major axis, eccentricity, orbital period, mean anomaly, inclination, longitude of the ascending node, and argument of perihelion. The semi-major axis and eccentricity describe an ellipse. The inclination, longitude of the ascending node, and argument of perihelion describe how the ellipse is oriented in three dimensions. The mean anomaly indicates where the planet is at the epoch; it’s the only time-dependent parameter. Finally, to calculate where the planet is at any given moment, we need the orbital period.

The epoch J2000 is noon on January 1, 2000. Here’s a Clojure map of Earth’s orbital parameters.


(def J2000 (js/Date. 2000 0 1 12 0 0))

(def Earth
  {:epoch J2000
   :semi-major-axis 149.598023e9 ; meters
   :eccentricity 0.0167086
   :period 365.256363004 ; Earth days
   :mean-anomaly-at-epoch 358.617 ; degrees
   :inclination 0.00005 ; degrees
   :longitude-of-ascending-node -11.26064 ; degrees
   :argument-of-perihelion 114.20783 ; degrees
   })

Converting these values into a heliocentric position has two major steps: turning the mean anomaly into the true anomaly, and orienting the ellipse in three-dimensional space.

In orbital mechanics, the term “anomaly” refers to an angle between a planet’s current position, the sun, and where the planet is when it’s closest to the sun. There are three of these anomalies: mean anomaly, eccentric anomaly, and true anomaly. We have mean anomaly, but we need true anomaly to calculate the actual heliocentric position of a planet. We’ll convert from mean anomaly to true anomaly using as a go-between the eccentric anomaly.

Orbital parameters typically give a planet’s mean anomaly at the epoch. The mean anomaly is the orbital angle of a planet in a hypothetical, perfectly circular orbit. In such an orbit, the planet would move at constant speed, so it’s simple to convert the mean anomaly at the epoch to the mean anomaly at some other date knowing nothing more than the planet’s orbital period. Here’s the calculation:


(defn mean-anomaly [mean-anomaly-at-epoch period epoch time]
  (+ mean-anomaly-at-epoch
     (/ (* 2 js/Math.PI (- time epoch))
        period)))

Unfortunately, planets aren’t actually in perfectly circular orbits, so the mean anomaly can’t be used directly for calculating a planet’s position. Once we have the current mean anomaly, we can calculate the eccentric anomaly.

The eccentric anomaly is the orbital angle of a planet supposing it’s in the proper elliptical orbit, but with the sun at the center of the ellipse rather than at one focus. It’s a useful intermediate step between the mean anomaly and the true anomaly.

Here’s the hardest part of the whole computation. While it’s easy to convert eccentric anomaly to mean anomaly, there’s no straight-forward way to go the other direction. We have to compute it numerically. Below is a conversion function that uses Newton’s method, with 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 (js/Math.abs (- x x')))
        (recur x')
        x'))))

(defn eccentric-anomaly [eccentricity mean-anomaly]
  (let [f #(- % (* eccentricity (js/Math.sin %)) mean-anomaly)
        f' #(- 1 (* eccentricity (js/Math.cos %)))]
    (newtons-method f f' mean-anomaly 1e-7)))

With the eccentric anomaly, we can calculate the true anomaly. The true anomaly is the planet’s orbital angle on the real ellipse, with the sun at one focus. The following function uses JavaScript’s Math.atan2 function, which takes two arguments rather than one fraction like Math.atan. Math.atan2 won’t be confused by a negative numerator and a negative denominator creating a positive fraction, so it will always give us an angle in the correct quadrant.


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

A planet’s distance from the sun depends on its orbit’s semi-major axis, eccentricity, and true anomaly:


(defn radius [semi-major-axis eccentricity true-anomaly]
  (/ (* semi-major-axis (- 1 (js/Math.pow eccentricity 2)))
     (+ 1 (* eccentricity (js/Math.cos true-anomaly)))))

We can now express a planet’s position in polar coordinates, using the radius and true anomaly. Here’s a function to convert those polar coordinates to three-dimensional Cartesian coordinates:


(defn polar->cartesian [radius angle]
  [(* radius (js/Math.cos angle))
   (* radius (js/Math.sin angle))
   0])

Notice that all the coordinates given by the above function will be in a two-dimensional plane. We still need to rotate the ellipse into position. Here’s Wikipedia’s diagram showing us how to use inclination, longitude of the ascending node, and the argument of perihelion, the three orbital parameters we haven’t yet factored in:

And here’s how the typical write-up suggests you do this sort of transformation:

These are the lengthy, non-intuitive expressions that never helped me understand what was really going on. Fortunately, there’s a better way to express this, and it helps to break it down step-by-step. We’ve already pulled two steps out of these equations: finding the polar coordinates and converting to Cartesian coordinates. But there are three more ideas still braided together that we can tease apart.

Let’s think about what we’re going to do to our ellipse to get the three-dimensional orientation Wikipedia shows. First, we’re going to rotate the orbit around the z-axis by the argument of perihelion. Then we’re going to rotate the orbit around the x-axis by the inclination. Finally, we’re going to rotate the orbit around the z-axis again, this time by the longitude of the ascending node. I have some practice visualizing these transformations, but I would love it if someone created a step-by-step animation.

Rotations like this can be expressed as transformation matrices, and that’s the key to making the math more transparent. I highly recommend Grant Sanderson’s Essence of Linear Algebra series for better understanding matrices as transformations.

We need two kinds of matrices, one that transforms a point by rotating it around the z-axis, and one that rotates a point around the x-axis. Here are two functions that each produce a rotation matrix:


(def sin js/Math.sin)
(def cos js/Math.cos)

(defn rotate-x [angle]
  [[1 0           0]
   [0 (cos angle) (- (sin angle))]
   [0 (sin angle) (cos angle)]])

(defn rotate-z [angle]
  [[(cos angle) (- (sin angle)) 0]
   [(sin angle) (cos angle)     0]
   [0           0               1]])

Here’s an implementation of matrix-vector multiplication that simply takes the dot product between the coordinate and each row of a matrix:


(defn dot-product [vector1 vector2]
  (reduce + (map * vector1 vector2)))

(defn transform [matrix coord]
  (map #(dot-product % coord) matrix))

Now we can put it all together. Here’s a function that takes a map of orbital parameters and a time, converts the parameters to the right units, carries out the operations we’ve defined, and outputs the position of a planet:


(def days->seconds #(* % 24 60 60))
(def date->seconds #(/ (.getTime %) 1000))
(def degrees->radians #(* % js/Math.PI (/ 180)))

(defn position [params time]
  (let [{:keys [epoch semi-major-axis eccentricity period mean-anomaly-at-epoch
                inclination longitude-of-ascending-node argument-of-perihelion]} params
        mean-anom (mean-anomaly (degrees->radians mean-anomaly-at-epoch)
                                (days->seconds period)
                                (date->seconds epoch)
                                (date->seconds time))
        eccentric-anom (eccentric-anomaly eccentricity mean-anom)
        true-anom (true-anomaly eccentricity eccentric-anom)
        radius (radius semi-major-axis eccentricity true-anom)]
    (->> (polar->cartesian radius true-anom)
      (transform (rotate-z (degrees->radians argument-of-perihelion)))
      (transform (rotate-x (degrees->radians inclination)))
      (transform (rotate-z (degrees->radians longitude-of-ascending-node))))))

Let’s try it. Where was Earth at the New Year?


(def date (js/Date. 2017 0 1)) ; January 1, 2017

(position Earth date)

Not very illuminating. Let’s draw it’s position relative to the sun:


(draw [{:color "yellow" :position [0 0 0]}
       {:color "dodgerblue" :position (position Earth date)}])

Let’s throw in a few more planets:


(def Mercury
  {:epoch J2000
   :semi-major-axis 57.90905e9
   :eccentricity 0.205630
   :period 87.969
   :mean-anomaly-at-epoch 174.796
   :inclination 7.005
   :longitude-of-ascending-node 48.331
   :argument-of-perihelion 29.124})

(def Venus
  {:epoch J2000
   :semi-major-axis 108.208e9
   :eccentricity 0.006772
   :period 224.701
   :mean-anomaly-at-epoch 50.115
   :inclination 3.39458
   :longitude-of-ascending-node 76.860
   :argument-of-perihelion 54.884})

(def Mars
  {:epoch J2000
   :semi-major-axis 227.9392e9
   :eccentricity 0.0934
   :period 686.971
   :mean-anomaly-at-epoch 19.373
   :inclination 1.850
   :longitude-of-ascending-node 49.558
   :argument-of-perihelion 286.502})

(draw [{:color "yellow" :position [0 0 0]}
       {:color "gray" :position (position Mercury date)}
       {:color "orange" :position (position Venus date)}
       {:color "dodgerblue" :position (position Earth date)}
       {:color "crimson" :position (position Mars date)}])

Is it right? Here’s what Wolfram Alpha shows:

Pretty good! Let’s check one more thing: the distance between Earth and Mars. Here’s what Wolfram Alpha says:


(defn dist [p1 p2]
  (js/Math.sqrt
    (reduce +
      (map (fn [a b] (js/Math.pow (- a b) 2))
           p1 p2))))

(dist (position Earth date) (position Mars date))

Not quite the same. How far off is it?


(let [real 2.462e11
      calculated (dist (position Earth date) (position Mars date))]
  (/ (js/Math.abs (- real calculated)) real))

Less than a percent off. The difference could be a quirk of the exact orbital parameters Wikipedia gives and what Wolfram Alpha uses. Since I’m not sending any actual spaceships to Mars, that error is low enough for me.

Here’s a more complete solar system mapper I’ve created, one that shows orbital paths, includes all the planets and the biggest asteroids, supports panning and zooming, and lets you see the distance between any two bodies at any given date. The source code is on GitHub.

Happy travels!