Plotting celestial objects with Hy and matplotlib

The theme of using Jupyter and Calysto Hy continue with this posting. This time I’m going to throw matplotlib into mix and see if I can make some neat looking graphs about planets and their trajectories.

I translated code by Paul Schlyter at http://www.stjarnhimlen.se/comp/ppcomp.html into Hy, in order to incorporate it into my little roguelike game. Games like Nethack have been aware of moon phases for the long time and modify the gameplay slightly based on that. I have similar idea, but I’m planning to also included planets and zodiac in the mix. This blog posting isn’t about the game though, but about plotting planet trajectories with matplotlib.

The code I’m using for Pyherc is written with Greek names for planets (I have sort of Greek theme going on with the game), but for this posting I’m using regular ones, thus we start by defining them:

(def sun 'helios)
(def moon 'selene)
(def mercury 'hermes)
(def venus 'aphrodite)
(def earth 'gaia)
(def mars 'ares)
(def jupiter 'dias)
(def saturn 'cronus)

Basic API is very simple and offers few functions to query heliocentric and geocentric position of an object or right ascension – declination pair for an object. Couple helper functions can be used to query for apparent angle between two objects, moon phase or in which zodiac some object is. All calculations are done for specific moment in time. But before we get started, lets import some libraries that we’ll need along the way.

(import [matplotlib :as mpl]
        [matplotlib [cbook]]
        [mpl_toolkits.mplot3d [Axes3D]]
        [matplotlib.pyplot :as plt]
        [datetime [datetime timedelta]])

So, lets get started calculating where Mars was on 1st of January and use several different coordinate systems.

(heliocentric-position mars (datetime 2016 1 1))
(-1.649022345024184, 0.16357386411721142, 0.044019515787414365)

Heliocentric position is calculated in coordinate system where Sun is at the origo (0, 0, 0). Distances are given in astronomical units (mean distance between Earth and Sun).

(geocentric-position mars (datetime 2016 1 1))
[-1.4784421312546756, -0.8048408970730346, 0.044019515787414365]

In geocentric system, Earth is in the origo and distances are again given in astronomical units.

(ra-decl-of mars (datetime 2016 1 1))
(207.0812659237103, -9.562456896419407)

Right ascension-declination pair is given in equatorial coordinate system (https://en.wikipedia.org/wiki/Equatorial_coordinate_system has better explanation than I could ever write). In this system Earth is again in origo and location is defined with two angles.

(house-of mars (datetime 2016 1 1))
'zygos'

This gives current zodiac house the object is in. It’s the most inaccurate way of describing a location of object in the library, but rather essential for my game programming purposes. Zygos is Libra in Greek (as far as I understand, I don’t actually speak Greek).

What does orbit of Mars look like then? Lets bust out matplotlib and graph it. First we’ll define some colours so in every graph each planet is consistently coloured with the same colour. All calculations are done for year 2016 and onward, thus we define basedate to be 1st of January, 2016. And in order to shift our view from API calls and manipulation of lists, we define couple helper functions.

(def colour-mercury "y")
(def colour-venus "g")
(def colour-earth "b")
(def colour-mars "r")
(def colour-jupiter "c")
(def colour-saturn "m")
(def colour-sun "y*")

(def basedate (datetime 2016 1 1))

(defn unpack-coordinates [coords]
    "split coordinate list into separate lists"
    (if (= (len (first coords)) 3)
      (, (list (map first coords))
         (list (map second coords))
         (list (map (fn [x] (get x 2)) coords)))
      (, (list (map first coords))
         (list (map second coords)))))

(defn heliocentric-orbit [body daycount]
  "calculate heliocentric orbit for given body and time"
  (unpack-coordinates (list (map (fn [x]
                                   (heliocentric-position body
                                                          (+ basedate (timedelta :days x))))
                                   (range daycount)))))

With these helpers at our disposal, plotting a heliocentric orbit of an planet should be easy. View is from above of ecliptic plane, looking towards the Sun.

(setv (, x y z) (heliocentric-orbit mars 365))

(.axis plt "equal")
(.legend plt :loc "upper right")

(.plot plt x y colour-mars :label "Mars")
(.plot plt 0 colour-sun)

(.show plt)

astro-1

That doesn’t look like a full orbit. Quick check reveals that Mars has orbit length of 687 days. Lets try again (otherwise same code, but this time with almost twice as long time).

(setv (, x y z) (heliocentric-orbit mars 687))

(.plot plt x y colour-mars :label "Mars")
(.plot plt 0 colour-sun)

(.axis plt "equal")
(.legend plt :loc "upper right")

(.show plt)

astro-2

Much better, that’s a full orbit around the Sun (marked by little yellow star roughly in the middle). This reveals how close to circle orbits actually are. For the longest time I was under impression that they would be much more elliptical, but that was just because most of the drawings I had seen exaggerated the shape to show it more clearly.

Lets try plotting all planets the system supports on the same graph. Here it’s actually starting to show that orbits aren’t completely round and centered (also, Mercury has slightly thicker line than Saturn, because of amount of orbits it completes while Saturn completes just one).

(def date-range 10760)

(setv (, x0 y0 z0) (heliocentric-orbit mercury date-range))
(setv (, x1 y1 z1) (heliocentric-orbit venus date-range))
(setv (, x2 y2 z2) (heliocentric-orbit earth date-range))
(setv (, x3 y3 z3) (heliocentric-orbit mars date-range))
(setv (, x4 y4 z4) (heliocentric-orbit jupiter date-range))
(setv (, x5 y5 z5) (heliocentric-orbit saturn date-range))

(.axis plt "equal")
(.axis plt [-11 11 -11 11])

(.plot plt x0 y0 colour-mercury :label "Mercury")
(.plot plt x1 y1 colour-venus :label "Venus")
(.plot plt x2 y2 colour-earth :label "Earth")
(.plot plt x3 y3 colour-mars :label "Mars")
(.plot plt x4 y4 colour-jupiter :label "Jupiter")
(.plot plt x5 y5 colour-saturn :label "Saturn")
(.plot plt 0 colour-sun)

(.legend plt :loc "upper right")

(.show plt)

astro-3

This graph is starting to show what kind of distances we’re talking about. Four innermost planets (Mercury, Venus, Earth and Mars) are so bunched up that it’s hard to tell which one is which on a quick glance. Then Jupiter and Saturn are suddenly considerable further away (and I’m not even plotting Uranus or Neptune here).

How would the same graph look like if centered on Earth (marked by blue circle)? For calculating geocentric orbit, we’ll define yet another helper function.

(defn geocentric-orbit [body daycount]
  "calculate geocentric orbit for given body and time"
  (unpack-coordinates (list (map (fn [x]
                                   (geocentric-position body
                                                        (+ basedate (timedelta :days x))))
                                 (range daycount)))))

We’re placing little blue dot at origo to represent Earth.

(def date-range 10760)

(setv (, x0 y0 z0) (geocentric-orbit mercury date-range))
(setv (, x1 y1 z1) (geocentric-orbit venus date-range))
(setv (, x2 y2 z2) (geocentric-orbit mars date-range))
(setv (, x3 y3 z3) (geocentric-orbit jupiter date-range))
(setv (, x4 y4 z4) (geocentric-orbit saturn date-range))

(.axis plt "equal")
(.axis plt [-11 11 -11 11])

(.plot plt x0 y0 colour-mercury :label "Mercury")
(.plot plt x1 y1 colour-venus :label "Venus")
(.plot plt x2 y2 colour-mars :label "Mars")
(.plot plt x3 y3 colour-jupiter :label "Jupiter")
(.plot plt x4 y4 colour-saturn :label "Saturn")
(.plot plt 0 "ob")

(.legend plt :loc "upper right")

(.show plt)

astro-4

That’s actually quite pretty. And quite complicated too. No wonder ancient astronomers had hard time reconciling their observations with geocentric model. You would have to build pretty elaborate system to model that. But the inner planets are kind of hard to see, so lets plot only them next (still using the values calculated in the previous step).

(setv (, x0 y0 z0) (geocentric-orbit mercury date-range))
(setv (, x1 y1 z1) (geocentric-orbit venus date-range))

(.axis plt "equal")

(.plot plt x0 y0 colour-mercury :label "Mercury")
(.plot plt x1 y1 colour-venus :label "Venus")
(.plot plt 0 "ob")

(.legend plt :loc "upper right")

(.show plt)

astro-5

That makes again quite beautiful looking pattern. You can see how sometimes Venus is further away than Mercury and sometimes it’s other way around. Venus also makes that five-looped spirograph pattern, while Mercury is covering more ground (or space rather).

Lets switch our location from far away in space and start looking how things look from our planet. I’m not going to take account things like exact position (longitude, latitude and height) of observer. Instead of that, we’re going to map planets against the stars using coordinates defined by right ascension and declination. In order to get our plot to look like planet’s trajectory across the sky, we’ll have to invert x-axis (so that 360 degrees is at left and 0 degrees is at right).

Another important thing to remember is that while planet’s trajectory is continuous, our graph won’t be. When right ascension gets larger that 360 degrees, it resets back to zero and vice versa. If we were to simply connect all calculated positions with straight lines, we might end up with nasty lines going all the way from one side to another. To get around this limitation, we’ll use scatter plot and just draw points at calculated locations. When amount of points is large enough, we get close enough representation of a line. And for that we need to define our colours once again, this time with modification that causes scatter plot to be drawn.

(defn ra-decl-path [body daycount]
  "calculate trajectory for given body and time"
  (unpack-coordinates (list (map (fn [x] (ra-decl-of body (+ basedate (timedelta :days x))))
                                 (range daycount)))))

(def scatter-mercury "y.")
(def scatter-venus "g.")
(def scatter-mars "r.")
(def scatter-jupiter "c.")
(def scatter-saturn "m.")

Actual plotting is just few lines of code:

(setv (, ra decl) (ra-decl-path mars 687))

(.axis plt [360 0 -90 90])

(.plot plt ra decl scatter-mars :label "Mars" :markersize 2)

(.legend plt :loc "upper right")

(.show plt)

astro-6

Because we’re plotting apparent movement of Mars, there’s funny little hook in the plot. This is called apparent retrograde motion and puzzled ancient astronomers quite a bit and quite elaborate models were constructed to explain this. But as soon as one switches from geocentric model to heliocentric model, apparent retrograde motion is not mystery anymore.

We can plot movement of all planets during a one year period and notice how they seem to be moving on a tilted path across the sky. This tilt is by-product of the coordinate system we’re using. Earth is tilted 23.5 degrees in relation to ecliptic plane (the plane the Sun appears to be orbiting Earth, or vice-versa, depending on how you look at things). Our coordinate system isn’t aligned with the ecliptic plane, but with the Earth’s equator, thus creating the tilt.

This same tilt is what causes changes in seasons. When Sun is at positive declination, northern hemisphere of Earth is getting more sunlight and there’s summer there. When declination of Sun as at the highest point (23.5 degrees) summer solstice occurs, day is the longest and night shortest. If you’re north enough (beyond polar circle, which is 23.5 from northpole), Sun doesn’t set at all. Eventually Sun will move towards zero declination and autumn solstice occurs. Day and night are equally long. As Sun reaches -23.5 degrees of declination, winter solstice occurs and night is the longest is northern hemisphere. Again, if you’re north side of polar circle, Sun doesn’t rise at all. Eventually Sun will move back to 0 degrees of declination, autumn solstice occurs and day and night are equally long again. Same things happens on the southern hemisphere, but in opposite order.

(def date-range 365)

(setv (, ra0 decl0) (ra-decl-path mercury date-range))
(setv (, ra1 decl1) (ra-decl-path venus date-range))
(setv (, ra2 decl2) (ra-decl-path mars date-range))
(setv (, ra3 decl3) (ra-decl-path saturn date-range))
(setv (, ra4 decl4) (ra-decl-path jupiter date-range))

(.axis plt [360 0 -90 90])

(.plot plt ra0 decl0 scatter-mercury :label "Mercury" :markersize 2)
(.plot plt ra1 decl1 scatter-venus :label "Venus" :markersize 2)
(.plot plt ra2 decl2 scatter-mars :label "Mars" :markersize 2)
(.plot plt ra3 decl3 scatter-jupiter :label "Jupiter" :markersize 2)
(.plot plt ra4 decl4 scatter-saturn :label "Saturn" :markersize 2)

(.legend plt :loc "upper right")

(.show plt)

astro-7

We can even use suitable star chart (I chose to use one from https://commons.wikimedia.org/wiki/File:CaldwellStarChart.svg) and plot trajectories of planets on top of that.

(setv image_file (.get-sample-data cbook "./1000px-CaldwellStarChart.svg.png"))
(setv image (.imread plt image-file))
(.figure plt :figsize (, 20 10))
(.imshow plt image :extent [360 0 -90 90])

(.axis plt [360 0 -90 90])

(.plot plt ra0 decl0 scatter-mercury :label "Mercury" :markersize 4)
(.plot plt ra1 decl1 scatter-venus :label "Venus" :markersize 4)
(.plot plt ra2 decl2 scatter-mars :label "Mars" :markersize 4)
(.plot plt ra3 decl3 scatter-jupiter :label "Jupiter" :markersize 4)
(.plot plt ra4 decl4 scatter-saturn :label "Saturn" :markersize 4)

(.legend plt :loc "upper right")

(.show plt)

astro-8

This made me think, how would one go about showing current location of planets? This is achieved with almost identical code. Instead of plotting multiple points in scatter plot, we just use one per planet.

(setv image-file (.get-sample-data cbook "./1000px-CaldwellStarChart.svg.png"))
(setv image (.imread plt image-file))
(.figure plt :figsize (, 20 10))
(.imshow plt image :extent [360 0 -90 90])

(.axis plt [360 0 -90 90])

(setv (, ra0 decl0) (ra-decl-of mercury (datetime.today)))
(setv (, ra1 decl1) (ra-decl-of venus (datetime.today)))
(setv (, ra2 decl2) (ra-decl-of mars (datetime.today)))
(setv (, ra3 decl3) (ra-decl-of jupiter (datetime.today)))
(setv (, ra4 decl4) (ra-decl-of saturn (datetime.today)))
(setv (, ra5 decl5) (ra-decl-of sun (datetime.today)))

(.plot plt ra0 decl0 scatter-mercury :label "Mercury" :markersize 10)
(.plot plt ra1 decl1 scatter-venus :label "Venus" :markersize 10)
(.plot plt ra2 decl2 scatter-mars :label "Mars" :markersize 10)
(.plot plt ra3 decl3 scatter-jupiter :label "Jupiter" :markersize 10)
(.plot plt ra4 decl4 scatter-saturn :label "Saturn" :markersize 10)
(.plot plt ra5 decl5 colour-sun :label "Sun" :markersize 10)

(.text plt ra0 decl0 "Mercury" :size "large")
(.text plt ra1 decl1 "Venus" :size "large")
(.text plt ra2 decl2 "Mars" :size "large")
(.text plt ra3 decl3 "Jupiter" :size "large")
(.text plt ra4 decl4 "Saturn" :size "large")
(.text plt ra5 decl5 "Sun" :size "large")

(.legend plt :loc "upper right")

(.show plt)

astro-9

As a final flourish, lets graph some 3d trajectories through space. Here we have six innermost planets of the solar system. In order to show the tilt of their orbits more clearly, I’m using different scale of y-axis. In reality the effect isn’t as pronounced as here.

(def date-range 10760)

(setv fig (.figure plt))
(setv ax (.gca fig :projection "3d"))
(.set-xlim3d ax -15 15)
(.set-ylim3d ax -15 15)
(.set-zlim3d ax -3 3)

(defn plot-projections [x y z c]
  "plot projections to three sides"
  (.plot ax x z c :zdir "y" :zs 15)
  (.plot ax y z c :zdir "x" :zs -15)
  (.plot ax x y c :zdir "z" :zs -3))

(setv (, x0 y0 z0) (heliocentric-orbit mercury date-range))
(setv (, x1 y1 z1) (heliocentric-orbit venus date-range))
(setv (, x2 y2 z2) (heliocentric-orbit earth date-range))
(setv (, x3 y3 z3) (heliocentric-orbit mars date-range))
(setv (, x4 y4 z4) (heliocentric-orbit jupiter date-range))
(setv (, x5 y5 z5) (heliocentric-orbit saturn date-range))

(.plot ax x0 y0 z0 "y")
(.plot ax x1 y1 z1 "g")
(.plot ax x2 y2 z2 "b")
(.plot ax x3 y3 z3 "r")
(.plot ax x4 y4 z4 "c")
(.plot ax x5 y5 z5 "m")

(plot-projections x0 y0 z0 "y")
(plot-projections x1 y1 z1 "g")
(plot-projections x2 y2 z2 "b")
(plot-projections x3 y3 z3 "r")
(plot-projections x4 y4 z4 "c")
(plot-projections x5 y5 z5 "m")
(.plot ax [0] [0] [0] "y*")

(.show plt)

astro-10

One can also use same routines to calculate moon phases (just apparent angle between Sun and Moon) or predict eclipses (relative location of Sun, Moon and Earth), although there are simpler routines for those. It’s also possible to check when solstices and equinoxes occur, find next conjuction of given planets and so on. Possibilities are many and more can be discovered by simply playing with data and charting it.

Advertisements

Leave a Reply

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

WordPress.com Logo

You are commenting using your WordPress.com 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 )

Google+ photo

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

Connecting to %s