#Evaluation of error in Euler method, Prof. Ellen Williams, 9/13/01

#Minor modifications, TLE, 9/03

from visual import *

#G is gravitational constant in N-m**2/kg**2, Re is earth radius in m

#Me is earth mass in kg

#radially falling object with acceleration -G*Me/r**2

G=6.67e-11

Re=6.37e6

Me=5.98e24

#Evaluate position and velocity vs time, h and v vs t for object falling from h0

#h0 and h in m, v in m/s, set initial conditions

#define h, v, t, and a as lists each with one initial element

h0=2.5e6

h=[h0+Re]

v=[0]

t=[0]

a=[-G*Me/h[0]**2]

#print a[0]

# dt is the time increment for Euler's eqn, units seconds

#repeat calculation with variable value of dt and evaluate change in t and v

dt=100.0

icount=0

#print " time height velocity accel'n"

#add this once the program works

while h[-1]>Re:

,,print icount

,,icount = icount+1

#Euler method----Note the order; would give different results if changed

,,h.append(h[-1]+v[-1]*dt) # this is h0 +v0*dt = h1, etc.

,,v.append(v[-1]+a[-1]*dt) # this is v0 +a0*dt = v1, etc.

,,a.append(-G*Me/h[-1]**2) # this is -G*Me/h1**2 = a1, etc.

,,t.append(t[-1]+dt) # t1 = t0 + dt, etc.

,,print "%10.0f %10.4E %10.4E %10.4E" % (t[-1], h[-1], v[-1], a[-1])

print " "

print "The array length and counter end values are %d and %d" % (len(h), count)

print "The value of dt is %10.2f" % (dt)

print "The time at which Re was reached "

print "%10.1f" % (t[-2] + (t[-1]-t[-2])*(h[-2]-Re)/(h[-2]-h[-1]))

print "The value of v at Re" print "%10.4E" % (v[-2] + (v[-1]-v[-2])*(h[-2]-Re)/(h[-2]-h[-1]))

#analytical result for h0=2.5e6 is 5.94e3 m/s #analytical result for h0=2.5e5 is 2.17e3 m/s