Earth Moon ISS System

RK8.txt

+C_Coef = [0, .5, .5, .5 + np.sqrt(21)/ 14, .5 + np.sqrt(21)/ 14, .5, .5 - np.sqrt(21)/ 14, .5 - np.sqrt(21)/ 14, .5, .5 + np.sqrt(21)/ 14, 1]
+
+B_Coef = [.05, 0, 0, 0, 0, 0, 0, 49/180, 16/45, 49/180, .05]
+
+a = np.zeros((12,12))
+a[2,1]=1/2
+a[3,1]=1/4
+a[3,2]=1/4
+a[4,1]=1/7
+a[4,2]=-1/14-3/98*21**(1/2)
+a[4,3]=3/7+5/49*21**(1/2)
+a[5,1]=11/84+1/84*21**(1/2)
+a[5,2]=0
+a[5,3]=2/7+4/63*21**(1/2)
+a[5,4]=1/12-1/252*21**(1/2)
+a[6,1]=5/48+1/48*21**(1/2)
+a[6,2]=0
+a[6,3]=1/4+1/36*21**(1/2)
+a[6,4]=-77/120+7/180*21**(1/2)
+a[6,5]=63/80-7/80*21**(1/2)
+a[7,1]=5/21-1/42*21**(1/2)
+a[7,2]=0
+a[7,3]=-48/35+92/315*21**(1/2)
+a[7,4]=211/30-29/18*21**(1/2)
+a[7,5]=-36/5+23/14*21**(1/2)
+a[7,6]=9/5-13/35*21**(1/2)
+a[8,1]=1/14
+a[8,2]=0
+a[8,3]=0
+a[8,4]=0
+a[8,5]=1/9-1/42*21**(1/2)
+a[8,6]=13/63-1/21*21**(1/2)
+a[8,7]=1/9
+a[9,1]=1/32
+a[9,2]=0
+a[9,3]=0
+a[9,4]=0
+a[9,5]=91/576-7/192*21**(1/2)
+a[9,6]=11/72
+a[9,7]=-385/1152-25/384*21**(1/2)
+a[9,8]=63/128+13/128*21**(1/2)
+a[10,1]=1/14
+a[10,2]=0
+a[10,3]=0
+a[10,4]=0
+a[10,5]=1/9
+a[10,6]=-733/2205-1/15*21**(1/2)
+a[10,7]=515/504+37/168*21**(1/2)
+a[10,8]=-51/56-11/56*21**(1/2)
+a[10,9]=132/245+4/35*21**(1/2)
+a[11,1]=0
+a[11,2]=0
+a[11,3]=0
+a[11,4]=0
+a[11,5]=-7/3+7/18*21**(1/2)
+a[11,6]=-2/5+28/45*21**(1/2)
+a[11,7]=-91/24-53/72*21**(1/2)
+a[11,8]=301/72+53/72*21**(1/2)
+a[11,9]=28/45-28/45*21**(1/2)
+a[11,10]=49/18-7/18*21**(1/2)
+A_Coef = a[1:, 1:]
+
+def k_func(k_step, func, t, state, delT, A_Coef, B_Coef, C_Coef, *args, **kwargs):
+    """RK integration recursive function.
+    
+    Args:
+        k_step: the current kfunction sub-step.
+        func: The ODE function to integrate.
+        t: The time of the initial condition state.
+        state: The value of the initial condition at t.
+        delT: The time step to integrate across.
+        A_Coef: The first set of Butcher Tableau Coefficients.
+        B_Coef: The second set of Butcher Tableau Coefficients.
+        C_Coef: The third set of Butcher Tableau Coefficients.
+        *args: Supplemental arguments for the ODE function func.
+        **kwargs: Supplemental argumentsfor ODE function func.
+        
+    Returns:
+        The 'k_step'th k value for the integration step.
+    """
+    
+    if k_step == 1: #call the derivative function when necessary
+        return func(t, state, *args, **kwargs)
+    else:
+        carry = 0 #otherwise build the k value as a sum of the previous k values and thier coefficients recursively
+        for i in range(1, k_step):
+            carry += A_Coef[k_step - 1, i-1] * k_func(i-1, func, t, state, delT, A_Coef, B_Coef, C_Coef, *args, **kwargs)
+        return func(t + C_Coef[k_step - 1] * delT, state + delT*carry, *args, **kwargs)
+
+def rk_step(func, t, state, delT, *args, **kwargs):
+    """Generalized Runge-Kutta integrator
+
+    Args:
+        func: The ODE function to integrate.
+        t: The time of the initial condition state.
+        state: The value of the initial condition at t.
+        delT: The time step to integrate across.
+        A_Coef: The first set of Butcher Tableau Coefficients.
+        B_Coef: The second set of Butcher Tableau Coefficients.
+        C_Coef: The third set of Butcher Tableau Coefficients.
+        *args: Supplemental arguments for the ODE function func.
+        **kwargs: Supplemental argumentsfor ODE function func.
+
+    Returns:
+        Updated value of state at time t + delT.
+    """
+    sum = 0 # sum of k functions and their coefficients.
+    for k_step in range(1, A_Coef.shape[1] + 1): #For each step according to the Butcher Tableau
+        sum += B_Coef[k_step - 1] * k_func(k_step, func, t, state, delT, A_Coef, B_Coef, C_Coef, *args, **kwargs) #Add contribution of each integration step
+    return state + delT * sum

test.py

@@ -2,27 +2,40 @@ import numpy as np
 from core_functions import *
 
 
-masses = [5.972e24, 419725]
+masses = [5.972e24, 419725, 7.34767309e22]
 
-u_0 = np.array([[0,0,0,0,0,0], [400000 + 6378100, 0, 0, 0, 7823.6, 0]])
+u_0 = np.array([[0,0,0,0,0,0], [414000 + 6378100, 0, 0, 0, 7662, 0], [363300e3, 0, 0, 0, 1081, 0]])
+u_0 = np.array([[0,0,0,0,0,0], [4214929.7, 0, 5326067, 0, 7662, 0], [361836242.8, 0, 32579493.6, 0, 1081, 0]])
 
-T = 5400
+T = 2360600
 
 delta_t = 10
 
 u, times = ivp_RK4(u_0, T, delta_t, masses)
 
-rx = []
-ry = []
-rz = []
+rxsat = []
+rysat = []
+rzsat = []
 
 for i in range(len(times)):
-    rx.append(u[i][1][0])
-    ry.append(u[i][1][1])
-    rz.append(u[i][1][2])
+    rxsat.append(u[i][1][0])
+    rysat.append(u[i][1][1])
+    rzsat.append(u[i][1][2])
+
+rxmoon = []
+rymoon = []
+rzmoon = []
+
+for i in range(len(times)):
+    rxmoon.append(u[i][2][0])
+    rymoon.append(u[i][2][1])
+    rzmoon.append(u[i][2][2])
 
 import matplotlib.pyplot as plt
 
-plt.plot(rx,ry)
+plt.rcParams['figure.figsize'] = [20, 5]
+plt.plot(rxmoon, rymoon, label = 'Moon')
+plt.plot(rxsat, rysat, label = 'Satellite')
 
+plt.legend()
 plt.show()