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