Mapping elliptical orbits

Calculating the position of a body in orbit of the sun only requires a handful of parameters, so in this post and the next I'll walk through the steps. This post will show how to orient an ellipse around the sun, while the next post will work out where on that ellipse the body is on an arbitrary date.

Wikipedia articles on objects in the solar system the sun often list their orbital parameters, and while terms like "longitude of the ascending node" and "argument of perihelion" sound highly technical, their use in orienting an orbit around the sun is straight-forward. If you want to understand the terminology better, I suggest Wikipedia's article on orbital elements. Here I'll focus mainly on implementation.

Here are the parameters we need for calculating orbital ellipses:

The planets' orbits (aside from Mercury) have low eccentricity and all but circular orbits, so I've included Pluto in the diagrams below to make it easier to see that we're working with ellipses.

The distance of an ellipse from one of its foci is given by the equation

where a is the ellipse's semi-major axis and e is its eccentricity. In Clojure,

(defn distance [semi-major-axis eccentricity]
  (fn [angle]
    (/ (* semi-major-axis (- 1 (* eccentricity eccentricity)))
       (+ 1 (* eccentricity (Math/cos (Math/toRadians angle)))))))

The above orbital ellipses can now be calculated in two dimensions:

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

(defn ellipse [{:keys [semi-major-axis eccentricity]}]
  (let [dist (distance semi-major-axis eccentricity)]
    (map (fn [angle]
           (polar->cartesian (dist angle) angle))
         (range 0 360))))

This isn't an accurate diagram of the solar system. For starters, the planets' closest approaches to the sun are not all on the same side. To rotate perihelia to the correct positions, we'll use matrix multiplication:

(defn dot-product [v1 v2]
  (reduce + (map * v1 v2)))

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

Here's a matrix to rotate around the Z-axis:

(defn rotate-z [angle]
  (let [theta (Math/toRadians angle)]
    [[(Math/cos theta) (- (Math/sin theta)) 0]
     [(Math/sin theta) (Math/cos theta)     0]
     [0                0                    1]]))

Building on ellipse, we rotate each coordinate around the Z-axis:

(defn orbit
  [{:keys [argument-of-perihelion]
    :as body}]
  (sequence
    (map (partial transform (rotate-z argument-of-perihelion)))
    (ellipse body)))

This still isn't quite right, since there's a second rotation around the Z-axis we still need to do. First, however, we need to add inclination, which requires rotation around the X-axis:

(defn rotate-x [angle]
  (let [theta (Math/toRadians angle)]
    [[1 0                0]
     [0 (Math/cos theta) (- (Math/sin theta))]
     [0 (Math/sin theta) (Math/cos theta)]]))

(defn orbit
  [{:keys [argument-of-perihelion
           inclination]
    :as body}]
  (sequence
    (comp
      (map (partial transform (rotate-z argument-of-perihelion)))
      (map (partial transform (rotate-x inclination))))
    (ellipse body)))

It's hard to see the change since we're looking top-down, but notice the gap where Pluto's orbit dips inside Neptune's: between the last map and this one, it's become noticeably larger, an effect of Pluto's orbit tipping to an inclination of 17° while Neptune's is less than 2°.

The last step to get the orbits in the right places is to rotate them around the Z-axis by the longitude of the ascending node:

(defn orbit
  [{:keys [argument-of-perihelion
           inclination
           longitude-of-ascending-node]
    :as body}]
  (sequence
    (comp
      (map (partial transform (rotate-z argument-of-perihelion)))
      (map (partial transform (rotate-x inclination)))
      (map (partial transform (rotate-z longitude-of-ascending-node))))
    (ellipse body)))

There we have it. Pluto's orbit looks the same as on Wolfram Alpha.

Locating the planets on these orbits requires some basic numerical approximation that we'll explore in the next post.