Calculating planetary motions

Celestial bodies and their trajectories across the skies have always fascinated people. Sun, moon and stars follow very predictable paths and one can predict certain events with great accuracy. Many religious events have been tied to timing of some special event (Easter for example, Sunday after first full moon after vernal equinox). So no wonder people have long kept track of these and tried to predict when things are about to happen.

Planets on the another hand seemed to be really tricky beasts. They travel independent of stars and can even reverse direction of their path on sky for a brief time. First geocentric models had all objects travel in circular paths around the Earth and failed to explain why this was the case. More complicated models were constructed, that could explain movement of some of the planets (but not all of them) with more complicated mechanisms. Only when it was realized that planets, including the Earth, were traveling around the Sun, these backwards (retrograde) movements started to make sense. However, ancients were able to do some pretty impressive feats of calculations with tools that they had in their disposal. One of the most famous ones is Antikythera mechanism.

But this blog post isn’t about ancient history though or even about specifically about calculating planetary positions, but what I’m planning to do with such a system (and little bit about the API).

How to compute planetary positions is a really good tutorial with examples that will walk through steps required to calculated position of planets and moon. I have translated the code in Hy and stuffed it in my game. At the time of writing this, I’m not using the code for anything in the game. The plan is to eventually use it to calculate in which zodiac each planet is and slightly adjust the game based on that. I don’t have very specific plans yet, but I’ll try to tie adjustments into mythology of zodiac and planets. Maybe when Ares (or Mars, but I’m using greek names for flavour) is in Tavros (Taurus), various large animals are stronger and more aggressive. Or when Hermes (Mercury) is in Zygos (Libray) luck is in favour of player and they will receive more bounty.

I wanted to keep the API as simple as possible. The main entry point for figuring out in which zodiac given planet (or sun or moon) is shown below:

(defn house-of [body d]
  "calculate zodiac house for given object"
  (ra-to-zodiac (normalize-angle (first (ra-decl-of body d)))))

So we calculate right ascension of the body, normalize the angle to be between 0 and 360 degrees and then convert this to zodiac. To calculate angle between two bodies (mainly used for figuring moon phases), I use following snippet:

(defn angle-between [body-1 body-2 d]
  "calculate angle between two bodies (RA)"
  (let [[(, ra1 decl1) (ra-decl-of body-1 d)]
        [(, ra2 decl2) (ra-decl-of body-2 d)]]
    (abs (- (normalize-angle ra1)
            (normalize-angle ra2)))))

We calculate right ascension for both bodies and calculate difference between normalized values. This value we don’t have to normalize again, as taking absolute value of the result is enough to ensure it’s between 0 and 360 degrees.

What about figuring out zodiacs? Since I had already decided that I only need one set of zodiacs, I could use global variable to hold necessary data:

(setv zodiacs [(, 0 30 'ihtheis)
               (, 30 60 'krios)
               (, 60 90 'tavros)
               (, 90 120 'didimoi)
               (, 120 150 'karkinos)
               (, 150 180 'leon)
               (, 180 210 'parthenos)
               (, 210 240 'zygos)
               (, 240 270 'skorpios)
               (, 270 300 'toksotis)
               (, 300 330 'aigokeros)
               (, 330 360 'ydrohoos)])

(defn ra-to-zodiac [ra]
  "convert RA to respective zodiac"
  (get (first (filter (fn [x]
                        (<= (first x) (normalize-angle ra) (second x)))

Similarly, orbital elements of planets are stored in a dictionary, where one can quickly locate them by given name of a planet. One tiny catch to know about this is that I’m storing them as functions, as orbital elements slowly change over time:

(setv orbits {})

(defn register-orbit [body func]
  (assoc orbits body func))

(defn new-orbit [name date long-ascending-node inclination arg-perihelion
                 semi-major-axis eccentricity mean-anomaly]
  "create data structure representing orbital elements"
  {:name name
   :date date
   :long-ascending-node (normalize-angle long-ascending-node)
   :inclination (normalize-angle inclination)
   :arg-perihelion (normalize-angle arg-perihelion)
   :semi-major-axis semi-major-axis
   :eccentricity eccentricity
   :mean-anomaly (normalize-angle mean-anomaly)})

(register-orbit 'hermes
                (fn [d]
                  "orbital elements of mercury/hermes on a given date"
                  (new-orbit :name 'hermes
                             :date d
                             :long-ascending-node (+ 48.3313 (* 3.24587E-5
                                                                (ast-date d)))
                             :arg-perihelion (+ 29.1241 (* 1.01444E-5
                                                           (ast-date d)))
                             :semi-major-axis 0.387098
                             :eccentricity (+ 0.205635 (* 5.59E-10
                                                          (ast-date d)))
                             :mean-anomaly (+ 168.6562 (* 4.0923344368
                                                          (ast-date d)))
                             :inclination (+ 7.0047 (* 5.00E-8
                                                       (ast-date d))))))

ast-date is a function used to calculate day number since 1st of January 2000. I added a little bit extra functionality in it, so that one can use both date and datetime objects with it. For objects that travel relatively slowly on the sky, time of the day doesn’t make that big difference. However, moon for example can travel roughly 12 degrees during a day, which is quite considerable amount.

(defn ast-date [d]
  "create date representation"
  (+ (* (. d year) 367)
     (- (// (* 7 (+ (. d year) (// (+ (. d month) 9) 12))) 4))
     (// (* (. d month) 275) 9)
     (. d day)
     (/ (if (hasattr d "hour")
          (. d hour)

I’m also normalizing all angles just in case. This again ensures that they’re between 0 and 360 degrees.

Speaking of angles. Trigonometric operations in Hy (and Python and most programming languages) usually use radians. But orbital elements are given in degrees and many calculations expected degrees. I could have been converting between degrees and radians when needed, but I knew that it’s easy to forget at somepoint and then wonder why the results are incorrect. Instead of that, I ended up writing little functions that did this calculation automatically for me. This way I could use degrees everywhere and computer could happily perform calculations using radians:

(defn sinᵒ [x]
  (math.sin (math.radians x)))

(defn cosᵒ [x]
  (math.cos (math.radians x)))

(defn atan2ᵒ [y x]
  (math.degrees (math.atan2 y x)))

Currently code can calculate positions of planets and the moon and find out which zodiac the position corresponds to. In addition it can calculate angle between two objects, which in turn can be used to calculate moon phases. It’s full moon when right ascension of moon and sun are 180 degrees apart. Likewise, when they’re same, it’s new moon.

Eventually I might add other functions that might come handy: checking if equinox or solstice falls on current date. All kinds of eclipses, conjuctions and such might be interesting too. Their effects will be quite tricky to test, so if/when I add those, I’ll also add method of setting game’s internal clock to specific date and time.

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s