From ea1b5a04f7d47e5f65d9faed1909688f3a05e98a Mon Sep 17 00:00:00 2001
From: Ceres <ceres@ceressees.dev>
Date: Mon, 13 Oct 2025 18:22:23 +0100
Subject: [PATCH] Added part 2

---
 Exercise 1/exercise1.py | 75 +++++++++++++++++++++++++++++++++++++++--
 1 file changed, 72 insertions(+), 3 deletions(-)

diff --git a/Exercise 1/exercise1.py b/Exercise 1/exercise1.py
index 274230a..66228e2 100644
--- a/Exercise 1/exercise1.py	
+++ b/Exercise 1/exercise1.py	
@@ -58,6 +58,75 @@ def part1(radFactor, orbits):
     ax.axvline(x=0, color='k', linewidth=0.5)
     plt.show()
 
+def part2(radFactor, orbits):
+    #Part one Diffeq
+    def f_part1(t, state, Me, Mm, G):
+        xm, ym, vx, vy, xp, yp, vpx, vpy = state
+
+        xpm = xp-xm
+        ypm = yp-ym
+
+        dxmdt = vx
+        dymdt = vy
+        dvxdt = -(Me*G*xm)/((xm**2+ym**2)**(3/2))
+        dvydt = -(Me*G*ym)/((xm**2+ym**2)**(3/2))
+        
+        dxpdt = vpx
+        dypdt = vpy
+        dvpxdt = -((Me*G*xp)/((xp**2+yp**2)**(3/2)))-((Mm*G*xpm)/((xpm**2+ypm**2)**(3/2)))
+        dvpydt = -((Me*G*yp)/((xp**2+yp**2)**(3/2)))-((Mm*G*ypm)/((xpm**2+ypm**2)**(3/2)))
+        
+        return (dxmdt, dymdt, dvxdt, dvydt, dxpdt, dypdt, dvpxdt, dvpydt)
+
+
+    # Initial Conditions
+    Me = 5.97*(10**24)
+    Mm = 7.35*(10**22)
+    G = 6.67*(10**-11)
+
+    t_min = 0
+    t_max = 2360620*int(orbits)
+    numpoints = 2000*int(orbits)
+    t = np.linspace(t_min, t_max, numpoints)
+
+    rm = 384400000
+    vm = sqrt((G*Me)/rm)*float(radFactor)
+
+    rpm = 10000000
+    vpm = sqrt((G*Mm)/rpm)
+    
+    xm0 = rm
+    ym0 = 0
+    vx0 = 0
+    vy0 = vm
+
+    xp0 = rpm+rm
+    yp0 = 0
+    vpx0 = 0
+    vpy0 = vm+vpm
+
+    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)
+
+
+    #Graph plotting
+    ax=plt.axes()    # This creates some axes, so that we
+    ax.set_aspect(1) # can set the aspect ratio to 1 i.e. 
+                    # x and y axes are scaled equally.                                
+    ax.set_xlabel("x coordinate (m)") # Must label axes (with
+    ax.set_ylabel("y coordinate (m)") # units) and give
+    ax.set_title("Orbit of moon around earth") # plot title.
+    ax.plot(results.y[0],results.y[1], label="Moon" ) # Make the plot
+    ax.plot(results.y[4],results.y[5], label="Probe")
+    ax.legend(loc='upper right')              # and add a key.
+    ax.axhline(y=0, color='k', linewidth=0.5)
+    ax.axvline(x=0, color='k', linewidth=0.5)
+    plt.show()
+
 
 MyInput = '0'
 while MyInput != 'q':
@@ -70,9 +139,9 @@ while MyInput != 'q':
         part1(radFactor, orbits)
     elif MyInput == '2':
         print('You have chosen part (2): earth-moon-probe system')
-        #
-        # put your code for part (2) here
-        #
+        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: ')
+        orbits = input(f'Enter the approximate number of orbits desired: ')
+        part2(radFactor, orbits)
     elif MyInput != 'q':
         print('This is not a valid choice')
 print('You have chosen to finish - goodbye.')
\ No newline at end of file