Merge branch 'main' of https://git.ceressees.dev/ceres-sees-all/Uni-Computing

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
 
-        dxmdt = vx
+        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.')