# Code for debunking a claim made in the video at http://www.youtube.com/watch?v=95zMdTvoqcQ
# The video was posted on March 8, 2011 and contains a prediction of a seismic event
# around March 11, 2011 (which turned out to correspond to an actual earthquake in Japan).
#
# The specific claim here is that both during the Chile earhquake on 2010/02/17 and
# the "predicted" earthquake on 2011/03/11, there is an "alignment" between Sun, Earth
# and comet C/2010 X1 (Elenin). This can be restated as follows:
#
# "There exists a magic angle between comet Elenin and the Sun (as seen from Earth)
#  when a seismic events occur"
#
# To verify this claim, we take a list of visible comets from the Minor Planet Center
# ( http://www.minorplanetcenter.org/iau/Ephemerides/Comets/Soft03Cmt.txt )
# and for each one, calculate the delta between Sun-comet angles on 2010/02/17
# and 2011/03/11.
#
# Output format:
#   delta, magic angle, comet name
#
# See also: http://www.bautforum.com/showthread.php/113532-Comet-Elenin-Cause-Of-Earthquakes
#
# This code is placed in public domain.
#
# Requires pyephem library for ephemeris calculation.
#

import ephem
import math

def getSeparation(date, object1, object2):
    object1.compute(date)
    object2.compute(date)
    return ephem.separation(object1, object2)

sun = ephem.Sun()

date1 = '2010/02/17'
date2 = '2011/03/11'

f = open('Soft03Cmt.txt')
hits = 0
comets = 0
for line in f:
    if line.find('#') == 0:
        continue

    comets = comets+1

    comet = ephem.readdb(line)
    sep1 = getSeparation(date1, sun, comet)
    sep2 = getSeparation(date2, sun, comet)
    delta = 180*(sep1-sep2)/math.pi

    # We only list comets with delta < 6 degrees, since in the original claim
    # there is actually a 5.7 degree difference in angle between the earthquakes.
    # I mean, Elenin should be the only comet for which it works, right? ;-)
    if abs(delta) < 6:
        hits = hits + 1
        print "%5.2f %5.1f %s" % ( abs(delta), 180*sep1/math.pi, comet.name )


print "%d comets, %d hits (%5.2f %%)" % (comets, hits, 100.0*hits/comets)