Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- from numpy import sqrt, array, cos, sin, pi, arange, array
- from pylab import plot, show, legend, xlabel, ylabel
- import time
- %matplotlib inline
- start_time = time.time()
- M = 1.989e30
- G = 6.674e-11
- # This program solves dr/dt = f(r,t)
- def f(r,t):
- # unpack variables
- x = r[0]
- y = r[1]
- vx = r[2]
- vy = r[3]
- radiusMagnitude = sqrt(x**2 + y**2)
- fx = vx
- fy = vy
- fvx = -G * M * x / (radiusMagnitude**3)
- fvy = -G * M * y / (radiusMagnitude**3)
- return array([fx, fy, fvx, fvy], float)
- # initial conditions
- t0 = 0.0 # initial time
- tf = 96 * 3.154e7 # 3.154e7 is # of seconds in years
- N = 75000
- h = (tf - t0) / N # timestep
- tpoints = arange(t0, tf, h) # array of time values
- x = 4e12 # initial x in meters
- y = 0.0 # intiial y
- vx = 0 # initial vx
- vy = 500 # initial vy m/s
- r = array([x, y, vx, vy], float) # initial r vector
- xpoints = []
- ypoints = []
- for t in tpoints:
- xpoints.append(r[0])
- ypoints.append(r[1])
- k1 = h * f(r, t)
- k2 = h * f(r + 0.5 * k1, t + 0.5 * h)
- k3 = h * f(r + 0.5 * k2, t + 0.5 * h)
- k4 = h * f(r + k3, t + h)
- r += 1/6 * (k1 + 2 * k2 + 2 * k3 + k4)
- plot (ypoints, xpoints, "b-", label="y vs x")
- xlabel("y pos")
- ylabel("x pos")
- legend()
- show()
- print("--- %s seconds taken to run ---" % (time.time() - start_time))
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement