Merge branch 'adbeck2' into 'main'

Created ISS Example

See merge request !2

core_functions.py

 import numpy as np
 
-def f_true(r1, r2, m2):
-
+# True Acceleration Function
+def f_true(u_ki, index, u_k, masses):
     G = 6.6743e-11
+    r1 = np.array([u_ki[0], u_ki[1], u_ki[2]])
+    
+    a = np.zeros(3)
+    for h in range(len(u_k)):
+        if h == index:
+            continue
+        m2 = masses[h]
+        r2 = np.array([u_k[h][0], u_k[h][1], u_k[h][2]])
+        a_temp = G*m2*(r2-r1) / (np.linalg.norm(r2-r1)**3)
+        a += a_temp
+
+    xdot = np.array([u_ki[3], u_ki[4], u_ki[5], a[0], a[1], a[2]])
+
+    return xdot
+
+# MUST CHANGE TO CALCULATE FOR EACH DISCRETE OBJECT
+def RK4(u_ki, index, u_k, delta_t, masses):
+
+    k1 = delta_t * f_true(u_ki, index, u_k, masses)
+    k2 = delta_t * f_true(u_ki + k1/2, index, u_k, masses)
+    k3 = delta_t * f_true(u_ki + k2/2, index, u_k, masses)
+    k4 = delta_t * f_true(u_ki + k3, index, u_k, masses)
+        
+    return u_ki + (1/6) * (k1 + 2*k2 + 2*k3 + k4)
 
-    a = G*m2*(r2-r1) / (np.linalg.norm(r2-r1)**3)
+def ivp_RK4(u_0, T, delta_t, masses):
+    u = [u_0]
+    
+    times = np.arange(0,T+delta_t,delta_t)
+    
+    for i in range(int(T/delta_t)):
+        previous = u[i]
 
-    return a
+        update = []
+        for j in range(len(previous)):
+            update_temp = RK4(previous[j], j, previous, delta_t, masses)  #UPDATES EACH OBJECT
+            update.append(update_temp)
 
-    
\ No newline at end of file
+        update = np.array(update)
+        u.append(update)
+        
+    u = np.array(u)
+    
+    return u, times
+  
\ No newline at end of file

test.py

+import numpy as np
+from core_functions import *
+
+
+masses = [5.972e24, 419725]
+
+u_0 = np.array([[0,0,0,0,0,0], [400000 + 6378100, 0, 0, 0, 7823.6, 0]])
+
+T = 5400
+
+delta_t = 10
+
+u, times = ivp_RK4(u_0, T, delta_t, masses)
+
+rx = []
+ry = []
+rz = []
+
+for i in range(len(times)):
+    rx.append(u[i][1][0])
+    ry.append(u[i][1][1])
+    rz.append(u[i][1][2])
+
+import matplotlib.pyplot as plt
+
+plt.plot(rx,ry)
+
+plt.show()