Cataclysmic Mutation

Machine Learning and Whatever Else

Debunking the “Supermoon” With Clojure

So some jackass realized that the moon doesn’t orbit in a perfect circle and decided that he alone in the annals of history was in possession of the information that Earthquakes and other doom might ensue, and since the Moon was closer than usual (more on that in a bit), something *bad* was about to happen. No big deal. There are about 7,000,000,000 people on planet Earth, and a good percentage of them are, even now, right as you’re reading this, drooling uncontrollably from their slack jaws. However, this happened to make what now passes for news just before a massive earthquake hit off the coast of Japan, and so now the idiots who have taught themselves to operate a television are convinced that they’re correct. Oh, dear.

There are so many problems here, it’s hard to know where to start. Possibly that the moon isn’t actually especially close to Earth right now, and couldn’t cause this earthquake even if it was. Anyway, Phil Plait has a nice takedown up right now, and he knows a lot more about planetary physics than I do anyway, so just go read what he has to say. I’m writing this for another reason. One of my Facebook friends wondered what the actual plot of earthquake frequency and magnitude versus lunar distance actually would look like. That sounded like a fun little exercise in programming.

There are basically two questions that need to be answered. First, when have there been major earthquakes? I went to Wikipedia and grabbed a list of the 24 largest earthquakes in known history (today’s Japanese quake comes in number six). This gives me the magnitudes, but also the dates the quakes occurred. The range of dates covered by these 24 mega-quakes spans from a 1575 Chilean quake right up to today’s Japanese quake. The second thing we need is an ability to calculate the distance between the Earth and moon on any given day in history. That’s a tougher problem, but a bit of Google-fu turned up the ELP2000-82B model of lunar ephemerides. That’s what we need here. However, it’s a very complex model, and I didn’t want to spend a week debugging scientific code, but fortunately, I found an implementation on Github.

A bit more reading of code and comments showed that the model needs Julian dates. A bit more searching to find the exact formula for this conversion, and a quick implementation of said formula in Clojure yields the following.

1
2
3
4
5
6
7
(defn gregorian-to-julian
  [year month day hour minute second]
  (let [a (/ (- 14 month) 12),
        y (+ year 4800 (- a)),
        m (+ month (* 12 a) (- 3)),
        jdn (+ day (/ (+ (* 153 m) 2) 5) (* 365 y) (/ y 4) (- (/ y 100)) (/ y 400) (- 32045))]
    (+ jdn (/ (- hour 12) 24) (/ minute 1440) (/ second 86400))))

Because this was intended to be quick and dirty, I just manually set up a simple data structure for my earthquake data, and rather than spend lots of time integrating my Clojure code with the C library that implemented the lunar model, I simply used the C code to compute the distances, and then copied them directly into Clojure as well.

1
2
3
4
5
6
7
8
9
(def *quakes*
     (list [[1960 5 22] 9.5] [[2004 12 26] 9.3] [[1964 3 27] 9.2]
           [[1952 11 4] 9.0] [[1868 8 13] 9.0] [[2011 3 11] 8.9]
           [[1700 1 26] 8.9] [[2010 2 27] 8.8] [[1906 1 31] 8.8]
           [[1833 11 25] 8.8] [[1965 2 4] 8.7] [[1755 11 1] 8.7]
           [[1730 7 8] 8.7] [[2005 3 28] 8.6] [[1957 3 9] 8.6]
           [[1950 8 15] 8.6] [[2007 9 12] 8.5] [[1963 10 13] 8.5]
           [[1938 2 1] 8.5] [[1923 2 3] 8.5] [[1922 11 11] 8.5]
           [[1751 5 24] 8.5] [[1687 10 20] 8.5] [[1575 12 16] 8.5]))

To get the Julian dates for each quake, it’s sufficient to just run a quick and dirty map over the quake data as follows.

1
2
3
4
(defn build-julian-dates
  [quake-list]
  (let [gregdates (map first quake-list)]
    (map #((comp int gregorian-to-julian) (nth % 0) (nth % 1) (nth % 2) 11 0 0) gregdates)))

Now we feed that data into the C code from Github to calculate the distance the moon was from the earth at noon on the day of the quake (I made that approximation rather than look up the precise time each quake occurred). After those distances have been computed, I just copied them back into a Clojure data structure.

1
2
3
4
5
(def *distances* (list 401457.7391 406456.7146 397480.5329 391134.7747
                       367611.4645 384653.3206 391978.4295 357983.2221 404404.0269
                       388230.3723 395956.3832 360256.3522 401928.1778 375390.8061
                       368042.0999 371760.8433 403062.0740 399405.1708 395555.1038
                       364759.5125 373865.8605 361623.8578 366941.6539 401339.3053))

With that, we have enough data to plot earthquake strength versus lunar distance for the 24 largest quakes in recorded history. Clojure has a really nice library available called Incanter that makes charting and statistical analysis very easy.

1
2
3
4
5
6
7
8
9
(defn make-plot
  []
  (let [plot (scatter-plot (map second *quakes*) *distances*
                           :title "Quake Strength vs. Lunar Distance"
                           :y-label "Distance to Moon in km"
                           :x-label "Magnitude")]
    (set-y-range plot 356380 406710) ;; min/max lunar distance between 1575 and 2011
    (view plot)
    (save plot "plot.png")))

The resulting plot is shown below. As you can see, there’s no relationship at all between the strength of the strongest quakes and the current distance from the Earth to the moon. Note that the Y-axis is scaled to the minimum and maximum distances the moon has been in about 450 years, so the quakes really are distributed all over the place.

One last thing. Although the graph is pretty obvious, we really shouldn’t try and use intuition as a substitute for actual analysis, so let’s go ahead and fit a linear model to the data. Again, using Incanter,

1
(:r-square (linear-model *distances* (map second *quakes*)))

we get an value of 0.1193, or basically no relationship at all.