Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- def fx(x):
- r = p / (1 + e*np.cos(x))
- return r
- import numpy as np
- import matplotlib.pyplot as plt
- pi, twopi = np.pi, 2*np.pi
- degs, rads = 180/pi, pi/180
- # Setup
- L = 9.11E+38 #L = angular momentum
- m = 3.28E+23 #m = mass of mercury
- M = 1.99E+30 #M = mass of sun
- a = 5.8E+10 #a = semi-major axis
- G = 6.674E-11 #G = Gravitationl constant
- k = G*M*m
- E = -k/(2*a) #E = energy
- #E = 0
- p = L**2 / (m*k)
- c = 1 + (2*E*L**2) / (m*k**2)
- e = np.sqrt(c) #e = eccentricity
- print(e,c)
- # Initialize
- phi = np.linspace(0, twopi, 100)
- radius = fx(phi)
- theta = degs * phi
- x, y = radius * np.cos(phi), radius * np.sin(phi)
- print('r =',radius)
- print('x =',x)
- print('y =',y)
- plt.plot(x,y)
- plt.show()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement