Skip to content

Commit

Permalink
Updated Asteroid_trajectory.py, works with the asteroid 2007 VW266
Browse files Browse the repository at this point in the history
  • Loading branch information
@Keith-Maxwell
Keith-Maxwell committed May 8, 2020
1 parent f0ffffc commit 7f86963
Show file tree
Hide file tree
Showing 4 changed files with 79 additions and 133 deletions.
195 changes: 70 additions & 125 deletions Asteroid_trajectory.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,15 +4,14 @@


class planet(object):
def __init__(self, m, a, i, e, Omega, omega, M0, t0=0):
def __init__(self, m, a, i, e, Omega, omega, M0):
self.m = m
self.a = a
self.i = np.radians(i)
self.e = e
self.Omega = np.radians(Omega)
self.omega = np.radians(omega)
self.M0 = np.radians(M0)
self.t0 = t0
self.T = np.sqrt((4 * np.pi ** 2) / (G * (m_sun + self.m)) * self.a ** 3)
self.w = 2 * np.pi / self.T

Expand All @@ -28,7 +27,7 @@ def orbitalparam2vectorList(self, timevector):

def completeOrbitalElem2Vector(self, t):
self.n = k / np.sqrt(self.a ** 3)
self.M = self.M0 + self.n * (t - self.t0)
self.M = self.M0 + self.n * (t)
self.E = self.newton(self.keplerEquation, self.M0)
self.bigX = self.a * (np.cos(self.E) - self.e)
self.bigY = self.a * np.sqrt(1 - self.e ** 2) * np.sin(self.E)
Expand All @@ -44,7 +43,8 @@ def completeOrbitalElem2Vector(self, t):
self.rotation1(-self.i)),
self.rotation3(-self.omega)),
np.array([[self.bigXdot], [self.bigYdot], [0]]))
return [self.position[0, 0], self.position[1, 0], self.position[2, 0]]
return [self.position[0, 0], self.position[1, 0], self.position[2, 0],
self.velocity[0, 0], self.velocity[1, 0], self.velocity[2, 0]]

def rotation1(self, theta):
return np.array([[1, 0, 0],
Expand All @@ -67,106 +67,76 @@ def keplerEquation(self, E):
return E - self.e * np.sin(E) - self.M


class twoBody(): # Only the Sun exerts it's influence on the body.
@classmethod
def func(cls, t, y): # using the state space representation of the equation of motion
y0 = y[3] # velocity on x
y1 = y[4] # velocity on y
y2 = y[5] # velocity on z
r = np.sqrt(y[0] ** 2 + y[1] ** 2 + y[2] ** 2) # norm of the vector r
y3 = -G * (m_sun + m_object) / (r ** 3) * \
y[0] # equation of motion on x
y4 = -G * (m_sun + m_object) / (r ** 3) * \
y[1] # equation of motion on y
y5 = -G * (m_sun + m_object) / (r ** 3) * \
y[2] # equation of motion on z
return np.array([y0, y1, y2, y3, y4, y5])

@classmethod
def RK4(cls, time_vector, initial_conditions, h):
# Definition of the Runge-Kutta method at the order 4
results = []
yin = initial_conditions
results.append(yin) # set the first value to the initial conditions
for t in time_vector[1:]:
k1 = cls.func(t, yin)
k2 = cls.func(t + h / 2, yin + h / 2 * k1)
k3 = cls.func(t + h / 2, yin + h / 2 * k2)
k4 = cls.func(t + h, yin + h * k3)
yin = yin + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
# each value calculated for given t is added to the list
results.append(yin)
return np.array(results)
class MultiBody():
''' This class groups together 3 functions :
- the Equation of Motion, which is the core of the program
- The Ruunge Kutta 4 integrator
- The Forward Backward method'''

@classmethod
def forward_backward(cls, time_vector, initial_conditions, h):
# allows for error computation. we just need to compute the difference at the starting point
# call RK4 in forward movement
y_forward = cls.RK4(time_vector, initial_conditions, h)
new_y0 = y_forward[-1] # new initial conditions
y_backward = cls.RK4(np.flip(time_vector),
new_y0, -h) # call RK4 backwards
return y_forward, y_backward

def func(cls, t, y, planets):
''' State space representation of the Equation of Perturbated Motion
Takes into account every selected planet.'''

class threeBody():
@classmethod
def func(cls, t, y, planet):
y0 = y[3] # velocity on x
y1 = y[4] # velocity on y
y2 = y[5] # velocity on z
r = np.sqrt(y[0] ** 2 + y[1] ** 2 + y[2] ** 2) # norm of the vector r
pos = planet.completeOrbitalElem2Vector(t)
# delta is the difference between r and position_Jupiter (distance between the asteroid and jupiter)
delta = np.sqrt((y[0] - pos[0]) ** 2 + (y[1] - pos[1]) ** 2 + (y[2] - pos[2]) ** 2)
# equation of motion on x
y3 = -G * (m_sun + m_object) / (r ** 3) * y[0] \
- G * planet.m * ((y[0] - pos[0]) / (delta ** 3) + pos[0] / np.linalg.norm(pos) ** 3)
# equation of motion on y
y4 = -G * (m_sun + m_object) / (r ** 3) * y[1] \
- G * planet.m * ((y[1] - pos[1]) / (delta ** 3) + pos[1] / np.linalg.norm(pos) ** 3)
# equation of motion on z
y5 = -G * (m_sun + m_object) / (r ** 3) * y[2] \
- G * planet.m * ((y[2] - pos[2]) / (delta ** 3) + pos[2] / np.linalg.norm(pos) ** 3)
y3 = -G * (m_sun + float(asteroid.m)) / (r ** 3) * y[0]
y4 = -G * (m_sun + float(asteroid.m)) / (r ** 3) * y[1]
y5 = -G * (m_sun + float(asteroid.m)) / (r ** 3) * y[2]
for planet in planets:
pos = planet.completeOrbitalElem2Vector(t)
# delta is the difference between 'r' and 'pos', it is the distance between asteroid & planet
delta = np.sqrt((y[0] - pos[0]) ** 2 + (y[1] - pos[1])
** 2 + (y[2] - pos[2]) ** 2)
# equation of motion on x
y3 += - G * planet.m * ((y[0] - pos[0]) / (delta ** 3) + pos[0] / np.linalg.norm(pos) ** 3)
# equation of motion on y
y4 += - G * planet.m * ((y[1] - pos[1]) / (delta ** 3) + pos[1] / np.linalg.norm(pos) ** 3)
# equation of motion on z
y5 += - G * planet.m * ((y[2] - pos[2]) / (delta ** 3) + pos[2] / np.linalg.norm(pos) ** 3)
return np.array([y0, y1, y2, y3, y4, y5])

@classmethod
def RK4(cls, time_vector, initial_conditions, h, planet):
# Definition of the Runge-Kutta method at the order 4
def RK4(cls, time_vector, initial_conditions, h, planets):
''' Runge kutta 4 integrator :
takes a time vector, an initial conditions vector, a step and a list of planets
as an input, and outputs the list of positions and accelerations'''

results = []
yin = initial_conditions
results.append(yin) # set the first value to the initial conditions
for i in range(1, len(time_vector[1:])):
k1 = cls.func(time_vector[i], yin, planet)
k2 = cls.func(time_vector[i] + h / 2, yin + h / 2 * k1, planet)
k3 = cls.func(time_vector[i] + h / 2, yin + h / 2 * k2, planet)
k4 = cls.func(time_vector[i] + h, yin + h * k3, planet)
for t in time_vector[1:]:
k1 = cls.func(t, yin, planets)
k2 = cls.func(t + h / 2, yin + h / 2 * k1, planets)
k3 = cls.func(t + h / 2, yin + h / 2 * k2, planets)
k4 = cls.func(t + h, yin + h * k3, planets)
yin = yin + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
# each value calculated for given t is added to the list
results.append(yin)
return np.array(results)

@classmethod
def forward_backward(cls, time_vector, initial_conditions, h, planet):
# allows for error computation. we just need to compute the difference at the starting point
y_forward = cls.RK4(time_vector, initial_conditions, h, planet) # call RK4 in forward movement
def forward_backward(cls, time_vector, initial_conditions, h, planets):
''' Calls the RK4 integrator a first time from the initial condition,
and a second time from the results, thus going back to the very
first starting point. By computing the difference between the results
and the initital conditions at the starting point, we obtain the
errors due to the integration'''

y_forward = cls.RK4(time_vector, initial_conditions, h, planets) # call RK4 forwards
new_y0 = y_forward[-1] # new initial conditions
y_backward = cls.RK4(np.flip(time_vector), new_y0, -h, planet) # call RK4 backwards
y_backward = cls.RK4(np.flip(time_vector), new_y0, - h, planets) # call RK4 backwards

return y_forward, y_backward


def extract_vectors(results):
# takes the list of outputs from RK4 and extracts a list of position and velocity vectors
r_list = []
r_dot_list = []
for i in range(len(results)):
r = []
r_dot = []
for j in range(0, 3):
r.append(results[i][j])
r_dot.append(results[i][j + 3])
r_list.append(r)
r_dot_list.append(r_dot)
'''takes the list of outputs from RK4 and extracts a list of position and velocity vectors'''

r_list = [[results[i][j] for j in range(0, 3)] for i in range(len(results))]
r_dot_list = [[results[i][j + 3] for j in range(0, 3)] for i in range(len(results))]
return r_list, r_dot_list


Expand Down Expand Up @@ -215,68 +185,43 @@ def plotOrbitalVariation(a, e, i, t):

# initial variables
m_sun = 1 # mass of the Sun
m_object = 0 # mass of the object studied
a = (5.2 ** 3 / (3) ** 2) ** (1 / 3) # semi-major axis
G = 0.000295824 # gravitation constant expressed in our own system of units
k = np.sqrt(G)
mu = G * m_sun

# definition of Jupiter orbital parameters
# Creation of the each planet from their orbital parameters at J2000
# mercury=planet(m_sun / 6023622.047, 0.38709893, 7.00487, 0.20563069, 48.33167, 77.45645, 126.464)
# venus=planet(m_sun / 408544.7263, 0.72333199, 3.39471, 0.00677323, 76.68069, 131.53298, -26.23394)
# earth=planet(m_sun / 332965.0462, 1, 0.00005, 0.01671022, 348.73936, 102.94719, -351.2222)
# mars=planet(m_sun / 3098854.873, 1.52366231, 1.85061, 0.09341233, 49.57854, 336.04084, -30.16606)
jupiter = planet(m_sun / 1047.348625, 5.202603, 1.303, 0.048498, 100.46, -86.13, 20.0)

# initial conditions of the asteroid [pos x, pos y, pos z, v x, v y, v z]
init_state = np.array([a, 0.0, 0.0, 0.0, k / np.sqrt(a), 0.0])
# saturn=planet(m_sun / 3498.926925, 9.5370703, 2.484, 0.054151, 113.72, 92.43, -156.2)
# uranus=planet(m_sun / 22906.30182, 19.19126393, 0.76986, 0.04716771, 74.22988, 170.96424, 68.0392)
# neptune=planet(, m_sun / 19418.13922, 30.06896348, 1.76917, 0.00858587, 131.72169, 44.97135, 128.19)

# integration parameters
tf = 40000 # final time
ti = 0 # starting time
tf = 50000 # Number of days
ti = 5055.5 # starting time (Epoch)
step = 1 # step wanted

# creation of the list containing each value of time
time = np.arange(ti, tf + step, step)

# =============================================================================================
# Start of the computation for 2 body problem
# =============================================================================================
'''
start_time = chrono.time()
forward, backward = twoBody.forward_backward(time, init_state, step)
# computation of the error between forward and backward
err_x = 150e6 * abs(forward[0][0] - backward[-1][0])
err_y = 150e6 * abs(forward[0][1] - backward[-1][1])
print(' Error on x : ', err_x, ' km\n', 'Error on y : ', err_y, ' km')
# TODO: Add automatic step adjustment in function of the error ?
time = np.arange(ti, tf + step + ti, step)

# plot of the trajectory
plt.plot(forward[:, 0], forward[:, 1], 'o', color='red', markersize=1)
plt.axis('equal')
plt.title('Unperturbed trajectory around the sun')
plt.show()
# get the position and velocity from previous results
#r_list, rdot_list = extract_vectors(forward)
# calculate the orbital parameters from the position and velocity
#a_list, e_list, i_list = orbital_parameters_list(r_list, rdot_list)
# plot the orbital parameters
#plotOrbitalVariation(a_list, e_list, i_list, time)
# initial conditions of the asteroid [pos x, pos y, pos z, v x, v y, v z]
asteroid = planet(0, 5.454, 108.358, 0.3896, 276.509, 226.107, 146.88)
init_state = np.array(asteroid.completeOrbitalElem2Vector(ti))

stop_time = chrono.time()
print("Computation time = ", stop_time-start_time, " seconds")
'''
# =============================================================================================
# Start of the computation for 3 body problem
# =============================================================================================
start_time = chrono.time()
forward2, backward2 = threeBody.forward_backward(time, init_state, step, jupiter)
forward2, backward2 = MultiBody.forward_backward(time, init_state, step, [jupiter])

# computation of the error between forward and backward
err_x2 = 150e6 * abs(forward2[0][0] - backward2[-1][0])
err_y2 = 150e6 * abs(forward2[0][1] - backward2[-1][1])
print(' Error on x : ', err_x2, ' km\n', 'Error on y : ', err_y2, ' km')
err_x = 150e6 * abs(forward2[0][0] - backward2[-1][0]) / 2
err_y = 150e6 * abs(forward2[0][1] - backward2[-1][1]) / 2
err_z = 150e6 * abs(forward2[0][2] - backward2[-1][2]) / 2
print(' Error on x : ', err_x, ' km\n', 'Error on y : ', err_y, ' km\n', 'Error on z : ', err_z, ' km')

# plot of the trajectory
jup_param = jupiter.orbitalparam2vectorList(time)
Expand All @@ -294,7 +239,7 @@ def plotOrbitalVariation(a, e, i, t):
a_list, e_list, i_list = orbital_parameters_list(r_list, rdot_list)

# plot the orbital parameters
plotOrbitalVariation(a_list, e_list, i_list, time[1:])
plotOrbitalVariation(a_list, e_list, i_list, time)

stop_time = chrono.time()
print("Computation time = ", stop_time - start_time, " seconds")
7 changes: 4 additions & 3 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ This project is part of my CIRI (Cours d’Initiation à la Recherche et à l’

The objective is to model the trajectory of an asteroid in the solar system subject to perturbations from Jupiter. We then study the variations of the orbital parameters such as the semi-major axis, the eccentricity and the inclination after a few centuries. This allows us to visualize the resonances.

We went further by adding all the other planets in the solar system, which we can choose from using a Graphical User Interface. The orbital parameters of the asteroid are user-defined, as well as the integration parameters.
We went further by adding all the other planets in the solar system, which we can choose using a Graphical User Interface. The orbital parameters of the asteroid are user-defined, as well as the integration parameters.

## How to use

Expand Down Expand Up @@ -47,7 +47,7 @@ If you want to re-run the simulation with other parameters, simply change the de
- Asteroid parameters :
- Circular orbit
- User defined orbit with parameters
- User defined orbit with position and velocity
- User defined orbit with position and velocity (bugs)
- Plot of the trajectory of every selected body in 3d
- Plot of the orbital parameters changes due to perturbation (only if a planet has been selected)

Expand All @@ -58,7 +58,8 @@ If you want to re-run the simulation with other parameters, simply change the de
- Option to choose from simplified or accurate planet trajectory model
- more integrators
- input proofing and displaying errors instead of crashing the program
- Logging
- Saving the results to a file for use in Excel for example
- Improve Performance
- re-write some functions in C++ to improve performance

### Contributions are welcome !
10 changes: 5 additions & 5 deletions main_GUI.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -509,9 +509,9 @@ def Start(self):
self.err_x = 150e6 * abs(self.forward[0][0] - self.backward[-1][0]) / 2
self.err_y = 150e6 * abs(self.forward[0][1] - self.backward[-1][1]) / 2
self.err_z = 150e6 * abs(self.forward[0][2] - self.backward[-1][2]) / 2
self.fwbwOutputLabel.setText(f"Error on x: \n{str(round(self.err_x, 2))} km\n\n \
Error on y: \n{str(round(self.err_y, 2))} km\n\n \
Error on z: \n{str(round(self.err_z, 2))} km")
self.fwbwOutputLabel.setText("Error on x: \n" + str(round(self.err_x, 2)) +
" km\n\nError on y: \n" + str(round(self.err_y, 2)) +
"km\n\nError on z: \n" + str(round(self.err_z, 2)) + "km")
self.results = self.forward # Store the results from the forward integration in a dedicated list

else: # if the user did not ask for a fw-bw integration, then only compute the results.
Expand Down Expand Up @@ -715,7 +715,7 @@ def get_init_state(self, t=0):
initial position and velocity vector. All the values are read from the inputs.'''

if ui.circularCheckBox.checkState(): # Circular initial orbit
self.a = float(ui.inputAsteroidSmAxis.text())
self.a = float(ui.inputAsteroidSMaxis.text())
self.init_state = np.array([self.a, 0, 0,
0, k / np.sqrt(self.a), 0])
return self.init_state
Expand All @@ -729,7 +729,7 @@ def get_init_state(self, t=0):
return self.init_state

else: # defined by orbital parameters
self.a = float(ui.inputAsteroidSmAxis.text())
self.a = float(ui.inputAsteroidSMaxis.text())
self.i = np.radians(float(ui.inputInc.text()))
self.e = float(ui.inputEcc.text())
self.Omega = np.radians(float(ui.inputLan.text()))
Expand Down
Binary file modified screen.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.

0 comments on commit 7f86963

Please sign in to comment.