Skip to content
Permalink
master
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?

Project-2/Project 2.ipynb

Go to file
 
 
Cannot retrieve contributors at this time
1468 lines (1468 sloc) 56.5 KB
Raw Blame
Open in GitHub Desktop
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# ME 3295-003: Orbital Mechanics - Project 2, Problem 1\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Stephanie Rodrigue and James Jordan"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from scipy.optimize import newton\n",
"from scipy.integrate import solve_ivp\n",
"import numpy as np\n",
"import math\n",
"from astropy.time import Time\n",
"from datetime import timedelta, time\n",
"from datetime import datetime, timezone\n",
"from datetime import date\n",
"from astropy import time\n",
"from astropy.coordinates import solar_system_ephemeris\n",
"solar_system_ephemeris.set(\"jpl\")\n",
"from poliastro.ephem import Ephem\n",
"from poliastro.bodies import Sun, Mars, Earth\n",
"from poliastro.twobody import Orbit\n",
"from poliastro.frames import Planes"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Introduction"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We launched a cargo ship with an ion propulsion drive from Earth to Mars. The cargo ship is initially in a circular orbit around the earth at the Earth-Moon $L_1$ point and departs on a parabolic trajectory using a rocket engine. Upon leaving the Earth’s sphere of influence, the ion drive was turned on and remained on until the ship reached the edge of Mars' sphere of influence. The ion drive was disabled, and the cargo ship was captured by Mars. The rocket was then engaged at the periapsis of the hyperbolic trajectory to circularize the orbit. During this trip we determined its orbital radius around Mars, the required mass of propellant for the rocket and the ion engine, and the date the cargo ship left $L_1$. This included the time to transit the spheres of influence of Earth and Mars and the time to transit to Mars. We also compared it to a standard impulsive Hohmann transfer’s transit time and propellant requirement. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Analysis and Results"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Non-Dimensional L1 Point"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$L_1$ is a collinear Lagrange point. To find the collinear Lagrange point, the assumption is made that $y=z=y^*=z^*=0$. In doing so, the possible solutions for $x^*$ and $y^*$ are \n",
"\n",
"$$0 = x^* - \\frac{1 - \\pi_2}{|x^* + \\pi_2|^3}(x^* + \\pi_2) - \\frac{\\pi_2}{|x^* - 1 + \\pi_2|^3}(x^* - 1 + \\pi_2)$$\n",
"$$y^* = 0$$\n",
"\n",
"$L_1$ lies between the Earth and Moon. The cell below defines a function for calculating the values of the collinear Lagrange points."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def collinear_lagrange(xstar, pi_2):\n",
" \"\"\"Calculate the resultant of the collinear Lagrange point equation.\n",
" \n",
" This is a function f(xstar, pi_2), where xstar is the nondimensional x coordinate\n",
" and pi_2 is the nondimensional mass ratio. The function should be passed to\n",
" scipy.optimize.newton (or another Newton solver) to find a value for xstar\n",
" that satsifies the equation, for a given value of pi_2.\n",
" \n",
" The solver will try different values of xstar until the return value is equal to zero.\n",
" \"\"\"\n",
" first_term = xstar\n",
" second_term = (1 - pi_2)/np.abs(xstar + pi_2)**3 * (xstar + pi_2)\n",
" third_term = pi_2 / np.abs(xstar - 1 + pi_2)**3 * (xstar - 1 + pi_2)\n",
" return first_term - second_term - third_term"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In order to design the trajectory, the parameters of the system must be defined. The mass of the Earth, 5.974E24 kg, is defined as $m_{earth}$ and the mass of the Moon, 73.48E21 kg, is defined as $m_{moon}$. The distance between the Earth and Moon, 384,400 km, is defined as $r_{em}$. This distance will be used to non-dimensionalize the equations of motion. $\\pi_2$ is a dimensionless ratio of the masses that is used when calculating the equations of motion. \n",
"\n",
"$$\\pi_2 = \\frac{m_2}{m_1 + m_2}$$"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"m_earth = 5.974E24 # kg\n",
"m_moon = 73.48E21 # kg\n",
"r_em = 384400 # km\n",
"pi_2 = m_moon/(m_earth + m_moon)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Collinear Lagrange function is employed using a newton solver to calculate the value of the non-dimensional Earth-Moon $L_1$ point. However, to put the location of the $L_1$ point into perspective, it can be dimensionalized by multiplying the value by the distance between the two bodies, $r_{em}$."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Earth-Moon L1 point = 321710.307 km\n"
]
}
],
"source": [
"non_dim_L_1 = newton(func=collinear_lagrange, x0=0, args=(pi_2,))\n",
"L_1 = non_dim_L_1*r_em\n",
"print('Earth-Moon L1 point =', round(L_1,3), 'km')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Spheres of Influence"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Close to a planet, that planet's gravity dominates the motion of a satellite. However, at some distance, the Sun's gravity starts to be more important. That distance is known as the sphere of influence (SOI). The SOI is given by \n",
"\n",
"$$r_{SOI} = R*(\\frac{m_p}{m_s})^\\frac{3}{2}$$\n",
"\n",
"where $R$ is the semimajor axis of the planet's orbit, $m_p$ is the mass of the desired planet, and $m_s$ is the mass of the sun. The SOI value was calculated for both Earth and Mars. "
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Earth sphere of influence = 924663.613 km\n",
"Mars sphere of influence = 577134.933 km\n"
]
}
],
"source": [
"m_sun = 1.989E30 # kg\n",
"R_earth = 149.6E6 # km\n",
"r_soi_e = R_earth*((m_earth/m_sun)**(2/5))\n",
"print('Earth sphere of influence =', round(r_soi_e,3),'km')\n",
"\n",
"R_mars = 227.9E6 # km\n",
"m_mars = 641.9E21 # kg\n",
"r_soi_m = R_mars*((m_mars/m_sun)**(2/5))\n",
"print('Mars sphere of influence =', round(r_soi_m,3),'km')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Specific Impulse"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Specific impulse is defined as $$I_{\\text{sp}} = \\frac{\\text{thrust}}{\\text{Sea-level weight rate of fuel consumption}}$$\n",
"\n",
"The type of propellant used in the rocket engine will vary the value of the specific impulse. In the Mars Reconnaissance Orbiter, the spacecraft used monopropellant hydrazine which has an $I_{sp}$ value of 230 seconds. On the spacecraft, there were six large thrusters, six medium thrusters, and eight small thrusters. Each large thruster was capable of producing 170 N of thrust [1]. Even using the large thruster, this propellant would not produce the necessary thrust. \n",
"\n",
"A more practical propellant for the rocket engine might be liquid oxygen/liquid hydrogen. Aerojet Rocketdyne has a spacecraft, RL-10C-1, that produces 101.82 kN of thrust and has a specific impulse value of 450 seconds [2][3]. \n",
"\n",
"The cargo ship is equipped with an ion propulsion drive. The ion propulsion drive produces small amounts of thrust. In this specific case, the ion propulsion drive is capable of producing 5 N. The specific impulse for an ion propulsion drive is > 3000 seconds [4]. For this project, it will be assumed that it is exactly equal to 5000 s. "
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"isp_ipd = 5000 # s\n",
"isp_re = 450 # s"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Circular Orbit"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"When defining the orbit of cargo ship around Earth, we assumed that it is in a circular orbit. Therefore, the location of the $L_1$ point will be equal to the radius, $r_1$ for calculations in this orbit. The gravitational parameter of Earth dominates during this trajectory, so we will employ $mu_earth = 398600$ km/s^2 when necessary.\n",
"\n",
"In a circular orbit, the velocity is defined as\n",
"\n",
"$$v = \\sqrt{\\frac{\\mu}{r}}$$\n",
"\n",
"This value will be necessary to calculate the $\\Delta v$ to leave the circular orbit and go onto the parabolic trajectory."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The velocity of the circular orbit = 1.113 km/s\n"
]
}
],
"source": [
"r_1 = L_1\n",
"mu_earth = 398600\n",
"\n",
"v_1 = np.sqrt(mu_earth/r_1)\n",
"print('The velocity of the circular orbit =', round(v_1,3), 'km/s')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Parabolic Trajectory"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As defined by a parabolic trajectory, $e = 1$. Objects that are very far from the Sun and have nearly zero velocity, but are captured by the Sun and swing around. Therefore, for $e = 1$, $r \\rightarrow \\infty$, and $v \\rightarrow 0$.\n",
"\n",
"Using the values previously defined as $mu_{earth}$ and $r_1$, we can calculate the velocity needed to escape the circular orbit, where\n",
"\n",
"$$v_{esc} = \\sqrt{\\frac{2 \\mu}{r}}$$\n",
"\n",
"To calculate the time to escape the circular trajectory is defined by \n",
"\n",
"$$t_{esc} = \\frac{M_p h^3}{\\mu^2}$$\n",
"\n",
"where $h$ is the magnitude of the angular momentum and $M_p$ is the mean anomaly of the parabolic trajectory.\n",
"\n",
"The magnitude of the angular momentum can be derived from \n",
"\n",
"$$h = r*v$$\n",
"\n",
"and the value of the mean anomaly for a parabolic trajectory is given by\n",
"\n",
"$$M_p = \\frac{1}{2}*tan(\\frac{\\theta}{2}) + \\frac{1}{6}*tan^3(\\frac{\\theta}{2})$$\n",
"\n",
"Theta is the value of the true anomaly, such that\n",
"\n",
"$$\\theta = cos^{-1}(\\frac{h^2}{e \\mu r} - \\frac{1}{e})$$\n",
"\n",
"When calculating the value of $\\theta$, $r$ is equivalent to $r_{SOI}$."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The time it takes to escape Earth's circular orbit = 10.523 days\n",
"The velocity leaving Earth's circular orbit = 1.574 km/s\n"
]
}
],
"source": [
"e = 1\n",
"\n",
"v_esc = np.sqrt((2*mu_earth)/r_1)\n",
"\n",
"h_esc = r_1*v_esc\n",
"\n",
"theta_esc = np.arccos(((h_esc**2)/(e*mu_earth*r_soi_e)) - (1/e))\n",
"\n",
"M_p = ((1/2)*np.tan(theta_esc/2)) + ((1/6)*((np.tan(theta_esc/2))**3))\n",
"\n",
"t_esc = (M_p*(h_esc**3))/(mu_earth**2)\n",
"\n",
"print(\"The time it takes to escape Earth's circular orbit =\", round(t_esc/(3600*24),3),'days')\n",
"print(\"The velocity leaving Earth's circular orbit =\", round(v_esc,3),'km/s')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Impulsive Maneuver to leave Earth's SOI"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The cargo ship utilizies an impulsive maneuver to leave Earth's sphere of influence. The rocket engine supplies the ship with 100 kN of thrust to leave the circular orbit on the parabolic trajectory. The $\\Delta v$ required to make this transition is given by \n",
"\n",
"$$\\Delta v = \\lvert (v_1 - v_{esc}) \\rvert$$\n",
"\n",
"where the values of $v_1$ and $v_{esc}$ were previously calculated. The value of $\\Delta v$ will be used in the calculation of, both, the amount of propellant consumed and the amount of time it takes to leave Earth's SOI.\n",
"\n",
"The amount of propellant consumed during the maneuver comes from the equation\n",
"\n",
"$$\\Delta m = m*(1-e^{\\frac{-\\Delta v}{I_{sp}g_0}})$$\n",
"\n",
"where $m$ is the mass of the cargo ship, $\\Delta v$ is the change in velocity to leave the orbit, $I_{sp}$ is the specific impulse, and $g_0$ is the gravity on Earth.\n",
"\n",
"Furthermore, the amount of time it takes to leave Earth's SOI is given by \n",
"\n",
"$$\\Delta t = \\frac{m \\Delta v}{thrust}$$\n",
"\n",
"where $m$ is the mass of the cargo ship and the thrust is that of the rocket engine."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The amount of propellant consumed during the impulsive maneuver is 2480.064 kg\n",
"The time it takes to leave Earth's sphere of influence is 1.921 mins\n"
]
}
],
"source": [
"thrust_re = 100 # kN\n",
"\n",
"delta_v1 = np.abs(v_1 - v_esc)\n",
"\n",
"mu_sun = 1.3271E11\n",
"\n",
"m_ship = 25000 # kg\n",
"\n",
"g_earth = 9.807E-3 # km/s**2\n",
"\n",
"m_consumed1 = m_ship*(1 - np.exp(-delta_v1/(isp_re*g_earth)))\n",
"print('The amount of propellant consumed during the impulsive maneuver is', round(m_consumed1,3), 'kg')\n",
"\n",
"delta_t1 = (m_ship*delta_v1)/(thrust_re)\n",
"print(\"The time it takes to leave Earth's sphere of influence is\", round(delta_t1/60,3), 'mins')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Non-Impulsive Maneuver to Mars SOI"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Upon leaving Earth's sphere of influence, the ion propulsion drive is turned on and remains on until the ship reaches the edge of Mars' sphere of influence. This is a non-impulsive maneuver.\n",
"\n",
"At the parabolic trajectory, from the perspective of Earth, the cargo ship has a $v_\\infty = 0$. Therefore, geocentrically, the velocity of the ship is equivalent to that of earth."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"v_earth = np.sqrt(mu_sun/R_earth)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The maneuver to reach Mars' SoI is performed by providing nonimpulsive thrust. This means that a very low thrust is provided over a sufficiently long time period to increase efficiency. Since the position is changing during the impulse, the equation of motion will have an additional force term on the right hand side, such that\n",
"\n",
"$$\\ddot{\\vec{r}} = -\\mu\\frac{\\vec{r}}{r^3} + \\frac{\\vec{F}}{m}$$\n",
"\n",
"Assuming that the force is a thrust provided in the same direction as the velocity, \n",
"\n",
"$$\\vec{F} = T\\frac{\\vec{v}}{v}$$\n",
"\n",
"where $T$ is the value of the thrust, $\\vec{v}$ is the velocity, and $v$ is the magnitude of the velocity. The mass of the spacecaft will decrease during this maneuver as the rocket motors eject propellant, where the change in mass is equal to \n",
"\n",
"$$\\frac{dm}{dt} = -\\frac{T}{I_{sp}g_0}$$\n",
"\n",
"The function nonimpulsive_maneuver takes the initial position, velocity, and mass vectors and stores them as an array, $Y_0$. In this case, the initial position is composed of the semimajor axis of Earth's orbit around the Sun in the x-direction. The initial velocity is equal to that of the speed of Earth's orbit around the sun in the y-direction. The mass of the ship is equal to the mass of the ship minus the amount of propellant ejected during the first maneuver. The thrust of the maneuver is supplied as well as the specific impulse, the gravitational parameter, and the sea-level acceleration of gravity."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"r_0 = np.array((R_earth, 0, 0)) # km\n",
"v_0 = np.array((0, v_earth, 0)) # km\n",
"m_0 = np.array((m_ship - m_consumed1)) # kg\n",
"Y_0 = np.hstack((r_0, v_0, m_0))\n",
"\n",
"thrust_ipd = 5.0E-3 # kN\n",
"T = thrust_ipd # kN\n",
"I_sp = isp_ipd # s\n",
"mu = mu_sun\n",
"r_2 = R_mars\n",
"g_0 = g_earth"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"def nonimpulsive_maneuver(t, Y, mu, T, I_sp, g_0, r_2):\n",
" r = np.sqrt(np.dot(Y[0:3], Y[0:3]))\n",
" v = np.sqrt(np.dot(Y[3:6], Y[3:6]))\n",
" m = Y[-1]\n",
" dY_dt = np.zeros(len(Y))\n",
" dY_dt[0:3] = Y[3:6]\n",
" dY_dt[3:6] = -mu * Y[0:3] / r**3 + T * Y[3:6] / (m * v)\n",
" dY_dt[-1] = - T / (I_sp * g_0)\n",
" return dY_dt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The reached_destination function determines whether the spacecraft has reached the specified $r_2$. In this case, $r_2$ is the value of Mars' semimajor axis. This will happen when the orbital radius is equal to $r_2$ or when the difference between the current radius and desired radius is zero. This will cause the integrations to stop.\n",
"\n",
"The function mass returns only the last element of the solution vector, equal to the current value of the mass. If the mass becomes zero before we have reached the final destination, it will cause the integration to terminate."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"def reached_destination(t, Y, mu, T, I_sp, g_0, r_2):\n",
" r_vec = Y[0:3]\n",
" r = np.sqrt(np.dot(r_vec, r_vec))\n",
" return r - r_2\n",
"\n",
"reached_destination.terminal = True"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"def mass(t, Y, mu, T, I_sp, g_0, r_2):\n",
" return Y[-1]\n",
"\n",
"mass.terminal = True"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The time of integration for these functions will be from 0 to the time of one year. The tolerances, rtol and atol, are set to 1E-12 and 1E-15, respectively. These values provide a stable answer and are about as small as possible to set while still impacting the answer. Employing the Runge-Kutta 7th order solver, DOP853, increases the solving speed. The events specify the reached_destination and mass functions from which we will find the time of the transit and the amount of propellant consumed."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"t_0 = 0\n",
"t_final = 3600*24*(365)\n",
"t_eval = np.linspace(t_0, t_final, int(1E6))\n",
"sol = solve_ivp(\n",
" nonimpulsive_maneuver,\n",
" t_span=(t_0, t_final),\n",
" y0=Y_0,\n",
" t_eval=t_eval,\n",
" events=(reached_destination, mass),\n",
" rtol=1E-12,\n",
" atol=1E-15,\n",
" method=\"DOP853\",\n",
" args=(mu, T, I_sp, g_0, r_2)\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The time it takes for the non-impulsive maneuver is 278.59 days\n"
]
}
],
"source": [
"delta_t2 = sol.t_events[0][0]\n",
"print('The time it takes for the non-impulsive maneuver is', round(delta_t2/(3600 * 24),3), 'days')"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The amount of propellant consumed during the non-impulsive maneuver is 2454.382 kg\n"
]
}
],
"source": [
"m = sol.y[-1]\n",
"m_consumed2 = m_0 - m[-1]\n",
"print('The amount of propellant consumed during the non-impulsive maneuver is', round(m_consumed2,3), 'kg')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Periapsis of Hyperbolic Trajectory"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Upon reaching the SOI of Mars, the ion drive is disabled and the cargo ship is captured by Mars. The 100 kN rocket engine will be engaged at periapsis of the hyperbolic trajectory. To find periapsis, we can use the postition and velocity vectors calculated to reach Mars' SOI. \n",
"\n",
"The position and velocity vecotrs are stored in sol.y. The final position and velocity can be found by calling the [-1] index of both vectors arrays. These values are equal to the arrival position and velocity velctors at Mars."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"r_vec = sol.y[0:3].T\n",
"r = np.sqrt(r_vec[:, 0]**2 + r_vec[:, 1]**2 + r_vec[:, 2]**2)\n",
"\n",
"v_vec = sol.y[3:6].T\n",
"v = np.sqrt(v_vec[:, 0]**2 + v_vec[:, 1]**2 + v_vec[:, 2]**2)\n",
"\n",
"r_m_arrival = r_vec[-1]\n",
"v_m_arrival = v_vec[-1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since we know the point on Mars', assumed circular, orbit around the Sun where the cargo ship arrives, we can determine the position and velocity vecotrs of Mars relative to the Sun. The point on Mars' orbit where we arrive is equal to the location of the cargo ship's position. The velocity of Mars is found using the equation for calculating the velocity of a circular orbit, as was done for $v_{earth}$\n",
"\n",
"$$v = \\sqrt{\\frac{\\mu}{r}}$$\n",
"\n",
"The position and velocity vectors of Mars relative to the Sun will be stored in arrays. The position vector is equal to the x and y components at the end of the integration. The velocity vector is composed of $v_{mars}$ in the y-direction."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"r_m_sun = r_m_arrival\n",
"\n",
"r_mag = np.sqrt(r_m_sun[0]**2 + r_m_sun[1]**2 + r_m_sun[2]**2)\n",
"\n",
"v_mars = np.sqrt(mu_sun/r_mag)\n",
"\n",
"v_m_sun = np.array([0,v_mars,0])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Taking the dot product of the velocity of Mars in reference to the Sun and the velocity of the cargo ship when it arrives at Mars will yield the vector magnitude, where\n",
"\n",
"$$v_{mag} = v_{m, arrival_x}v_{m, sun_x} + v_{m, arrival_y}v_{m, sun_y} + v_{m, arrival_z}v_{m, sun_z}$$\n",
"\n",
"Using the vector magnitude, we can calculate the true anomaly of the asymptote, where\n",
"\n",
"$$\\theta_\\infty = 180 - cos(v_{mag})$$\n",
"\n",
"The value of $\\theta_\\infty$ can be used to find the eccentricity of the hyperbolic orbit, where\n",
"\n",
"$$e = \\frac{-1}{cos(\\theta_\\infty)}\n",
"\n",
"To calculate the angular momentum, we also need the value of $v_\\infty$. To find $v_\\infty$, we can take the vector difference between the final velocity of the cargo ship upon arrival at Mars and the velocity of Mars relative to the Sun to calculate the velocity of the cargo ship relative to Mars. The magnitude of the velocity vector is equal to $v_\\infty$. \n",
"\n",
"With the values of $\\mu_{mars}$, $v_\\infty$, and $e$, we can find the magnitude of the angular momentum. The angular momentum is equal to \n",
"\n",
"$$h = \\frac{\\mu}{v_\\infty}\\sqrt{e^2 - 1}$$"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"v_ship_mars = np.abs(v_m_arrival - v_m_sun)\n",
"\n",
"v_dot = np.dot(v_m_arrival,v_m_sun)\n",
"\n",
"theta_inf = 180 - np.degrees(np.cos(v_dot))\n",
"\n",
"e_mars = -1/np.cos(np.radians(theta_inf))\n",
"\n",
"v_ship_mars_mag = np.sqrt(v_ship_mars[0]**2 + v_ship_mars[1]**2 + v_ship_mars[2]**2)\n",
"\n",
"v_m_inf = v_ship_mars_mag\n",
"\n",
"mu_mars = 42828\n",
"\n",
"h_mag = (mu_mars/v_m_inf)*(np.sqrt(e_mars**2 - 1))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using the eccentricity and magnitude of the angular momentum, we can find the periapsis point on the hyperbolic trajectory where\n",
"\n",
"$$r_p = \\frac{h^2}{\\mu}\\frac{1}{1+e}$$\n",
"\n",
"The gravitational parameter that dominates here is that of Mars, equal to 42828 km^3/s^2."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Periapsis of the hyperbolic trajectory = 5.331 km\n"
]
}
],
"source": [
"r_p = ((h_mag**2)/mu_mars)*(1/(1+e_mars))\n",
"print('Periapsis of the hyperbolic trajectory =', round(r_p,3), 'km')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Impulsive Manuever to Circularize the Orbit"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Upon reaching periapsis of the hyperbolic trajectory, the 100 kN rocket engine is implusively engaged to circularize the orbit. To calculate the velocity of the cargo ship as it orbits cirularly around Mars, we can use the equation\n",
"\n",
"$$v = \\sqrt{\\frac{\\mu}{R}}$$\n",
"\n",
"The change in velocity between the hyperbolic trajectory to the circular orbit is needed to calculate the change in time and the amount of propellant consumed, where\n",
"\n",
"$$\\Delta v = \\lvert (v_3 - v_{m, inf}) \\rvert$$\n",
"\n",
"The amount of propellant consumed is given by \n",
"\n",
"$$\\Delta m = m*(1-e^{\\frac{-\\Delta v}{I_{sp}g_0}})$$\n",
"\n",
"where $m$ is the mass of the cargo ship minus the amount of fuel consumed by the first two maneuvers, $\\Delta v$ is the change in velocity upon entering the orbit, $I_{sp}$ is the specific impulse, and $g_0$ is the gravity on Mars.\n",
"\n",
"Furthermore, the amount of time it takes to enter Mars' circular orbit is given by \n",
"\n",
"$$\\Delta t = \\frac{m \\Delta v}{thrust}$$\n",
"\n",
"where $m$ is the mass of the cargo ship minus the amount of fuel consumed by the first two manuevers and the thrust is that of the rocket engine."
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The amount of propellant consumed during the impulsive maneuver is 20065.554 kg\n",
"The time it takes to enter Mars' sphere of influence is 2.41 hrs\n"
]
}
],
"source": [
"g_mars = 3.711E-3 # km/s**2\n",
"\n",
"v_3 = np.sqrt(mu_mars/r_p)\n",
"\n",
"delta_v3 = np.abs(v_3 - v_m_inf)\n",
"\n",
"m_consumed3 = (m_ship - m_consumed1 - m_consumed2)*(1 - np.exp(-delta_v3/(isp_re*g_mars)))\n",
"print('The amount of propellant consumed during the impulsive maneuver is', round(m_consumed3,3), 'kg')\n",
"\n",
"delta_t3 = ((m_ship - m_consumed1 - m_consumed2)*delta_v3)/(thrust_re)\n",
"print(\"The time it takes to enter Mars' sphere of influence is\", round(delta_t3/3600,3), 'hrs')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Total Amount of Propellant Consumed"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The total amount of propellant consumed is the sum of the propellant used during each maneuver\n",
"\n",
"$$m_{total} = m_1 + m_2 + m_3$$"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The total amount of propellant consumed is 25000.0 kg\n"
]
}
],
"source": [
"m_consumed_total = m_consumed1 + m_consumed2 + m_consumed3\n",
"print('The total amount of propellant consumed is', round(m_consumed_total,3),'kg')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Total Transit Time"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The total transit time is the sum of the time taken for each maneuver\n",
"\n",
"$$t_{total} = t_{esc} + t_1 + t_2 + t_3$$"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The total transit time is 289.214 days\n"
]
}
],
"source": [
"t_total = t_esc + delta_t1 + delta_t2 + delta_t3\n",
"print('The total transit time is', round(t_total/(3600*24),3), 'days')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Comparison to Impulsive Hohmann Transfer"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The cargo ship's journey from Earth to Mars can also be modeled as a standard impulsive Hohmann transfer. \n",
"\n",
"The difference in velocity of the cargo ship leaving Earth can be found from\n",
"\n",
"$$\\Delta V_D = \\sqrt{\\frac{mu_{sun}}{R_1}}(\\sqrt{\\frac{2R_2}{R_1 + R_2}}-1)$$\n",
"\n",
"In this equation, $R_1$ is the distance between the Earth-Moon $L_1$ point and the Sun and $R_2$ is the distance between Mars and the Sun. \n",
"\n",
"The corresponding propellant used for this maneuver is given by\n",
"\n",
"$$\\Delta m = m*(1-e^{\\frac{-\\Delta v}{I_{sp}g_0}})$$\n",
"\n",
"where $m$ is the mass of the cargo ship, $\\Delta v$ is the value of $\\Delta V_D$, $I_{sp}$ is the specific impulse of the rocket engine, and $g_0$ is the gravity on Earth."
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"r_2 = 223600000\n",
"r_1 = 147260000 + L_1"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The velocity leaving Earth's SOI = 2.928 km/s\n",
"The amount of propellant it takes to leave Earth's SOI = 12122.834 kg\n"
]
}
],
"source": [
"delta_vd = np.sqrt(mu_sun/r_1)*(np.sqrt((2*r_2)/(r_1+r_2))-1)\n",
"\n",
"delta_m1 = m_ship * (1 - np.exp(-delta_vd/(isp_re*g_earth))) \n",
"\n",
"print(\"The velocity leaving Earth's SOI =\", round(delta_vd,3),'km/s')\n",
"print(\"The amount of propellant it takes to leave Earth's SOI =\", round(delta_m1,3),'kg')"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The velocity to reach Mars = 2.928 km/s\n",
"The amount of propellant it takes to reach Mars = 10646.621 kg\n"
]
}
],
"source": [
"delta_v2 = np.sqrt(mu_sun/r_1) * (math.sqrt((2*r_2)/(r_1+r_2))-1)\n",
"\n",
"delta_m2 = (m_ship - delta_m1) * (1 - np.exp(-delta_v2/(isp_re*g_mars)))\n",
"\n",
"print('The velocity to reach Mars =', round(delta_v2, 3), 'km/s')\n",
"print(\"The amount of propellant it takes to reach Mars =\", round(delta_m2,3),'kg')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"After departing from Earth, the difference in the velocity of the cargo ship entering Mars can be found using the equation\n",
"\n",
"$$\\Delta V_A = \\sqrt{\\frac{mu_{sun}}{R_2}}(1 - \\sqrt{\\frac{2R_1}{R_1 + R_2}})$$\n",
"\n",
"Similar to $\\Delta V_D$, the values of $R_1$ and $R_2$ are the distance between the Earth-Moon $L_1$ point and the Sun and Mars and the Sun, respectively.\n",
"\n",
"With the propellant consumed for this maneuver using the same equation as above,\n",
"\n",
"$$\\Delta m = m*(1-e^{\\frac{-\\Delta v}{I_{sp}g_0}})$$\n",
"\n",
"except the value of $g_0$ is now equal to that of Mars instead of Earth and the $\\Delta v$ is equal to that of $\\Delta V_A$."
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The velocity to circularize Mars = 2.637 km/s\n",
"The amount of propellant it takes to circularize Mars = 1770.817 kg\n"
]
}
],
"source": [
"delta_va = np.sqrt(mu_sun/r_2) * (1 - np.sqrt((2*r_1)/(r_1+r_2)))\n",
"\n",
"delta_m3 = (m_ship - delta_m1 - delta_m2) * (1 - np.exp(-delta_va/(isp_re*g_mars))) \n",
"\n",
"print(\"The velocity to circularize Mars =\", round(delta_va,3),'km/s')\n",
"print(\"The amount of propellant it takes to circularize Mars =\", round(delta_m3,3),'kg')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The transit time for the Hohmann transfer is given by \n",
"\n",
"$$t = \\frac{\\pi}{\\sqrt{\\mu_{sun}}}(\\frac{R_{earth} + R_{mars}}{2})^{\\frac{3}{2}}$$\n",
"\n",
"where $R_{Earth}$ is the semimajor axis of Earth's orbit and $R_{Mars}$ is the semimajor axis of Mars' orbit. The semimajor axis of both orbits is assumed to be the radius of the planets' circular orbits around the Sun.\n",
"\n",
"The transit time represents the time it takes the cargo ship to reach Mars."
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The time it takes to reach Mars = 258.83 days\n"
]
}
],
"source": [
"transit_t = (np.pi/np.sqrt(mu_sun)) * (((R_earth + R_mars)/2)**(3/2))\n",
"print('The time it takes to reach Mars =', round(transit_t/(3600*24), 3), 'days')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The amount of time it takes to leave Earth's SOI and circularize Mars is given by \n",
"\n",
"$$\\Delta t = \\frac{m \\Delta v}{thrust}$$\n",
"\n",
"For the time it takes to leave Earth, $m$ is the mass of the cargo ship, $\\Delta v$ is the change in velocity when departing from Earth, and the thrust is that supplied by the rocket engine.\n",
"\n",
"For the time it takes to circularize Mars, $m$ is the mass of the cargo ship, $\\Delta v$ is the change in velocity when arriving at Mars, and the thrust is that supplied by the rocket engine. "
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The time it takes to leave Earth's SOI = 12.199 mins\n"
]
}
],
"source": [
"delta_t1 = (m_ship * delta_vd)/thrust_re\n",
"print(\"The time it takes to leave Earth's SOI =\", round(delta_t1/60,3), 'mins')"
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The time it takes to circularize Mars = 0.98 mins\n"
]
}
],
"source": [
"delta_t3 = ((m_ship - delta_m1 - delta_m2) * delta_va)/thrust_re\n",
"print(\"The time it takes to circularize Mars =\", round(delta_t3/60,3), 'mins')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The total time of the Hohmann transfer is given by the sum of the time for each part,\n",
"\n",
"$$t = \\Delta t_1 + \\Delta t_3 + t_{transit}$$\n",
"\n",
"The total propellant consumed during the Hohmann transfer is given by the sum of the amount of propellant used for each maneuver,\n",
"\n",
"$$m = \\Delta m_1 + \\Delta m_2 + \\Delta m_3$$"
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The total time = 258.839 days\n",
"The total propellant consumed = 24540.272 kg\n"
]
}
],
"source": [
"total_time = delta_t3 + delta_t1 + transit_t\n",
"\n",
"total_propellant = delta_m1 + delta_m3 + delta_m2\n",
"\n",
"print('The total time =', round(total_time/(3600*24),3), 'days')\n",
"print('The total propellant consumed =', round(total_propellant,3),'kg')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Discussion"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We calculated the Lagrange point to be at 321710.307 km,Earth’s sphere of influence at 924663.613 km and Mar’s sphere of influence at 577134.933 km. Through research we came to the conclusion that the Isp of the rocket would be 5000 s [4]. While the velocity of the circular orbit was 1.113 km/s with the velocity required to leave Earth 1.574 km/s. It took 10.523 days to leave Earth's circular orbit and 1.921 minutes and 2480.064 kg of fuel to leave Earth’s SOI. While the non-impulsive maneuver to reach Mars’ SOI took 2454.382 kg of propellant and 278.59 days. The amount of propellant required to enter Mars is 20065.554 kg and it only takes 2.41 hrs. The total time taken was 289.214 days and 25000 kgof propellant was used.\n",
"\n",
"While using the Hohmann method our velocity to leave Earth and travel to Mars was 2.928 km/s and the velocity to circularize Mars was 2.637 km/s. The propellant needed to depart Earth’s SOI was 12122.834 kg, as well as 10646.621 kg to reach Mars and 1770.817 kg to circularize Mars with a total of 24540.272 kg of propellant needed. While the total time for the trip would be 258.839 days. \n",
"\n",
"The overall difference between the two methods is 459.728 kg of propellant and 19.860 days. So the Hohmann method resulted in an overall decrease of propellant and a significant decrease in flight time. However, the Hohmann transfer used 9642.77 kg of propellant more than the standard method in the initial maneuvers."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Future Calculations"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Planetary Ephemerides"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Determining the state vector of Earth and Mars at a given time is necessary to design the realistic interplanetray mission of the cargo ship. The planetary ephemerides allows the calculation of the phase angle between planets to determine the window of time to reach Mars.\n",
"\n",
"Julius Caesar introduced the Julian calendar consisting of a year with an average length of 365.25 days, slightly longer than the actual solar year. Pope Greory XIII introduced the Gregorian calendar which corrects the average length of the year to 365.2425 days by skipping the leap year under certain conditions. The Julian calendar, however, was adopted before the Gregorian calendar. Therefore a formula is needed to convert a given Gregorian date into a Julian day. The given Gregorian date is composed of three parts: the month number $M$, the day of the month $D$, and the year $Y$. These components are used to define constant that allow for the calculation of the Julian day. The constants are\n",
"\n",
"$$A = INT\\frac{M-14}{12}$$\n",
"$$B = 1461(Y+4800+A)$$\n",
"$$C = 367(M-2-12A)$$\n",
"$$E = INT\\frac{Y+4900+A}{100}$$\n",
"\n",
"with a final Julian date equal to\n",
"\n",
"$$JDN = INT\\frac{B}{4} + INT\\frac{C}{12} - INT\\frac{3E}{4} + D - 32075$$\n",
"\n",
"The function below converts a given Gregorian day, month, and year to the Julian day number."
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
"def gregorian_to_julian_day_number(month, day, year):\n",
" \"\"\"Convert the given proleptic Gregorian date to the equivalent Julian Day Number.\"\"\"\n",
" if month < 1 or month > 12:\n",
" raise ValueError(\"month must be between 1 and 12, inclusive\")\n",
" if day < 1 or day > 31:\n",
" raise ValueError(\"day must be between 1 and 31, inclusive\")\n",
" A = int((month - 14) / 12)\n",
" B = 1461 * (year + 4800 + A)\n",
" C = 367 * (month - 2 - 12 * A)\n",
" E = int((year + 4900 + A) / 100)\n",
" JDN = int(B / 4) + int(C / 12) - int(3 * E / 4) + day - 32075\n",
" return JDN"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Julian date, as described in the function below, is a decimal number that includes the Julian day in the whole number and the fraction of the day towards the next in the decimal."
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [],
"source": [
"def gregorian_to_julian_date(dt):\n",
" \"\"\"Convert a Gregorian date to a Julian Date.\"\"\"\n",
" JDN = gregorian_to_julian_day_number(dt.month, dt.day, dt.year)\n",
" JDT = JDN + (dt.hour - 12) / 24 + dt.minute / 1_440 + dt.second / 86_400 + dt.microsecond / 86_400_000_000\n",
" return JDT"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The functions below relate the given Julian day to the Gregorian date."
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [],
"source": [
"def julian_day_number_to_gregorian(jdn):\n",
" \"\"\"Convert the Julian Day Number to the proleptic Gregorian Month, Day, Year.\"\"\"\n",
" L = jdn + 68569\n",
" N = int(4 * L / 146_097)\n",
" L = L - int((146097 * N + 3) / 4)\n",
" I = int(4000 * (L + 1) / 1_461_001)\n",
" L = L - int(1461 * I / 4) + 31\n",
" J = int(80 * L / 2447)\n",
" day = L - int(2447 * J / 80)\n",
" L = int(J / 11)\n",
" month = J + 2 - 12 * L\n",
" year = 100 * (N - 49) + I + L\n",
" return year, month, day\n",
"\n",
"def julian_date_to_gregorian(jd):\n",
" \"\"\"Convert a decimal Julian Date to the equivalent proleptic Gregorian day and time.\"\"\"\n",
" jdn = int(jd)\n",
" if jdn < 1_721_426:\n",
" raise ValueError(\"Julian Day Numbers less than 1,721,426 are not supported, because\"\n",
" \"Python's date class cannot represent years before AD 1.\")\n",
" year, month, day = julian_day_number_to_gregorian(jdn)\n",
" offset = timedelta(days=(jd % 1), hours=+12)\n",
" dt = datetime(year=year, month=month, day=day, tzinfo=timezone.utc)\n",
" return dt + offset"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can use the Julian dates to calculate the ephemerides for Mars and Earth using the simplified formulas from JPL. The simplified formulas give the Keplerian orbital elements for each planet. The Keplerian elements include the semimajor axis, eccentricity, inclination, true anomaly, argument of periapsis, and longitude of the ascending node. The values produced is given with its absolute vale at a reference point in time, plus the rate of change of that element. The rate of change are given per century, the reference point being J2000 epoch. After calculating the Julian date using the above functions, we can compus the number of Julian centuries that have elapsed, where\n",
"\n",
"$$T = \\frac{JDT - 2451545}{36525}$$\n",
"\n",
"For each planet, Mars and Earth, we must identify the values of the semimajor axis $a$, eccentricity $e$, inclination $i$, mean longitude $L$, longitude of perihelion $\\bar\\omega$, and longitude of the ascending node $\\Omega$. These values are listed on the JPL website in two rows, the first listing the absolute value at J2000 epoch and the second listing the rate of change.\n",
"\n",
"Using the provided values, the true anomaly can be calculated, where\n",
"\n",
"$$M_e = L - \\bar\\omega$$\n",
"\n",
"The mean anomaly allows us to calculate the eccentric anomaly $E$, where\n",
"\n",
"$$M_e = E - esin(E)$$\n",
"\n",
"can be solved using Kepler's equation.\n",
"\n",
"The eccentric anomaly can be related to the true anomaly by\n",
"\n",
"$$\\theta = 2tan^{-1}(\\sqrt{\\frac{1+e}{1-e}}tan(\\frac{E}{2}))$$\n",
"\n",
"The argument of perihelion $\\omega$ can be calculated, where\n",
"\n",
"$$\\omega = \\bar\\omega - \\Omega$$\n",
"\n",
"These parameters define the location and motion of their respective planets, as shown in the cells below for Mars and Earth."
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"𝑎 = 2.2794E+08 km\n",
"𝑒 = 0.09341\n",
"𝑖 = 1.85°\n",
"𝜃 = 81.44°\n",
"𝜔 = -73.35°\n",
"𝛺 = 49.50°\n"
]
}
],
"source": [
"mars = {\n",
" \"a\": 1.52371034, \"a_dot\": 0.00001847,\n",
" \"e\": 0.09339410, \"e_dot\": 0.00007882,\n",
" \"i\": 1.84969142, \"i_dot\": -0.00813131,\n",
" \"L\": -4.55343205, \"L_dot\": 19140.30268499,\n",
" \"long_peri\": -23.94362959, \"long_peri_dot\": 0.44441088,\n",
" \"long_node\": 49.55953891, \"long_node_dot\": -0.29257343,\n",
"}\n",
"T_eph_m = datetime(month=12, day=16, year=2020, hour=21, minute=30, tzinfo=timezone.utc)\n",
"JDT_m = gregorian_to_julian_date(T_eph_m)\n",
"T_m = (JDT_m - 2_451_545) / 36_525\n",
"\n",
"a_m = (mars[\"a\"] + mars[\"a_dot\"] * T_m) * 149_597_870.7\n",
"e_m = mars[\"e\"] + mars[\"e_dot\"] * T_m\n",
"i_m = np.radians(mars[\"i\"] + mars[\"i_dot\"] * T_m)\n",
"L_m = np.radians(mars[\"L\"] + mars[\"L_dot\"] * T_m)\n",
"long_peri_m = np.radians(mars[\"long_peri\"] + mars[\"long_peri_dot\"] * T_m)\n",
"long_node_m = np.radians(mars[\"long_node\"] + mars[\"long_node_dot\"] * T_m)\n",
"\n",
"M_e_m = (L_m - long_peri_m) % (2 * np.pi)\n",
"def kepler(E, M_e, e):\n",
" \"\"\"Kepler's equation, to be used in a Newton solver.\"\"\"\n",
" return E - e * np.sin(E) - M_e\n",
"\n",
"def d_kepler_d_E(E, M_e, e):\n",
" \"\"\"The derivative of Kepler's equation, to be used in a Newton solver.\n",
" \n",
" Note that the argument M_e is unused, but must be present so the function\n",
" arguments are consistent with the kepler function.\n",
" \"\"\"\n",
" return 1 - e * np.cos(E)\n",
"\n",
"E_m = newton(func=kepler, fprime=d_kepler_d_E, x0=np.pi, args=(M_e_m, e_m))\n",
"theta_m = (2 * np.arctan(np.sqrt((1 + e_m) / (1 - e_m)) * np.tan(E_m / 2))) % (2 * np.pi)\n",
"\n",
"omega_m = long_peri_m - long_node_m\n",
"\n",
"print(f\"𝑎 = {a_m:.5G} km\", f\"𝑒 = {e_m:.5F}\", f\"𝑖 = {np.degrees(i_m):.2F}°\", f\"𝜃 = {np.degrees(theta_m):.2F}°\",\n",
" f\"𝜔 = {np.degrees(omega_m):.2F}°\", f\"𝛺 = {np.degrees(long_node_m):.2F}°\", sep=\"\\n\")"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"𝑎 = 1.496E+08 km\n",
"𝑒 = 0.01670\n",
"𝑖 = -0.00°\n",
"𝜃 = 342.11°\n",
"𝜔 = 103.01°\n",
"𝛺 = 0.00°\n"
]
}
],
"source": [
"earth = {\n",
" \"a\": 1.00000261, \"a_dot\": 0.00000562,\n",
" \"e\": 0.01671123, \"e_dot\": -0.00004392,\n",
" \"i\": -0.00001531, \"i_dot\": -0.01294668,\n",
" \"L\": 100.46457166, \"L_dot\": 35999.37244981,\n",
" \"long_peri\": 102.93768193, \"long_peri_dot\": 0.32327364,\n",
" \"long_node\": 0, \"long_node_dot\": 0,\n",
"}\n",
"T_eph_e = datetime(month=12, day=16, year=2020, hour=21, minute=30, tzinfo=timezone.utc)\n",
"JDT_e = gregorian_to_julian_date(T_eph_e)\n",
"T_e = (JDT_e - 2_451_545) / 36_525\n",
"\n",
"a_e = (earth[\"a\"] + earth[\"a_dot\"] * T_e) * 149_597_870.7\n",
"e_e = earth[\"e\"] + earth[\"e_dot\"] * T_e\n",
"i_e = np.radians(earth[\"i\"] + earth[\"i_dot\"] * T_e)\n",
"L_e = np.radians(earth[\"L\"] + earth[\"L_dot\"] * T_e)\n",
"long_peri_e = np.radians(earth[\"long_peri\"] + earth[\"long_peri_dot\"] * T_e)\n",
"long_node_e = np.radians(earth[\"long_node\"] + earth[\"long_node_dot\"] * T_e)\n",
"\n",
"M_e_e = (L_e - long_peri_e) % (2 * np.pi)\n",
"def kepler(E, M_e, e):\n",
" \"\"\"Kepler's equation, to be used in a Newton solver.\"\"\"\n",
" return E - e * np.sin(E) - M_e\n",
"\n",
"def d_kepler_d_E(E, M_e, e):\n",
" \"\"\"The derivative of Kepler's equation, to be used in a Newton solver.\n",
" \n",
" Note that the argument M_e is unused, but must be present so the function\n",
" arguments are consistent with the kepler function.\n",
" \"\"\"\n",
" return 1 - e * np.cos(E)\n",
"\n",
"E_e = newton(func=kepler, fprime=d_kepler_d_E, x0=np.pi, args=(M_e_e, e_e))\n",
"theta_e = (2 * np.arctan(np.sqrt((1 + e_e) / (1 - e_e)) * np.tan(E_e / 2))) % (2 * np.pi)\n",
"\n",
"omega_e = long_peri_e - long_node_e\n",
"\n",
"print(f\"𝑎 = {a_e:.5G} km\", f\"𝑒 = {e_e:.5F}\", f\"𝑖 = {np.degrees(i_e):.2F}°\", f\"𝜃 = {np.degrees(theta_e):.2F}°\",\n",
" f\"𝜔 = {np.degrees(omega_e):.2F}°\", f\"𝛺 = {np.degrees(long_node_e):.2F}°\", sep=\"\\n\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The departure date function, ideally, would take the input value of $\\phi_0$ and solve for the Julian date. The initial phase angle, $\\phi_0$, is the difference between the true anomalies of Earth and Mars (obtained using the simplified formulas for ephemerides). \n",
"\n",
"$$\\phi_0 = \\theta_{mars} - \\theta_{earth}$$\n",
"\n",
"The Julian date can then be converted to its corresponding Gregorian value. \n",
"\n",
"The issue encountered is that the value of $\\phi_0$ is not iterable. In the future, more extensive research would allow us to correctly identify the iterable value of $\\phi_0$ to calculate the departure date."
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [
{
"ename": "TypeError",
"evalue": "departure_date() argument after * must be an iterable, not numpy.float64",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mTypeError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m<ipython-input-37-fc6614719253>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m 10\u001b[0m \u001b[0mphi_0\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtheta_m\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mtheta_e\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 11\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 12\u001b[1;33m \u001b[0mjd\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnewton\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfunc\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mdeparture_date\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mx0\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mJDT_m\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0margs\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mphi_0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 13\u001b[0m \u001b[0mgregorian\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mjulian_date_to_gregorian\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mjd\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m~\\anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py\u001b[0m in \u001b[0;36mnewton\u001b[1;34m(func, x0, fprime, args, tol, maxiter, fprime2, x1, rtol, full_output, disp)\u001b[0m\n\u001b[0;32m 324\u001b[0m \u001b[0mp1\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mx0\u001b[0m \u001b[1;33m*\u001b[0m \u001b[1;33m(\u001b[0m\u001b[1;36m1\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0meps\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 325\u001b[0m \u001b[0mp1\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0meps\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mp1\u001b[0m \u001b[1;33m>=\u001b[0m \u001b[1;36m0\u001b[0m \u001b[1;32melse\u001b[0m \u001b[1;33m-\u001b[0m\u001b[0meps\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 326\u001b[1;33m \u001b[0mq0\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mfunc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mp0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 327\u001b[0m \u001b[0mfuncalls\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 328\u001b[0m \u001b[0mq1\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mfunc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mp1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;31mTypeError\u001b[0m: departure_date() argument after * must be an iterable, not numpy.float64"
]
}
],
"source": [
"def departure_date(jd, phi_0):\n",
" epoch = Time(jd, format=\"jd\")\n",
" mars = Ephem.from_body(Mars, epoch.tdb, attractor=Sun, plane=Planes.EARTH_ECLIPTIC)\n",
" ss_mars = Orbit.from_ephem(Sun, Mars, epoch)\n",
" earth = Ephem.from_body(Earth, epoch.tdb, attractor=Sun, plane=Planes.EARTH_ECLIPTIC)\n",
" ss_earth = Orbit.from_ephem(Sun, Earth, epoch)\n",
" return ss_mars.nu.value - ss_earth.nu.value - phi_0\n",
"\n",
"\n",
"phi_0 = theta_m - theta_e\n",
"\n",
"jd = newton(func=departure_date, x0=JDT_m, args=(phi_0))\n",
"gregorian = julian_date_to_gregorian(jd)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Conclusions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The total time taken using our first method was 278.699 daysand 25000 kg of propellant was used. While using the Hohmann method a total of 24540.272 kg of propellant was needed and the total time for the trip would be 258.839 days. The overall difference between the two methods is 459.728 kg of propellant and 19.860 days. So the Hohmann method resulted in an overall decrease of propellant and a significant decrease in flight time. Yet, the Hohmann transfer requires a mcuh larger initial use of propellant, requiring 9642.77 kg of propellant more than the standard method. The tradeoffs in our design choice are that having a high isp uses more fuel. Some alternatives we could do to improve our design would be to use more engines to make circularizing Mars easier. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## References"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[1]“RL10 Engine,” RL10 Engine | Aerojet Rocketdyne, 2020. [Online]. Available: https://www.rocket.com/space/liquid-engines/rl10-engine. [Accessed: 14-Dec-2020].\n",
"\n",
"[2]“Comparison of orbital rocket engines,” Wikipedia, 08-Dec-2020. [Online]. Available: https://en.wikipedia.org/wiki/Comparison_of_orbital_rocket_engines. [Accessed: 14-Dec-2020].\n",
"\n",
"[3]H. D. Curtis, Orbital mechanics for engineering students: revised fourth edition. Amsterdam: Elsevier, 2021.\n",
"\n",
"[4]“Propulsion,” NASA. [Online]. Available: https://mars.nasa.gov/mro/mission/spacecraft/parts/propulsion/. [Accessed: 14-Dec-2020]. "
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}