Advertisement
Guest User

Untitled

a guest
Apr 7th, 2020
174
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 1.45 KB | None | 0 0
  1. from numpy import sqrt, array, cos, sin, pi, arange, array
  2. from pylab import plot, show, legend, xlabel, ylabel
  3. import time
  4.  
  5. %matplotlib inline
  6.  
  7. start_time = time.time()
  8.  
  9. M = 1.989e30
  10. G = 6.674e-11
  11.  
  12. # This program solves dr/dt = f(r,t)
  13. def f(r,t):
  14.    
  15.     # unpack variables
  16.     x = r[0]
  17.     y = r[1]
  18.     vx = r[2]
  19.     vy = r[3]
  20.    
  21.     radiusMagnitude = sqrt(x**2 + y**2)
  22.    
  23.     fx = vx
  24.     fy = vy
  25.     fvx = -G * M * x / (radiusMagnitude**3)
  26.     fvy = -G * M * y / (radiusMagnitude**3)
  27.    
  28.     return array([fx, fy, fvx, fvy], float)
  29.  
  30. # initial conditions
  31. t0 = 0.0                          # initial time
  32. tf = 96 * 3.154e7 # 3.154e7 is # of seconds in years
  33. N = 75000
  34. h = (tf - t0) / N # timestep
  35. tpoints = arange(t0, tf, h) # array of time values
  36.  
  37. x = 4e12                          # initial x in meters
  38. y = 0.0                          # intiial y
  39. vx = 0                # initial vx
  40. vy = 500       # initial vy m/s
  41. r = array([x, y, vx, vy], float) # initial r vector
  42. xpoints = []
  43. ypoints = []
  44.  
  45. for t in tpoints:
  46.     xpoints.append(r[0])
  47.     ypoints.append(r[1])
  48.    
  49.     k1 = h * f(r, t)
  50.     k2 = h * f(r + 0.5 * k1, t + 0.5 * h)
  51.     k3 = h * f(r + 0.5 * k2, t + 0.5 * h)
  52.     k4 = h * f(r + k3, t + h)
  53.     r += 1/6 * (k1 + 2 * k2 + 2 * k3 + k4)
  54.    
  55. plot (ypoints, xpoints, "b-", label="y vs x")
  56. xlabel("y pos")
  57. ylabel("x pos")
  58. legend()
  59. show()
  60.  
  61. print("--- %s seconds taken to run ---" % (time.time() - start_time))
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement