"""
Sun-earth colision time
jap, 11/2006
"""
from math import *


def g(r):
    "Acceleration of gravity, using units where G=1, Msun=1"
    return -1./r**2

t,dt = 0.,1E-4
rsol = 0.0047 # Approximate value for sun radius, in AU units
r,v = 1.0,0


def days(t):
    """
    Converts the time unit conforming to the choice of other units above into
    days. Note: If G=1, Msun=1, the earth orbital period around the
    sun is 2*pi time units, which mens that 1 time unit = 365.257/(2*pi) days!
    """
    return 365.257*t/(2*pi)


while r > rsol:
    a = g(r)
    v += a*dt
    r += v*dt
    t += dt

print "Collision time",days(t)
