"""
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)