diff --git a/Exercise 1/exercise1.py b/Exercise 1/exercise1.py index 66228e2..5c0e797 100644 --- a/Exercise 1/exercise1.py +++ b/Exercise 1/exercise1.py @@ -4,38 +4,40 @@ from scipy.integrate import solve_ivp import matplotlib.pyplot as plt -def part1(radFactor, orbits): +def part1(velocityFactor, orbits): #Part one Diffeq def f_part1(t, state, Me, Mm, G): - xm, ym, vx, vy = state + xm, ym, vx, vy = state #all input values - dvxdt = -(Me*G*xm)/((xm**2+ym**2)**(3/2)) - dxmdt = vx - dvydt = -(Me*G*ym)/((xm**2+ym**2)**(3/2)) + dxmdt = vx #functions for each diffeq dymdt = vy + dvxdt = -(Me*G*xm)/((xm**2+ym**2)**(3/2)) + dvydt = -(Me*G*ym)/((xm**2+ym**2)**(3/2)) + - return (dxmdt, dymdt, dvxdt, dvydt) + return (dxmdt, dymdt, dvxdt, dvydt) #return differentiated values # Initial Conditions - Me = 5.97*(10**24) + Me = 5.97*(10**24) #Constants Mm = 7.35*(10**22) G = 6.67*(10**-11) t_min = 0 - t_max = 2360620*int(orbits) - numpoints = 2000*int(orbits) + t_max = 2360620*int(orbits) #Time for one orbit * number of orbits + numpoints = 2000*int(orbits) #To keep good clarity when plotting t = np.linspace(t_min, t_max, numpoints) - r = 384400000 - v = sqrt((G*Me)/r)*float(radFactor) + r = 384400000 #Distance or orbit + v = sqrt((G*Me)/r)*float(velocityFactor)#Balancing centrepetal force and gravitational force, + #then multiplying by a scalar to obtain an elliptical orbit - xm0 = r - ym0 = 0 + xm0 = r #Initial position and velocity of moon + ym0 = 0 #Set so moon is travelling only in the y direction initially vx0 = 0 vy0 = v - rtol = 1e-5 + rtol = 1e-5 #Tolerences set to balance computation speed with accuracy atol = 1e-8 @@ -58,7 +60,7 @@ def part1(radFactor, orbits): ax.axvline(x=0, color='k', linewidth=0.5) plt.show() -def part2(radFactor, orbits): +def part2(velocityFactor, orbits): #Part one Diffeq def f_part1(t, state, Me, Mm, G): xm, ym, vx, vy, xp, yp, vpx, vpy = state @@ -90,7 +92,7 @@ def part2(radFactor, orbits): t = np.linspace(t_min, t_max, numpoints) rm = 384400000 - vm = sqrt((G*Me)/rm)*float(radFactor) + vm = sqrt((G*Me)/rm)*float(velocityFactor) rpm = 10000000 vpm = sqrt((G*Mm)/rpm) @@ -108,7 +110,6 @@ def part2(radFactor, orbits): rtol = 1e-6 atol = 1e-9 - #Solver results = solve_ivp(f_part1, (t_min,t_max), (xm0, ym0, vx0, vy0, xp0, yp0, vpx0, vpy0), args=(Me, Mm, G), t_eval=t, atol=atol, rtol=rtol) @@ -134,14 +135,14 @@ while MyInput != 'q': print('You entered the choice: ',MyInput) if MyInput == '1': print('You have chosen part (1): simulation of a lunar orbit') - radFactor = input(f'Enter the velocity scalar. A value of 1 will result in a circular orbit, and anything else will give an elliptical orbit: ') + velocityFactor = input(f'Enter the velocity scalar. A value of 1 will result in a circular orbit, and anything else will give an elliptical orbit: ') orbits = input(f'Enter the approximate number of orbits desired: ') - part1(radFactor, orbits) + part1(velocityFactor, orbits) elif MyInput == '2': print('You have chosen part (2): earth-moon-probe system') - radFactor = input(f'Enter the velocity scalar. A value of 1 will result in a circular orbit, and anything else will give an elliptical orbit: ') + velocityFactor = input(f'Enter the velocity scalar. A value of 1 will result in a circular orbit, and anything else will give an elliptical orbit: ') orbits = input(f'Enter the approximate number of orbits desired: ') - part2(radFactor, orbits) + part2(velocityFactor, orbits) elif MyInput != 'q': print('This is not a valid choice') print('You have chosen to finish - goodbye.') \ No newline at end of file