diff --git a/GEO_orbital_simulation.ipynb b/GEO_orbital_simulation.ipynb
new file mode 100644
index 0000000..86a7cb5
--- /dev/null
+++ b/GEO_orbital_simulation.ipynb
@@ -0,0 +1,669 @@
+{
+ "nbformat": 4,
+ "nbformat_minor": 0,
+ "metadata": {
+ "colab": {
+ "provenance": [],
+ "authorship_tag": "ABX9TyNehg1/9TgLhaE6psc6cCCX",
+ "include_colab_link": true
+ },
+ "kernelspec": {
+ "name": "python3",
+ "display_name": "Python 3"
+ },
+ "language_info": {
+ "name": "python"
+ }
+ },
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "view-in-github",
+ "colab_type": "text"
+ },
+ "source": [
+ "![\"Open](\"https://colab.research.google.com/assets/colab-badge.svg\")
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "# Intro & Goal of this program\n",
+ "\n",
+ "This program simulates the orbital profile of a payload carried by a two-stage rocket from a LEO orbit to insertion in a GEO orbit. To achieve this, the spacecraft performs a classic Hohmann transfer maneuver. The code implements realistic perturbing forces, such as the gravitational influence of the Moon, the Sun, J-type forces, etc. The goal of this program is to determine what is the associated delta-V cost and fuel mass needed to insert a payload in a GEO orbit, given specific launch parameters such as launch latitude, initial LEO altitude, etc."
+ ],
+ "metadata": {
+ "id": "EfZO5Rd4jXec"
+ }
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {
+ "colab": {
+ "base_uri": "https://localhost:8080/"
+ },
+ "id": "BjpU7NeZhLM-",
+ "outputId": "5bdfe54d-e654-45a9-d2f2-987df4745667"
+ },
+ "outputs": [
+ {
+ "output_type": "stream",
+ "name": "stdout",
+ "text": [
+ "Requirement already satisfied: hapsira in /usr/local/lib/python3.12/dist-packages (0.18.0)\n",
+ "Requirement already satisfied: astropy<7.0 in /usr/local/lib/python3.12/dist-packages (6.1.7)\n",
+ "Requirement already satisfied: astroquery in /usr/local/lib/python3.12/dist-packages (from hapsira) (0.4.11)\n",
+ "Requirement already satisfied: jplephem in /usr/local/lib/python3.12/dist-packages (from hapsira) (2.24)\n",
+ "Requirement already satisfied: matplotlib<3.8 in /usr/local/lib/python3.12/dist-packages (from hapsira) (3.7.5)\n",
+ "Requirement already satisfied: numba in /usr/local/lib/python3.12/dist-packages (from hapsira) (0.60.0)\n",
+ "Requirement already satisfied: numpy in /usr/local/lib/python3.12/dist-packages (from hapsira) (1.26.4)\n",
+ "Requirement already satisfied: pandas in /usr/local/lib/python3.12/dist-packages (from hapsira) (2.2.2)\n",
+ "Requirement already satisfied: plotly in /usr/local/lib/python3.12/dist-packages (from hapsira) (5.24.1)\n",
+ "Requirement already satisfied: pyerfa in /usr/local/lib/python3.12/dist-packages (from hapsira) (2.0.1.5)\n",
+ "Requirement already satisfied: scipy in /usr/local/lib/python3.12/dist-packages (from hapsira) (1.16.3)\n",
+ "Requirement already satisfied: astropy-iers-data>=0.2024.10.28.0.34.7 in /usr/local/lib/python3.12/dist-packages (from astropy<7.0) (0.2026.3.30.0.54.34)\n",
+ "Requirement already satisfied: PyYAML>=3.13 in /usr/local/lib/python3.12/dist-packages (from astropy<7.0) (6.0.3)\n",
+ "Requirement already satisfied: packaging>=19.0 in /usr/local/lib/python3.12/dist-packages (from astropy<7.0) (26.0)\n",
+ "Requirement already satisfied: contourpy>=1.0.1 in /usr/local/lib/python3.12/dist-packages (from matplotlib<3.8->hapsira) (1.3.3)\n",
+ "Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.12/dist-packages (from matplotlib<3.8->hapsira) (0.12.1)\n",
+ "Requirement already satisfied: fonttools>=4.22.0 in /usr/local/lib/python3.12/dist-packages (from matplotlib<3.8->hapsira) (4.62.1)\n",
+ "Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.12/dist-packages (from matplotlib<3.8->hapsira) (1.5.0)\n",
+ "Requirement already satisfied: pillow>=6.2.0 in /usr/local/lib/python3.12/dist-packages (from matplotlib<3.8->hapsira) (11.3.0)\n",
+ "Requirement already satisfied: pyparsing>=2.3.1 in /usr/local/lib/python3.12/dist-packages (from matplotlib<3.8->hapsira) (3.3.2)\n",
+ "Requirement already satisfied: python-dateutil>=2.7 in /usr/local/lib/python3.12/dist-packages (from matplotlib<3.8->hapsira) (2.9.0.post0)\n",
+ "Requirement already satisfied: requests>=2.19 in /usr/local/lib/python3.12/dist-packages (from astroquery->hapsira) (2.32.4)\n",
+ "Requirement already satisfied: beautifulsoup4>=4.8 in /usr/local/lib/python3.12/dist-packages (from astroquery->hapsira) (4.13.5)\n",
+ "Requirement already satisfied: html5lib>=0.999 in /usr/local/lib/python3.12/dist-packages (from astroquery->hapsira) (1.1)\n",
+ "Requirement already satisfied: keyring>=15.0 in /usr/local/lib/python3.12/dist-packages (from astroquery->hapsira) (25.7.0)\n",
+ "Requirement already satisfied: pyvo>=1.5 in /usr/local/lib/python3.12/dist-packages (from astroquery->hapsira) (1.8.1)\n",
+ "Requirement already satisfied: llvmlite<0.44,>=0.43.0dev0 in /usr/local/lib/python3.12/dist-packages (from numba->hapsira) (0.43.0)\n",
+ "Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.12/dist-packages (from pandas->hapsira) (2025.2)\n",
+ "Requirement already satisfied: tzdata>=2022.7 in /usr/local/lib/python3.12/dist-packages (from pandas->hapsira) (2025.3)\n",
+ "Requirement already satisfied: tenacity>=6.2.0 in /usr/local/lib/python3.12/dist-packages (from plotly->hapsira) (9.1.4)\n",
+ "Requirement already satisfied: soupsieve>1.2 in /usr/local/lib/python3.12/dist-packages (from beautifulsoup4>=4.8->astroquery->hapsira) (2.8.3)\n",
+ "Requirement already satisfied: typing-extensions>=4.0.0 in /usr/local/lib/python3.12/dist-packages (from beautifulsoup4>=4.8->astroquery->hapsira) (4.15.0)\n",
+ "Requirement already satisfied: six>=1.9 in /usr/local/lib/python3.12/dist-packages (from html5lib>=0.999->astroquery->hapsira) (1.17.0)\n",
+ "Requirement already satisfied: webencodings in /usr/local/lib/python3.12/dist-packages (from html5lib>=0.999->astroquery->hapsira) (0.5.1)\n",
+ "Requirement already satisfied: SecretStorage>=3.2 in /usr/local/lib/python3.12/dist-packages (from keyring>=15.0->astroquery->hapsira) (3.5.0)\n",
+ "Requirement already satisfied: jeepney>=0.4.2 in /usr/local/lib/python3.12/dist-packages (from keyring>=15.0->astroquery->hapsira) (0.9.0)\n",
+ "Requirement already satisfied: jaraco.classes in /usr/local/lib/python3.12/dist-packages (from keyring>=15.0->astroquery->hapsira) (3.4.0)\n",
+ "Requirement already satisfied: jaraco.functools in /usr/local/lib/python3.12/dist-packages (from keyring>=15.0->astroquery->hapsira) (4.4.0)\n",
+ "Requirement already satisfied: jaraco.context in /usr/local/lib/python3.12/dist-packages (from keyring>=15.0->astroquery->hapsira) (6.1.2)\n",
+ "Requirement already satisfied: charset_normalizer<4,>=2 in /usr/local/lib/python3.12/dist-packages (from requests>=2.19->astroquery->hapsira) (3.4.6)\n",
+ "Requirement already satisfied: idna<4,>=2.5 in /usr/local/lib/python3.12/dist-packages (from requests>=2.19->astroquery->hapsira) (3.11)\n",
+ "Requirement already satisfied: urllib3<3,>=1.21.1 in /usr/local/lib/python3.12/dist-packages (from requests>=2.19->astroquery->hapsira) (2.5.0)\n",
+ "Requirement already satisfied: certifi>=2017.4.17 in /usr/local/lib/python3.12/dist-packages (from requests>=2.19->astroquery->hapsira) (2026.2.25)\n",
+ "Requirement already satisfied: cryptography>=2.0 in /usr/local/lib/python3.12/dist-packages (from SecretStorage>=3.2->keyring>=15.0->astroquery->hapsira) (43.0.3)\n",
+ "Requirement already satisfied: more-itertools in /usr/local/lib/python3.12/dist-packages (from jaraco.classes->keyring>=15.0->astroquery->hapsira) (10.8.0)\n",
+ "Requirement already satisfied: cffi>=1.12 in /usr/local/lib/python3.12/dist-packages (from cryptography>=2.0->SecretStorage>=3.2->keyring>=15.0->astroquery->hapsira) (2.0.0)\n",
+ "Requirement already satisfied: pycparser in /usr/local/lib/python3.12/dist-packages (from cffi>=1.12->cryptography>=2.0->SecretStorage>=3.2->keyring>=15.0->astroquery->hapsira) (3.0)\n"
+ ]
+ }
+ ],
+ "source": [
+ "#Install necessary packages for Google Colab. Note: hapsira is the modern version of poliastro to handle Python 3.12\n",
+ "!pip install hapsira \"astropy<7.0\""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "# Code Imports\n",
+ "This code utilizes several Python libraries for different functionalities:\n",
+ "\n",
+ "1. poliastro: defines celestial bodies and physical constants, calculates the evolution of the orbit using Cowell Propagator, and introduces the effects of solar radiation and other forces\n",
+ "\n",
+ "2. astropy: handling time objects, astronomical units, etc."
+ ],
+ "metadata": {
+ "id": "_sNzUodbmjgF"
+ }
+ },
+ {
+ "cell_type": "code",
+ "source": [
+ "from astropy import units as u\n",
+ "from hapsira.twobody import Orbit\n",
+ "from hapsira.core.perturbations import (\n",
+ " J2_perturbation, J3_perturbation, third_body,\n",
+ " radiation_pressure, atmospheric_drag_exponential\n",
+ ")\n",
+ "from hapsira.ephem import build_ephem_interpolant, Ephem\n",
+ "from hapsira.bodies import Earth, Moon, Sun\n",
+ "from hapsira.maneuver import Maneuver\n",
+ "from hapsira.constants import rho0_earth, H0_earth\n",
+ "\n",
+ "import numpy as np\n",
+ "from astropy.time import Time\n",
+ "from astropy.coordinates import CartesianRepresentation, CartesianDifferential, SkyCoord, ICRS, GCRS\n",
+ "from astropy.coordinates import solar_system_ephemeris\n",
+ "from hapsira.twobody.propagation import CowellPropagator\n",
+ "from hapsira.core.propagation import func_twobody"
+ ],
+ "metadata": {
+ "id": "pChUe3tU9j3G"
+ },
+ "execution_count": 15,
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "# Interpolating values\n",
+ "\n",
+ "Because this program is meant to run multiple times with different launch parameters, we provide a list of apporixmate delta-V losses during takeoff to reach a specific LEO altitude. Interpolation is then used if any delta-V values is needed that is in between any of the provided LEO altitudes. This saves us the trouble of having to make complex calculations for the delta-V losses due to gravity and atmospheric drag for a specific altitude, although more sophisticated models can include this to increase accuracy"
+ ],
+ "metadata": {
+ "id": "_rNn-mNi92Dz"
+ }
+ },
+ {
+ "cell_type": "code",
+ "source": [
+ "\n",
+ "# ===============================================================\n",
+ "# Simulation and Interpolation Configuration\n",
+ "# ===============================================================\n",
+ "# Example interpolation data for gravitational and atmospheric losses based on LEO altitudes\n",
+ "LEO_ALTITUDES_KM = [200, 300, 400, 500, 600, 700, 800]\n",
+ "DELTA_V_GRAV_KMS = [1520, 1540, 1560, 1580, 1600, 1620, 1640] # in m/s (example values)\n",
+ "DELTA_V_ATM_KMS = [320, 280, 250, 230, 210, 190, 170] # in m/s (example values)\n"
+ ],
+ "metadata": {
+ "id": "eOVI25ac_Guv"
+ },
+ "execution_count": 16,
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "# Physical Constants and Ephemeris Setup\n",
+ "\n",
+ "Because we will be computing drag, solar radiation, and other types of forces, we define a few constants with the necessary units for future calculations. We also set a starting epoch and initialize two vectors containing the positions of the Sun and the Moon from the starting time until one day later for every half an hour. This is more computationally efficient than calculating the exact position of the celestial bodies for every time step of the simulation."
+ ],
+ "metadata": {
+ "id": "M4UynXSc_KbF"
+ }
+ },
+ {
+ "cell_type": "code",
+ "source": [
+ "# ===============================================================\n",
+ "# Global Constants and Ephemeris Setup\n",
+ "# ===============================================================\n",
+ "# Earth-related constants\n",
+ "EARTH_RADIUS_KM = Earth.R.to(u.km).value\n",
+ "EARTH_K = Earth.k.to(u.km**3 / u.s**2).value\n",
+ "DRAG_COEFFICIENT = 2.2 # dimensionless\n",
+ "\n",
+ "# Atmospheric model constants (Earth)\n",
+ "SCALE_HEIGHT_KM = H0_earth.to(u.km).value\n",
+ "SEA_LEVEL_DENSITY = rho0_earth.to(u.kg / u.km**3).value\n",
+ "GLOBAL_AREA_TO_MASS = 0 # Placeholder value; updated later in calculation\n",
+ "\n",
+ "#Solar radiation constants\n",
+ "rad_coeff = 1.3\n",
+ "Wdivc_s = 3.828e20/299792.458\n",
+ "\n",
+ "# Define time epoch and set solar system ephemerides\n",
+ "EPOCH = Time(2451543.597, format=\"jd\", scale=\"tdb\")\n",
+ "solar_system_ephemeris.set(\"de432s\")\n",
+ "\n",
+ "# Build ephemeris interpolants for Moon and Sun (half-hour resolution over one day)\n",
+ "start_time = EPOCH\n",
+ "end_time = EPOCH + 1 * u.day\n",
+ "time_step = 0.5 * u.hr\n",
+ "\n",
+ "# Calculate the number of points needed to include both start and end times at the given step.\n",
+ "# For 1 day (24h) and 0.5h step: (24/0.5) + 1 = 48 + 1 = 49 points.\n",
+ "num_points = int(round((end_time - start_time).to(u.hr).value / time_step.to(u.hr).value)) + 1\n",
+ "times = start_time + np.arange(num_points) * time_step\n",
+ "\n",
+ "moon_ephem = build_ephem_interpolant(\n",
+ " Moon,\n",
+ " times, # Changed to positional argument\n",
+ " attractor=Earth # Explicitly set attractor\n",
+ ")\n",
+ "sun_ephem = build_ephem_interpolant(\n",
+ " Sun,\n",
+ " times, # Changed to positional argument\n",
+ " attractor=Earth # Explicitly set attractor\n",
+ ")\n"
+ ],
+ "metadata": {
+ "id": "KMyHUYvR_Ktu"
+ },
+ "execution_count": 17,
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "# Initializing the LEO orbit\n",
+ "\n",
+ "Here, we initialize the Orbit object that defines our starting LEO orbit. The altitude and the inclination are user-defined, while all other parameters are set to 0."
+ ],
+ "metadata": {
+ "id": "vJHHI11HA-Hm"
+ }
+ },
+ {
+ "cell_type": "code",
+ "source": [
+ "def create_leo_orbit(leo_altitude_km, inclination_deg):\n",
+ " \"\"\"\n",
+ " Create a circular LEO orbit given altitude and inclination.\n",
+ "\n",
+ " Parameters:\n",
+ " leo_altitude_km (float): Altitude above the Earth’s surface in km.\n",
+ " inclination_deg (float): Orbit inclination in degrees.\n",
+ "\n",
+ " Returns:\n",
+ " Orbit: A poliastro Orbit object created from classical orbital elements.\n",
+ " \"\"\"\n",
+ " semi_major_axis = (leo_altitude_km + 6378) * u.km\n",
+ " eccentricity = 0.0 * u.one\n",
+ " inclination = inclination_deg * u.deg\n",
+ " raan = 0 * u.deg\n",
+ " arg_periapsis = 0 * u.deg\n",
+ " true_anomaly = 0 * u.deg\n",
+ "\n",
+ " return Orbit.from_classical(Earth, semi_major_axis, eccentricity,\n",
+ " inclination, raan, arg_periapsis, true_anomaly)"
+ ],
+ "metadata": {
+ "id": "HDeyr9oXA-wS"
+ },
+ "execution_count": 18,
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "# Computing acceleration vector\n",
+ "\n",
+ "This function simply computes the accelerations caused by each force on the spacecraft at a specific time. To do this, we use poliastro's built-in functions for each force, the physical constants we defined earlier, and the arrays containing the position of the celestial bodies, as appropiate. It then adds up all the accelerations to get a final acceleration vector and returns this value."
+ ],
+ "metadata": {
+ "id": "ntzKSHj5BY_c"
+ }
+ },
+ {
+ "cell_type": "code",
+ "source": [
+ "# ===============================================================\n",
+ "# Perturbation and Propagation Functions\n",
+ "# ===============================================================\n",
+ "def perturbing_acceleration(t0, state, gravitational_parameter, J2, J3, radius,\n",
+ " drag_coeff, area_to_mass, scale_height, density0, Wdivc_s, rad_coeff, third_body_ephem):\n",
+ " \"\"\"\n",
+ " Compute total perturbing acceleration including J2, J3, atmospheric drag,\n",
+ " and third-body gravitational effects.\n",
+ " \"\"\"\n",
+ " accel_J2 = J2_perturbation(t0, state, gravitational_parameter, J2, radius)\n",
+ " accel_drag = atmospheric_drag_exponential(t0, state, gravitational_parameter, radius,\n",
+ " drag_coeff, area_to_mass, scale_height, density0)\n",
+ " accel_moon = third_body(t0, state, gravitational_parameter,\n",
+ " k_third=Moon.k.to(u.km**3 / u.s**2).value,\n",
+ " perturbation_body=moon_ephem)\n",
+ " accel_sun = third_body(t0, state, gravitational_parameter,\n",
+ " k_third=Sun.k.to(u.km**3 / u.s**2).value,\n",
+ " perturbation_body=third_body_ephem)\n",
+ " accel_J3 = J3_perturbation(t0, state, gravitational_parameter, J3, radius)\n",
+ " accel_rad = radiation_pressure(t0, state, gravitational_parameter, radius, rad_coeff, area_to_mass, Wdivc_s, third_body_ephem)\n",
+ "\n",
+ " return accel_J2 + accel_drag + accel_moon + accel_sun + accel_J3"
+ ],
+ "metadata": {
+ "id": "fbSfX8CKBZUV"
+ },
+ "execution_count": 19,
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "# Computing the derivative of the state vector\n",
+ "\n",
+ "This function computes the 2-body Keplerian effect on the spacecraft and adds the total acceleration vector from the other forces to obtain the derivative of the state vector at a specific time. This derivative is a vector containing the velocity and acceleration of the spacecraft in the x, y and z directions. It is used as an indicator of how the position and velocity of the object is chnaging at that point in time."
+ ],
+ "metadata": {
+ "id": "KOiVl6oODUso"
+ }
+ },
+ {
+ "cell_type": "code",
+ "source": [
+ "def total_dynamics(t0, state, gravitational_parameter):\n",
+ " \"\"\"\n",
+ " Compute the total time derivative of the state vector (Keplerian dynamics\n",
+ " plus perturbations).\n",
+ " \"\"\"\n",
+ " du_kep = func_twobody(t0, state, gravitational_parameter)\n",
+ " ax, ay, az = perturbing_acceleration(\n",
+ " t0, state, gravitational_parameter,\n",
+ " Earth.J2.value, Earth.J3.value,\n",
+ " EARTH_RADIUS_KM, DRAG_COEFFICIENT, GLOBAL_AREA_TO_MASS,\n",
+ " SCALE_HEIGHT_KM, SEA_LEVEL_DENSITY, Wdivc_s, rad_coeff, sun_ephem\n",
+ " )\n",
+ " perturbation = np.array([0, 0, 0, ax, ay, az])\n",
+ " return du_kep + perturbation"
+ ],
+ "metadata": {
+ "id": "r7ISbywRDVLt"
+ },
+ "execution_count": 20,
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "# Optimizing the rocket's mass\n",
+ "\n",
+ "In this simulation, it is assumed that the two-stage rocket provides all the delta-V needed to go from Earth to a GTO orbit, where the payload leaves and circularizes into a GEO orbit. This function looks at the most efficient way to split this delta-V between the two stages of the rocket to minimize its mass.\n",
+ "\n",
+ "Based on typical applications and previous literature, we assume the delta-V provided by stage 1 is between 25% to 40% of the total. We test each of the 16 possible scenarios (e.g. going 1% up each time from 25% to 40%) to see which configuration is the most efficient. To do this, we calculate the mass of each stage and add them all up to get the total initial mass. The mass of each stage is calculated using the following equation:\n",
+ "The modified Tsiolkovsky rocket equation accounting for structural mass fraction ($\\epsilon$) is:\n",
+ "\n",
+ "$$M_{stage} = \\frac{M_{stages} \\left( e^{\\frac{\\Delta V}{V_e}} - 1 \\right)}{1 - \\epsilon e^{\\frac{\\Delta V}{V_e}}}$$\n",
+ "\n",
+ "This is a modified version of Tsiolkovsky rocket equation that accounts for structural mass fraction $\\epsilon$. $V_e$ is the exhaust velocity of such stage. Δ $V$ is the delta-V for that stage. ${M_{stages}}$ represents sum of the masses of each stage that comes after the current stage.\n",
+ "\n",
+ "The derivation of this equation can be found on this website: https://space.geometrian.com/calcs/opt-multi-stage.php"
+ ],
+ "metadata": {
+ "id": "YV1wTXoPFAr_"
+ }
+ },
+ {
+ "cell_type": "code",
+ "source": [
+ "def optimize_rocket_mass(exhaust_vel1, exhaust_vel2, satellite_mass, mag_deltaV, deltaV_takeoff):\n",
+ " \"\"\"\n",
+ " Optimize the distribution of delta-V between stage 1 and stage 2 to minimize\n",
+ " the total rocket mass. The function computes stage masses over a range of splits.\n",
+ "\n",
+ " Parameters:\n",
+ " exhaust_vel1 (float): Exhaust velocity for the first stage.\n",
+ " exhaust_vel2 (float): Exhaust velocity for the second stage.\n",
+ " satellite_mass (float): Dry mass of the satellite.\n",
+ " mag_deltaV (float): Magnitude of the computed delta-V from orbit propagation (km/s).\n",
+ " deltaV_takeoff (float): Delta-V required for take-off (km/s).\n",
+ "\n",
+ " Returns:\n",
+ " tuple: Contains stage1_mass_loss, stage2_mass_loss, wet_stage2_mass,\n",
+ " deltaV_stage1, deltaV_stage2 corresponding to the optimal split.\n",
+ " \"\"\"\n",
+ " total_required_dv = (mag_deltaV * 1000.0) + deltaV_takeoff\n",
+ " start_fraction = 0.25\n",
+ "\n",
+ " # Generate candidate splits for stage 1 delta-V allocation (16 candidates)\n",
+ " deltaV_stage1_candidates = [start_fraction * total_required_dv + i * 0.01 * total_required_dv for i in range(16)]\n",
+ " deltaV_stage2_candidates = [total_required_dv - dv1 for dv1 in deltaV_stage1_candidates]\n",
+ "\n",
+ " stage1_mass_losses = [None] * 16\n",
+ " stage2_mass_losses = [None] * 16\n",
+ " total_rocket_masses = [None] * 16\n",
+ "\n",
+ " # Compute mass requirements for each candidate split using the rocket equation with\n",
+ " # additional mass fractions (0.06 for stage 1, 0.04 for stage 2)\n",
+ " for i in range(16):\n",
+ " stage2_mass_losses[i] = (satellite_mass *\n",
+ " (np.e**(deltaV_stage2_candidates[i] / exhaust_vel2) - 1.0)\n",
+ " / (1.0 - (np.e**(deltaV_stage2_candidates[i] / exhaust_vel2) * 0.04))\n",
+ " )\n",
+ " stage1_mass_losses[i] = ((satellite_mass + stage2_mass_losses[i]) *\n",
+ " (np.e**(deltaV_stage1_candidates[i] / exhaust_vel1) - 1.0)\n",
+ " / (1.0 - (np.e**(deltaV_stage1_candidates[i] / exhaust_vel1) * 0.06))\n",
+ " )\n",
+ " total_rocket_masses[i] = stage1_mass_losses[i] + stage2_mass_losses[i] + satellite_mass\n",
+ "\n",
+ " # Find candidate index with the minimum total rocket mass\n",
+ " optimal_index = np.argmin(total_rocket_masses)\n",
+ "\n",
+ " return (\n",
+ " stage1_mass_losses[optimal_index] * 0.06,\n",
+ " stage2_mass_losses[optimal_index] * 0.04,\n",
+ " stage2_mass_losses[optimal_index],\n",
+ " deltaV_stage1_candidates[optimal_index],\n",
+ " deltaV_stage2_candidates[optimal_index]\n",
+ " )\n"
+ ],
+ "metadata": {
+ "id": "9Y1lW1ZPFBEU"
+ },
+ "execution_count": 21,
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "# Propagating the orbit anc calculating delta-V and fuel mass\n",
+ "\n",
+ "This function performs all the necessary maneuvers to insert the satellite into GEO at the desired inclination. It begins by computing the velocity of the spacecraft at LEO based on its altitude using the Vis-Viva equation:\n",
+ "\n",
+ "$$v = \\sqrt{\\mu \\left( \\frac{2}{r} - \\frac{1}{a} \\right)}$$\n",
+ "\n",
+ "where $\\mu$ is Newton's gravitational constant times mass of the Earth, $r$ is the distance from the earth, and $a$ is the semi-major axis of the orbit (1 for circular orbits).\n",
+ "\n",
+ "This same equation is used to compute the expected velocity of the spacecraft in a GTO orbit at the same distance $r$ from the Earth. The delta-V for insertion into GTO us simply the difference between these two velocities. Trigonometry is used to apply the burn in a direction where the inclination is unchanged. The orbit is then propagated by half of the period of the current GTO orbit\n",
+ "\n",
+ "Next, the code performs a combined plane change maneuver where both the inclination and the altitude of the orbit are changed to insert the payload into GEO with a desired inclination. We use the law of cosines to find the magnitude of this delta-V and the law of sines to find how this delta-V should be split between the x, y and z components. The maneuver is performed and the fuel mass for this insertion is calculated using Tsiolkovsky rocket equation.\n",
+ "\n"
+ ],
+ "metadata": {
+ "id": "mJC9RKZgNbfj"
+ }
+ },
+ {
+ "cell_type": "code",
+ "source": [
+ "def propagate_and_compute_transfer_dv(leo_orbit, exhaust_vel1, exhaust_vel2, exhaust_vel_prop,\n",
+ " satellite_mass, deltaV_takeoff, inc_change_percent):\n",
+ " \"\"\"\n",
+ " Propagate the orbit from LEO to a transfer orbit and compute the transfer delta-V\n",
+ " components along with required propellant loss.\n",
+ "\n",
+ " Parameters:\n",
+ " leo_orbit (Orbit): Initial LEO orbit.\n",
+ " exhaust_vel1 (float): Exhaust velocity for the first stage.\n",
+ " exhaust_vel2 (float): Exhaust velocity for the second stage.\n",
+ " exhaust_vel_prop (float): Propulsive exhaust velocity.\n",
+ " satellite_mass (float): Satellite mass.\n",
+ " deltaV_takeoff (float): Delta-V required for take-off.\n",
+ " inc_change_percent (float): Percentage factor for inclination change.\n",
+ "\n",
+ " Returns:\n",
+ " tuple: (magnitude_deltaV, magnitude_plane_change_deltaV, total_propellant_loss)\n",
+ " \"\"\"\n",
+ " # Compute LEO orbital parameters\n",
+ " leo_semi_major_km = leo_orbit.a.to(u.km).value\n",
+ " leo_velocity_kms = np.sqrt((Earth.k.value / 1000.0) / leo_semi_major_km) / 1000.0\n",
+ "\n",
+ " # Define GTO parameters from LEO altitude\n",
+ " gto_semi_major_km = (42164.0 + leo_semi_major_km) / 2.0\n",
+ " gto_perigee_velocity = np.sqrt((Earth.k.value/1000.0) * ((2.0/leo_semi_major_km) - (1.0/gto_semi_major_km))) / 1000.0\n",
+ " transfer_deltaV = gto_perigee_velocity - leo_velocity_kms\n",
+ "\n",
+ " # Delta-V components (projected along orbital plane components)\n",
+ " dv_components = [\n",
+ " 0,\n",
+ " transfer_deltaV * np.cos(np.deg2rad(leo_orbit.inc.value)),\n",
+ " transfer_deltaV * np.sin(np.deg2rad(leo_orbit.inc.value))\n",
+ " ]\n",
+ " dv_vector = np.array(dv_components) * u.km / u.s\n",
+ "\n",
+ " # Apply first maneuver to raise orbit from LEO to GTO\n",
+ " impulse_maneuver = Maneuver.impulse(dv_vector)\n",
+ " gto_orbit = leo_orbit.apply_maneuver(impulse_maneuver)\n",
+ "\n",
+ " # Compute magnitude of applied delta-V (in km/s)\n",
+ " mag_deltaV = np.linalg.norm(dv_vector.to(u.km/u.s).value)\n",
+ "\n",
+ " # Optimize rocket mass based on stages using an external function below\n",
+ " stage_results = optimize_rocket_mass(exhaust_vel1, exhaust_vel2, satellite_mass,\n",
+ " mag_deltaV, deltaV_takeoff)\n",
+ " stage1_mass_loss, stage2_mass_loss, wet_stage2_mass, deltaV_stage1, deltaV_stage2 = stage_results\n",
+ "\n",
+ " # Compute updated area-to-mass ratio for drag calculations\n",
+ " local_area_to_mass = ((((np.pi) * (1.85**2)) * (u.m**2)) /\n",
+ " ((satellite_mass + stage2_mass_loss) * u.kg)\n",
+ " ).to(u.km**2 / u.kg).value\n",
+ "\n",
+ " # Propagate the orbit until apogee (half period propagation)\n",
+ " time_to_apogee = (gto_orbit.period.value / 2.0) * u.s\n",
+ " propagated_orbit = gto_orbit.propagate(time_to_apogee, method=CowellPropagator(f=total_dynamics))\n",
+ "\n",
+ " # Calculate velocity for GEO based on norm of propagated position vector\n",
+ " norm_radius = np.linalg.norm(propagated_orbit.r.value)\n",
+ " geo_velocity = np.sqrt(Earth.k.value / 1000.0 * ((2.0/norm_radius) - (1.0/42164.0))) / 1000.0\n",
+ "\n",
+ " # Cancel x-component of velocity due to perturbations with a second burn\n",
+ " dv_cancel = np.array([-propagated_orbit.v.value[0], 0, 0]) * (u.m / u.s)\n",
+ " impulse_cancel = Maneuver.impulse(dv_cancel)\n",
+ " propagated_orbit_corrected = propagated_orbit.apply_maneuver(impulse_cancel)\n",
+ "\n",
+ " # Compute delta-V for second burn (aiming to circularize GEO orbit)\n",
+ " mag_deltaV_plane = np.sqrt(np.linalg.norm(propagated_orbit_corrected.v.value)**2 +\n",
+ " geo_velocity**2 -\n",
+ " 2 * geo_velocity *\n",
+ " np.linalg.norm(propagated_orbit_corrected.v.value) *\n",
+ " np.cos(propagated_orbit_corrected.inc.value *\n",
+ " (inc_change_percent/100)))\n",
+ "\n",
+ " # Determine flight path angle based on circularization maneuver requirements\n",
+ " flight_path_angle_deg = (180 - np.degrees(np.arcsin((geo_velocity *\n",
+ " np.sin(propagated_orbit_corrected.inc.value * (inc_change_percent/100)) /\n",
+ " mag_deltaV_plane))))\n",
+ "\n",
+ " # Calculate the components of the delta-V vector using different approaches based on inclination\n",
+ " if leo_orbit.inc.value >= 59 and inc_change_percent == 100:\n",
+ " vector1_magnitude = np.abs(mag_deltaV_plane * np.cos(np.radians(flight_path_angle_deg - 90)))\n",
+ " else:\n",
+ " vector1_magnitude = np.abs(mag_deltaV_plane * np.cos(np.radians(flight_path_angle_deg)))\n",
+ "\n",
+ " angle_vector1 = np.arctan2(propagated_orbit_corrected.v.value[2],\n",
+ " propagated_orbit_corrected.v.value[1])\n",
+ " vector1 = [0, vector1_magnitude * np.cos(angle_vector1),\n",
+ " vector1_magnitude * np.sin(angle_vector1)]\n",
+ "\n",
+ " if leo_orbit.inc.value >= 59 and inc_change_percent == 100:\n",
+ " vector2_magnitude = np.abs(mag_deltaV_plane * np.sin(np.radians(flight_path_angle_deg - 90)))\n",
+ " else:\n",
+ " vector2_magnitude = np.abs(mag_deltaV_plane * np.sin(np.radians(flight_path_angle_deg)))\n",
+ "\n",
+ " # Calculate perpendicular vector to vector1 in the orbital plane\n",
+ " vector1_perp = [0, propagated_orbit_corrected.v.value[2], -propagated_orbit_corrected.v.value[1]]\n",
+ " angle_vector2 = np.arctan2(vector1_perp[2], vector1_perp[1])\n",
+ " vector2 = [0, vector2_magnitude * np.cos(angle_vector2),\n",
+ " vector2_magnitude * np.sin(angle_vector2)]\n",
+ "\n",
+ " # Combine vectors to produce the required GEO plane-change maneuver\n",
+ " if leo_orbit.inc.value >= 59 and inc_change_percent == 100:\n",
+ " combined_dv_vector = [0, -(vector1[2] + vector2[2]),\n",
+ " vector1[1] + vector2[1]]\n",
+ " else:\n",
+ " combined_dv_vector = [0, -(vector1[1] + vector2[1]),\n",
+ " -(vector1[2] + vector2[2])]\n",
+ "\n",
+ " impulse_plane_change = Maneuver.impulse(np.array(combined_dv_vector) * (u.km/u.s))\n",
+ " geo_orbit = propagated_orbit_corrected.apply_maneuver(impulse_plane_change)\n",
+ "\n",
+ " # Total combined delta-V for plane change maneuver (including cancellation burn)\n",
+ " combined_plane_dv = mag_deltaV_plane + np.abs(propagated_orbit.v.value[0])\n",
+ "\n",
+ " # Propellant loss calculations using the rocket equation (staged)\n",
+ " loss_stage3 = satellite_mass - (satellite_mass / np.e**((mag_deltaV_plane * 1000.0) / exhaust_vel_prop))\n",
+ " loss_stage2 = (satellite_mass + stage2_mass_loss) * (np.e**((mag_deltaV * 1000.0) / (348.0 * 9.8)) - 1.0)\n",
+ " loss_stage1 = (satellite_mass + stage2_mass_loss + loss_stage2) * (np.e**((deltaV_stage2 - (mag_deltaV * 1000.0)) / (348.0 * 9.8)) - 1.0)\n",
+ " loss_initial = (satellite_mass + wet_stage2_mass + stage1_mass_loss) * (np.e**(deltaV_stage1 / (300.0 * 9.8)) - 1.0)\n",
+ " total_propellant_loss = loss_initial + loss_stage1 + loss_stage2 + loss_stage3\n",
+ "\n",
+ " return mag_deltaV, combined_plane_dv, total_propellant_loss"
+ ],
+ "metadata": {
+ "id": "T8OKwYHTNb11"
+ },
+ "execution_count": 22,
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "# Simulating the orbit\n",
+ "\n",
+ "This function utilizes all the functions described above to obtain the total fuel mass and delta-V for all maneuevers. It also computes an approximate of the delta-V required to insert the rocket into LEO orbit, accounting for drag, gravity, earth's boost, etc."
+ ],
+ "metadata": {
+ "id": "lYtITy71UGlW"
+ }
+ },
+ {
+ "cell_type": "code",
+ "source": [
+ "# ===============================================================\n",
+ "# Simulation Functions\n",
+ "# ===============================================================\n",
+ "def simulate_orbit(sim_params):\n",
+ " \"\"\"\n",
+ " Perform a single orbit simulation.\n",
+ "\n",
+ " Parameters:\n",
+ " sim_params (tuple): Contains the following parameters:\n",
+ " - inclination_deg: LEO orbital inclination in degrees.\n",
+ " - leo_altitude_km: LEO altitude above Earth surface (km).\n",
+ " - exhaust_vel_stage1: Exhaust velocity for stage 1.\n",
+ " - exhaust_vel_stage2: Exhaust velocity for stage 2.\n",
+ " - exhaust_vel_propulsive: Propulsive exhaust velocity.\n",
+ " - satellite_mass: Dry mass of satellite.\n",
+ " - inclination_change_percent: Percent value for inclination change.\n",
+ "\n",
+ " Returns:\n",
+ " tuple: (inclination_deg, inclination_change_percent, leo_altitude_km,\n",
+ " exhaust_vel_stage1, exhaust_vel_stage2, exhaust_vel_propulsive,\n",
+ " satellite_mass, f_total, total_deltaV)\n",
+ " \"\"\"\n",
+ " # Unpack simulation parameters\n",
+ " (inclination_deg, leo_altitude_km,\n",
+ " exhaust_vel_stage1, exhaust_vel_stage2, exhaust_vel_propulsive,\n",
+ " satellite_mass, inc_change_percent) = sim_params\n",
+ "\n",
+ " # Create an initial LEO orbit object using given altitude and inclination\n",
+ " leo_orbit = create_leo_orbit(leo_altitude_km, inclination_deg)\n",
+ "\n",
+ " # Compute delta-V components from orbital mechanics and perturbations\n",
+ " transfer_dv = np.sqrt(Earth.k.value / (6378000 + (leo_altitude_km * 1000)))\n",
+ " potential_dv = (np.sqrt(Earth.k.value / 6378000) *\n",
+ " np.sqrt(2 - (6378000 / (6378000 + (leo_altitude_km * 1000)))) -\n",
+ " np.sqrt(Earth.k.value / (6378000 + (leo_altitude_km * 1000))))\n",
+ " grav_loss_dv = np.interp(leo_altitude_km, LEO_ALTITUDES_KM, DELTA_V_GRAV_KMS)\n",
+ " atm_loss_dv = np.interp(leo_altitude_km, LEO_ALTITUDES_KM, DELTA_V_ATM_KMS)\n",
+ " earth_boost_dv = (7.2921159e-5) * Earth.R.to(u.m).value * np.cos(leo_orbit.inc.to(u.rad).value)\n",
+ "\n",
+ " # Total take-off delta-V is the sum of various components minus the boost provided by Earth's rotation\n",
+ " deltaV_takeoff = transfer_dv + grav_loss_dv + potential_dv + atm_loss_dv - earth_boost_dv\n",
+ "\n",
+ " # Compute additional delta-V contributions by propagating the orbit and other calculations.\n",
+ " mag_deltaV, mag_deltaV_plane, total_propulsive_loss = propagate_and_compute_transfer_dv(\n",
+ " leo_orbit, exhaust_vel_stage1, exhaust_vel_stage2, exhaust_vel_propulsive,\n",
+ " satellite_mass, deltaV_takeoff, inc_change_percent\n",
+ " )\n",
+ "\n",
+ " total_deltaV = deltaV_takeoff + mag_deltaV * 1000 + mag_deltaV_plane * 1000\n",
+ "\n",
+ " # Return complete simulation results\n",
+ " return (inclination_deg, inc_change_percent, leo_altitude_km,\n",
+ " exhaust_vel_stage1, exhaust_vel_stage2, exhaust_vel_propulsive,\n",
+ " satellite_mass, total_propulsive_loss, total_deltaV)"
+ ],
+ "metadata": {
+ "id": "Ldu_y1EJUGCD"
+ },
+ "execution_count": 23,
+ "outputs": []
+ }
+ ]
+}
\ No newline at end of file
diff --git a/Jupyter_tutorial_ARC.ipynb b/Jupyter_tutorial_ARC.ipynb
new file mode 100644
index 0000000..e482f89
--- /dev/null
+++ b/Jupyter_tutorial_ARC.ipynb
@@ -0,0 +1,802 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "view-in-github",
+ "colab_type": "text"
+ },
+ "source": [
+ ""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "ZJ2210Cdwwal"
+ },
+ "source": [
+ "# Analyzing orbital perturbations of TESS orbit\n",
+ "\n",
+ " The goal of this tutorial is to simulate the orbit of the Transiting Exoplanet Survey Satellite (TESS) under the effect of different perturbing forces (solar radiation, gravitational forces, etc.). The simulation will be performed for approximately 8.5 years. Therefore, we aim to understand how we can minimize computational processing time by compromising some degree of accuracy in our simulation. This is achieved by ignoring the effect of each perturbing force and computing the evolution of the fractional error of the position and velocity of TESS over time. The JPL Horizons database is utilized to obtain our reference values for the position and velocity of TESS for the timespan of the simulation. By understanding which perturbing forces can be dismissed while maintaining a low fractional error, we can generalize the results of this simulation to others which have a similar orbit and thus save computational time.\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true,
+ "id": "sLFju77hwW6w"
+ },
+ "outputs": [],
+ "source": [
+ "#Install necessary packages for Google Colab. Note: hapsira is the modern version of poliastro to handle Python 3.12\n",
+ "!pip install hapsira \"astropy<7.0\" jplephem"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "fqRzoMKXwu_r"
+ },
+ "source": [
+ "# Code Imports\n",
+ "This code utilizes several Python libraries for different functionalities:\n",
+ "1. poliastro: defines celestial bodies and physical constants, calculates the evolution of the orbit using Cowell Propagator, and introduces the effects of solar radiation and other forces\n",
+ "\n",
+ "2. astropy: imports the JPL horizons database and handles coordinate system conversions to accurately compare the simulation results\n",
+ "\n",
+ "3. numpy: performs math computations\n",
+ "\n",
+ "4. matplotlib: plots evolution of fractional error for position and velocity over time in an organized graph with axes labels and legend\n",
+ "\n",
+ "5. gc: clean up data\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "NI5Hi5Tj7kF9"
+ },
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt\n",
+ "import astropy.units as u\n",
+ "from astropy.time import Time\n",
+ "from astropy.coordinates import solar_system_ephemeris\n",
+ "from hapsira.bodies import Earth, Moon, Sun, Jupiter, Venus\n",
+ "from hapsira.twobody import Orbit\n",
+ "from hapsira.twobody.propagation import CowellPropagator\n",
+ "from hapsira.core.perturbations import third_body, radiation_pressure\n",
+ "from hapsira.core.propagation import func_twobody\n",
+ "from hapsira.ephem import build_ephem_interpolant, Ephem\n",
+ "from hapsira.constants import rho0_earth, H0_earth\n",
+ "import time\n",
+ "import gc\n",
+ "from astropy.coordinates import (\n",
+ " CartesianRepresentation, CartesianDifferential, SkyCoord, ICRS, GCRS\n",
+ ")\n",
+ "\n",
+ "# ---------------------- Simulation Configuration ----------------------\n",
+ "# Total simulation duration in days.\n",
+ "SIMULATION_DURATION_DAYS = 3050"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "34eJzUCH-WOg"
+ },
+ "source": [
+ "# Environment Initialization\n",
+ "\n",
+ "Ths function defines the starting time of the simulation as June 1st, 2018. It then creates an evenly-spaced position vector for each celestial body with its current position every 1 hour from the beginning to the end of the simulation. This is more effective than looking for the exact position of each celestial body for each small step in the simulation since it would significantly increase the computational time."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "4EjAFUmaAaGw"
+ },
+ "outputs": [],
+ "source": [
+ "def initialize_environment():\n",
+ " \"\"\"\n",
+ " Set up the solar system ephemeris and build ephemeris interpolants.\n",
+ " \"\"\"\n",
+ " # Choose ephemeris data source (JPL in this case)\n",
+ " solar_system_ephemeris.set(\"jpl\")\n",
+ "\n",
+ " # 1. Define simulation epoch as a pure Time object\n",
+ " simulation_start_epoch = Time(\"2018-06-01T00:00:00\", scale=\"utc\", format=\"isot\")\n",
+ "\n",
+ " # 2. Calculate the exact end Time object by adding the days directly to the Time object (NO .jd!)\n",
+ " simulation_end_epoch = simulation_start_epoch + SIMULATION_DURATION_DAYS * u.day\n",
+ "\n",
+ " # Generate an array of epochs at 1-hour intervals for the interpolation\n",
+ " time_span_delta = simulation_end_epoch - simulation_start_epoch\n",
+ " num_steps = int(time_span_delta.to(u.hr).value) # Calculate number of steps to ensure including the end epoch\n",
+ " epochs_for_interpolant = simulation_start_epoch + np.arange(num_steps + 1) * u.hr\n",
+ "\n",
+ " # 3. Build ephemeris interpolants, specifying Earth as the attractor\n",
+ " moon_ephem_interp = build_ephem_interpolant(\n",
+ " Moon, epochs_for_interpolant, attractor=Earth\n",
+ " )\n",
+ " sun_ephem_interp = build_ephem_interpolant(\n",
+ " Sun, epochs_for_interpolant, attractor=Earth\n",
+ " )\n",
+ " jupiter_ephem_interp = build_ephem_interpolant(\n",
+ " Jupiter, epochs_for_interpolant, attractor=Earth\n",
+ " )\n",
+ " venus_ephem_interp = build_ephem_interpolant(\n",
+ " Venus, epochs_for_interpolant, attractor=Earth\n",
+ " )\n",
+ "\n",
+ " # Clean up memory\n",
+ " gc.collect()\n",
+ "\n",
+ " return simulation_start_epoch, moon_ephem_interp, sun_ephem_interp, jupiter_ephem_interp, venus_ephem_interp"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "xPL8KK8kAhek"
+ },
+ "source": [
+ "# Setting up physical constants\n",
+ "\n",
+ "This function mainly computes the required physical constants to perform the orbital simulations and converts them to the correct units"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "bzZ33Af7BMOW"
+ },
+ "outputs": [],
+ "source": [
+ "# ---------------------- Physical Constants ----------------------\n",
+ "def define_physical_constants():\n",
+ " \"\"\"\n",
+ " Retrieve and compute key physical constants required for the simulation.\n",
+ "\n",
+ " Returns:\n",
+ " dict: A dictionary containing various physical parameters.\n",
+ " \"\"\"\n",
+ " earth_radius_km = Earth.R.to(u.km).value\n",
+ " gravitational_param = Earth.k.to(u.km**3 / u.s**2).value\n",
+ "\n",
+ " # Atmospheric drag coefficient (negligible for TESS orbit, hence set as a constant)\n",
+ " drag_coefficient = 2.2\n",
+ "\n",
+ " # Area-to-mass ratio, converting from square meters/kg to square km/kg\n",
+ " area_to_mass = ((3) * (u.m**2) / (350 * u.kg)).to_value(u.km**2 / u.kg)\n",
+ "\n",
+ " # Atmospheric density at Earth's surface and scale height\n",
+ " surface_density = rho0_earth.to(u.kg / u.km**3).value\n",
+ " scale_height = H0_earth.to(u.km).value\n",
+ "\n",
+ " # Radiation pressure coefficient and solar flux-to-speed constant ratio\n",
+ " radiation_coefficient = 1.3\n",
+ " solar_flux_over_c = (3.828e20) / 299792.458\n",
+ "\n",
+ " return {\n",
+ " \"earth_radius\": earth_radius_km,\n",
+ " \"grav_param\": gravitational_param,\n",
+ " \"drag_coefficient\": drag_coefficient,\n",
+ " \"area_to_mass\": area_to_mass,\n",
+ " \"surface_density\": surface_density,\n",
+ " \"scale_height\": scale_height,\n",
+ " \"radiation_coefficient\": radiation_coefficient,\n",
+ " \"solar_flux_over_c\": solar_flux_over_c,\n",
+ " }"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "bMhDPYC3BUr9"
+ },
+ "source": [
+ "# Computing Acceleration Vector\n",
+ "\n",
+ "This function creates a 3D vector representing the acceleration in the x, y, and z directions for our spacecraft. It then checks whether it is desired to calculate the effect of a perturbing force by checking the force_flags dictionary. If it is listed as \"True,\" it calls the corresponding function from poliastro (e.g. \"third_body\" for celestial bodies) to add the acceleration produced by such force. In this way, we can bypass the effect of a force to analyze how accurate the simulation can still be. The function returns an acceleration vector as the addition of all the accelerations produced by the desired perturbing forces"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "_iUBJbRcD4l0"
+ },
+ "outputs": [],
+ "source": [
+ "# ---------------------- Perturbation Acceleration Functions ----------------------\n",
+ "def compute_perturbation_acceleration(t0, state, grav_param, constants, sun_interp, moon_interp,\n",
+ " jupiter_interp, venus_interp, force_flags):\n",
+ " \"\"\"\n",
+ " Computes the total perturbation acceleration due to third-body effects and radiation pressure.\n",
+ "\n",
+ " Args:\n",
+ " t0 (float): Initial time.\n",
+ " state (ndarray): Current state vector of the spacecraft.\n",
+ " grav_param (float): Gravitational parameter of the primary body.\n",
+ " constants (dict): Physical constants (e.g., Earth radius, radiation coefficients).\n",
+ " sun_interp: Ephemeris interpolant for the Sun.\n",
+ " moon_interp: Ephemeris interpolant for the Moon.\n",
+ " jupiter_interp: Ephemeris interpolant for Jupiter.\n",
+ " venus_interp: Ephemeris interpolant for Venus.\n",
+ " force_flags (dict): Dictionary indicating which forces to include in the simulation.\n",
+ "\n",
+ " Returns:\n",
+ " ndarray: Combined acceleration vector due to all selected perturbation forces.\n",
+ " \"\"\"\n",
+ " total_acceleration = np.zeros(3)\n",
+ "\n",
+ " # Moon's third-body effect\n",
+ " if force_flags.get(\"Moon\", True):\n",
+ " total_acceleration += third_body(t0, state, grav_param,\n",
+ " k_third=Moon.k.to(u.km**3/u.s**2).value,\n",
+ " perturbation_body=moon_interp)\n",
+ " # Sun's third-body effect (using sun_interp as the star body)\n",
+ " if force_flags.get(\"Sun\", True):\n",
+ " total_acceleration += third_body(t0, state, grav_param,\n",
+ " k_third=Sun.k.to(u.km**3/u.s**2).value,\n",
+ " perturbation_body=sun_interp)\n",
+ " # Jupiter's third-body effect\n",
+ " if force_flags.get(\"Jupiter\", True):\n",
+ " total_acceleration += third_body(t0, state, grav_param,\n",
+ " k_third=Jupiter.k.to(u.km**3/u.s**2).value,\n",
+ " perturbation_body=jupiter_interp)\n",
+ " # Venus's third-body effect\n",
+ " if force_flags.get(\"Venus\", True):\n",
+ " total_acceleration += third_body(t0, state, grav_param,\n",
+ " k_third=Venus.k.to(u.km**3/u.s**2).value,\n",
+ " perturbation_body=venus_interp)\n",
+ " # Solar radiation pressure effect\n",
+ " if force_flags.get(\"Radiation\", True):\n",
+ " total_acceleration += radiation_pressure(\n",
+ " t0, state, grav_param,\n",
+ " constants[\"earth_radius\"],\n",
+ " constants[\"radiation_coefficient\"],\n",
+ " constants[\"area_to_mass\"],\n",
+ " constants[\"solar_flux_over_c\"],\n",
+ " sun_interp\n",
+ " )\n",
+ " return total_acceleration"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "gQVWwLQnD--r"
+ },
+ "source": [
+ "# Computing the derivative of the state vector\n",
+ "\n",
+ "In order to use the Cowell Propagator from poliastro, this function computes how fast the position and velocity of the spacecraft is changing at a specific point in time. It calculates the velocity and acceleration vector based on a 2-body Keplerian motion model between the Earth and the spacecraft. Then, it adds all the third-body effects to this vector to obtain a final vector describing how both the position and velocity are changing through time at that instant."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "WA3qxCw4KLXE"
+ },
+ "outputs": [],
+ "source": [
+ "def perturbation_derivatives(t0, state, grav_param, constants, sun_interp, moon_interp,\n",
+ " jupiter_interp, venus_interp, force_flags):\n",
+ " \"\"\"\n",
+ " Computes the time derivatives of the state vector including both Keplerian and perturbation accelerations.\n",
+ "\n",
+ " Args:\n",
+ " t0 (float): Initial time.\n",
+ " state (ndarray): Current state vector.\n",
+ " grav_param (float): Gravitational parameter for the primary body.\n",
+ " constants (dict): Dictionary of physical constants.\n",
+ " sun_interp, moon_interp, jupiter_interp, venus_interp: Ephemeris interpolants.\n",
+ " force_flags (dict): Flags indicating which perturbing forces are active.\n",
+ "\n",
+ " Returns:\n",
+ " ndarray: Time derivative of the state vector.\n",
+ " \"\"\"\n",
+ " # Two-body (Keplerian) derivative\n",
+ " keplerian_deriv = func_twobody(t0, state, grav_param)\n",
+ "\n",
+ " # Compute total perturbation acceleration.\n",
+ " accel_perturb = compute_perturbation_acceleration(\n",
+ " t0, state, grav_param, constants, sun_interp, moon_interp, jupiter_interp, venus_interp, force_flags\n",
+ " )\n",
+ "\n",
+ " # Append the perturbation acceleration (affecting velocity components) to the Keplerian derivative.\n",
+ " return keplerian_deriv + np.array([0, 0, 0, accel_perturb[0], accel_perturb[1], accel_perturb[2]])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "eiUzvvn5KQDl"
+ },
+ "source": [
+ "# Initializing the orbit\n",
+ "\n",
+ "This function defines the Orbit object needed to start the simulation. To do so, it extracts the the state vector of TESS at the starting time of our simulation from the JPL Horizons database. It performs some necessary unit conversions and translates the vector from ICRS (Sun-centered) to GCRS (Earth-centered) coordinate system, as poliastro uses GCRS as default for orbit propagations."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "_AaJQNupKg2Y"
+ },
+ "outputs": [],
+ "source": [
+ "# ---------------------- Orbit Initialization ----------------------\n",
+ "def initialize_orbit(simulation_epoch):\n",
+ " \"\"\"\n",
+ " Create the initial orbit state from spacecraft ephemeris data (using TESS as reference).\n",
+ " The state vector is transformed from the ICRS frame to the GCRS frame.\n",
+ "\n",
+ " Args:\n",
+ " simulation_epoch (Time): The epoch of the simulation start.\n",
+ "\n",
+ " Returns:\n",
+ " Orbit: The initial Orbit object.\n",
+ " \"\"\"\n",
+ " # Retrieve TESS ephemeris data from Horizons system.\n",
+ " tess_ephem = Ephem.from_horizons(\"TESS\", epochs=simulation_epoch)\n",
+ "\n",
+ " # Get position and velocity vectors in km and km/s respectively.\n",
+ " position_icrf = np.array(tess_ephem.rv()[0].to(u.km).value).reshape(-1)\n",
+ " velocity_icrf = np.array(tess_ephem.rv()[1].to(u.km / u.s).value).reshape(-1)\n",
+ "\n",
+ " # Create initial orbit in ICRF frame from position and velocity.\n",
+ " orbit_icrf = Orbit.from_vectors(Earth, position_icrf * u.km, velocity_icrf * (u.km / u.s), epoch=simulation_epoch)\n",
+ "\n",
+ " # Extract ICRF vectors\n",
+ " r_icrf, v_icrf = orbit_icrf.rv()\n",
+ "\n",
+ " # Transform state vector from ICRS/ICRF to the GCRS frame.\n",
+ " cart_representation = CartesianRepresentation(*r_icrf)\n",
+ " cart_differentials = CartesianDifferential(*v_icrf)\n",
+ " icrs_coord = SkyCoord(cart_representation.with_differentials(cart_differentials),\n",
+ " frame=ICRS, obstime=simulation_epoch)\n",
+ " gcrs_coord = icrs_coord.transform_to(GCRS(obstime=simulation_epoch))\n",
+ "\n",
+ " # Retrieve transformed position and velocity components.\n",
+ " r_gcrs = gcrs_coord.data.represent_as(CartesianRepresentation).xyz\n",
+ " v_gcrs = gcrs_coord.data.differentials[\"s\"].d_xyz\n",
+ "\n",
+ " # Clean up temporary variables.\n",
+ " del cart_representation, cart_differentials, icrs_coord, gcrs_coord\n",
+ " gc.collect()\n",
+ "\n",
+ " return Orbit.from_vectors(Earth, r_gcrs, v_gcrs, epoch=simulation_epoch)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "PRC-BO0IMUfx"
+ },
+ "source": [
+ "# Orbit propagation\n",
+ "\n",
+ "This function simply propagates the starting orbit using poliastro's Cowell Propagator over the specified time interval. It returns the final orbit object after the simulation, from which we can extract the state vector. In principle, this function is called for each time step during the simulation to obtain a continuous evolution of the position and velocity of the spacecraft"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "vmGWMSUdMXER"
+ },
+ "outputs": [],
+ "source": [
+ "# ---------------------- Orbit Propagation & Data Processing ----------------------\n",
+ "def propagate_orbit(orbit_obj, time_interval, perturbation_function):\n",
+ " \"\"\"\n",
+ " Propagate the orbit over a specified time interval using the Cowell propagator.\n",
+ "\n",
+ " Args:\n",
+ " orbit_obj (Orbit): The initial Orbit object.\n",
+ " time_interval (Quantity): Duration over which to propagate.\n",
+ " perturbation_function (function): Function to compute perturbations.\n",
+ "\n",
+ " Returns:\n",
+ " Orbit: Final orbit after propagation.\n",
+ " \"\"\"\n",
+ " propagated_orbit = orbit_obj.propagate(\n",
+ " time_interval, method=CowellPropagator(rtol=3e-14, f=perturbation_function)\n",
+ " )\n",
+ " gc.collect()\n",
+ " return propagated_orbit"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "uSBQa23XNQiR"
+ },
+ "source": [
+ "# Creating evenly-spaced time intervals\n",
+ "\n",
+ "This function creates epoch objects which are evenly spaced. The number of intervals is user-defined and can affect the how precise our simulation is. A good value should be chosen so that we can have enough precision while still running the program in a reasonable amount of time. Thus, 100 intervals were chosen."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "yTf9vpShNREu"
+ },
+ "outputs": [],
+ "source": [
+ "def create_simulation_epochs(start_epoch, end_epoch, num_intervals):\n",
+ " \"\"\"\n",
+ " Generate an array of epochs spanning from the start to the end epoch.\n",
+ " \"\"\"\n",
+ " # Safely convert the TimeDelta into raw seconds\n",
+ " total_seconds = (end_epoch - start_epoch).to(u.s).value\n",
+ "\n",
+ " # Generate the checkpoints and add them back to the start time\n",
+ " epochs = start_epoch + np.linspace(0, total_seconds, num_intervals + 1) * u.s\n",
+ "\n",
+ " gc.collect()\n",
+ " return epochs"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "NX7MXI5aOFr6"
+ },
+ "source": [
+ "# Extracting data for error analysis\n",
+ "\n",
+ "This function calls the Cowell Propagator for each time interval to obtain the Orbit object for each time step. Then, it extracts the position and velocity vector from each object. These values are converted from GCRS to the ICRS coordinate system for comparison with the JPL Horizons database"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "JcRYLKiZOEcQ"
+ },
+ "outputs": [],
+ "source": [
+ "def simulate_spacecraft_trajectory(orbit_obj, perturbation_function, epoch_array):\n",
+ " \"\"\"\n",
+ " Propagate the spacecraft orbit for each epoch in the provided array and convert the results\n",
+ " to ICRS frame state vectors.\n",
+ "\n",
+ " Args:\n",
+ " orbit_obj (Orbit): Initial orbit.\n",
+ " perturbation_function (function): Function to compute perturbation accelerations.\n",
+ " epoch_array (Quantity): Array of epochs at which to propagate the orbit.\n",
+ "\n",
+ " Returns:\n",
+ " tuple: (r_state, v_state) state arrays in the ICRS frame.\n",
+ " \"\"\"\n",
+ " start_time = time.perf_counter()\n",
+ " propagated_orbits = [\n",
+ " orbit_obj.propagate(ep, method=CowellPropagator(rtol=3e-14, f=perturbation_function))\n",
+ " for ep in epoch_array\n",
+ " ]\n",
+ " end_time = time.perf_counter()\n",
+ " print(\"Propagation time: {:.3f} seconds\".format(end_time - start_time))\n",
+ "\n",
+ " positions, velocities, epochs_list = [], [], []\n",
+ " for orbit_instance in propagated_orbits:\n",
+ " positions.append(orbit_instance.r)\n",
+ " velocities.append(orbit_instance.v)\n",
+ " epochs_list.append(orbit_instance.epoch)\n",
+ "\n",
+ " # Convert lists to numpy arrays\n",
+ " pos_array = np.array(positions)\n",
+ " vel_array = np.array(velocities)\n",
+ " epochs_array = Time(epochs_list)\n",
+ "\n",
+ " # Convert coordinates from GCRS to ICRS for output and analysis.\n",
+ " conversion_start = time.perf_counter()\n",
+ " sky_coordinates = SkyCoord(\n",
+ " CartesianRepresentation(\n",
+ " pos_array[:, 0] * u.km,\n",
+ " pos_array[:, 1] * u.km,\n",
+ " pos_array[:, 2] * u.km,\n",
+ " differentials=CartesianDifferential(\n",
+ " vel_array[:, 0] * (u.km / u.s),\n",
+ " vel_array[:, 1] * (u.km / u.s),\n",
+ " vel_array[:, 2] * (u.km / u.s)\n",
+ " )\n",
+ " ),\n",
+ " frame=GCRS,\n",
+ " obstime=epochs_array\n",
+ " )\n",
+ " icrs_coordinates = sky_coordinates.transform_to(ICRS)\n",
+ " r_state = icrs_coordinates.cartesian.xyz.T.value\n",
+ " v_state = icrs_coordinates.cartesian.differentials[\"s\"].d_xyz.T.value\n",
+ "\n",
+ " # Clean up and free memory.\n",
+ " del positions, velocities, pos_array, vel_array, epochs_list, epochs_array, sky_coordinates, icrs_coordinates\n",
+ " gc.collect()\n",
+ "\n",
+ " conversion_end = time.perf_counter()\n",
+ " print(\"Coordinate conversion time: {:.3f} seconds\".format(conversion_end - conversion_start))\n",
+ " return r_state, v_state"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "WDbOY6Y8PP0E"
+ },
+ "source": [
+ "# Retrieving TESS data from JPL Horizons\n",
+ "\n",
+ "At each time step, this function retrieves the position and velocity vectors from JPL Horizons for TESS"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "bmRkh3r_PQae"
+ },
+ "outputs": [],
+ "source": [
+ "def retrieve_actual_spacecraft_data(epoch_array, num_intervals):\n",
+ " \"\"\"\n",
+ " Obtain reference spacecraft data from Horizons (using TESS as reference) for error comparison.\n",
+ "\n",
+ " Args:\n",
+ " epoch_array (Quantity): Array of simulation epochs.\n",
+ " num_intervals (int): Number of intervals (subdivisions used during simulation).\n",
+ "\n",
+ " Returns:\n",
+ " tuple: (r_actual, v_actual) state arrays from reference data.\n",
+ " \"\"\"\n",
+ " tess_ephem_list = [Ephem.from_horizons(\"TESS\", epochs=ep) for ep in epoch_array]\n",
+ " actual_positions = np.array([ephem.rv()[0].to(u.km).value for ephem in tess_ephem_list])\n",
+ " actual_velocities = np.array([ephem.rv()[1].to(u.km / u.s).value for ephem in tess_ephem_list])\n",
+ "\n",
+ " # Ensure that the arrays have the expected dimensions.\n",
+ " r_actual = actual_positions.reshape(num_intervals + 1, 3)\n",
+ " v_actual = actual_velocities.reshape(num_intervals + 1, 3)\n",
+ " del tess_ephem_list\n",
+ " gc.collect()\n",
+ " return r_actual, v_actual"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "8vTnimfdPo_m"
+ },
+ "source": [
+ "# Computing and plotting fractional error over time\n",
+ "\n",
+ "These functions calculate and plot the fractional error as a function of time for the duration of the simulation. The fractional error for the position and velocity vectors are computed as the difference between the Euclidean norm between the simulation vector and JPL Horizons vector dividided by the JPL Horizons vector. This is done for each time step using vectorized operations with the help of numpy. The results are plotted in two different figures: one for the position and one for the velocity vector"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "FlquX_O1PpeE"
+ },
+ "outputs": [],
+ "source": [
+ "# ---------------------- Error Analysis ----------------------\n",
+ "def calculate_fractional_errors(simulated_positions, reference_positions, simulated_velocities, reference_velocities):\n",
+ " \"\"\"\n",
+ " Calculate the fractional errors between simulated and reference position and velocity vectors.\n",
+ "\n",
+ " Args:\n",
+ " simulated_positions (ndarray): Simulated positions.\n",
+ " reference_positions (ndarray): Reference positions.\n",
+ " simulated_velocities (ndarray): Simulated velocities.\n",
+ " reference_velocities (ndarray): Reference velocities.\n",
+ "\n",
+ " Returns:\n",
+ " tuple: Fractional errors for positions and velocities.\n",
+ " \"\"\"\n",
+ " # Position error calculations\n",
+ " position_error_norm = np.linalg.norm(simulated_positions - reference_positions, axis=1)\n",
+ " reference_position_norm = np.linalg.norm(reference_positions, axis=1)\n",
+ " fractional_error_position = position_error_norm / reference_position_norm\n",
+ "\n",
+ " # Velocity error calculations\n",
+ " velocity_error_norm = np.linalg.norm(simulated_velocities - reference_velocities, axis=1)\n",
+ " reference_velocity_norm = np.linalg.norm(reference_velocities, axis=1)\n",
+ " fractional_error_velocity = velocity_error_norm / reference_velocity_norm\n",
+ "\n",
+ " return fractional_error_position, fractional_error_velocity\n",
+ "\n",
+ "# ---------------------- Plotting Functions ----------------------\n",
+ "def plot_fractional_errors(time_years, position_errors_dict, velocity_errors_dict):\n",
+ " \"\"\"\n",
+ " Plot and save the fractional error comparisons for positions and velocities.\n",
+ "\n",
+ " Args:\n",
+ " time_years (ndarray): Time axis in years.\n",
+ " position_errors_dict (dict): Dictionary of position error arrays keyed by simulation case label.\n",
+ " velocity_errors_dict (dict): Dictionary of velocity error arrays keyed by simulation case label.\n",
+ " \"\"\"\n",
+ " # Plot position fractional errors.\n",
+ " fig_pos, ax_pos = plt.subplots(figsize=(10, 6))\n",
+ " for label, error_array in position_errors_dict.items():\n",
+ " ax_pos.plot(time_years, error_array, 'o-', label=label)\n",
+ " ax_pos.set_title('Fractional Error in Position Magnitude')\n",
+ " ax_pos.set_xlabel('Time Since Start (years)')\n",
+ " ax_pos.set_ylabel('Fractional Error')\n",
+ " ax_pos.legend()\n",
+ " fig_pos.savefig(\"position_magnitude_error.png\", dpi=130)\n",
+ "\n",
+ " # Plot velocity fractional errors.\n",
+ " fig_vel, ax_vel = plt.subplots(figsize=(10, 6))\n",
+ " for label, error_array in velocity_errors_dict.items():\n",
+ " ax_vel.plot(time_years, error_array, 'o-', label=label)\n",
+ " ax_vel.set_title('Fractional Error in Velocity Magnitude')\n",
+ " ax_vel.set_xlabel('Time Since Start (years)')\n",
+ " ax_vel.set_ylabel('Fractional Error')\n",
+ " ax_vel.legend()\n",
+ " fig_vel.savefig(\"velocity_magnitude_error.png\", dpi=130)\n",
+ "\n",
+ " plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "-_FqFJ3-RMho"
+ },
+ "source": [
+ "# Main program function\n",
+ "\n",
+ "This is the main function of the program that handles all other functions. Its main purpose is to update the force_flags dictionary each time with a new force turned \"off\" and plot the results for that corresponding case"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "DyzoIjwGRM_F",
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "# ---------------------- Main Simulation Routine ----------------------\n",
+ "def main_simulation():\n",
+ " # Initialize simulation environment and physical constants.\n",
+ " sim_epoch, moon_interp, sun_interp, jupiter_interp, venus_interp = initialize_environment()\n",
+ " physical_constants = define_physical_constants()\n",
+ "\n",
+ " # Initialize the spacecraft orbit using reference (TESS) data.\n",
+ " initial_orbit = initialize_orbit(sim_epoch)\n",
+ "\n",
+ " # Define the number of intervals for simulation sampling.\n",
+ " NUM_INTERVALS = 100\n",
+ "\n",
+ " # Propagate orbit without any modifications to define final simulation end epoch.\n",
+ " final_orbit = propagate_orbit(\n",
+ " initial_orbit,\n",
+ " SIMULATION_DURATION_DAYS * u.day,\n",
+ " lambda t0, state, k: perturbation_derivatives(\n",
+ " t0, state, k, physical_constants, sun_interp, moon_interp, jupiter_interp,\n",
+ " venus_interp, {\"Moon\": True, \"Sun\": True, \"Jupiter\": True, \"Venus\": True, \"Radiation\": True}\n",
+ " )\n",
+ " )\n",
+ "\n",
+ " # Generate an array of epochs between the start and final orbit epochs.\n",
+ " simulation_epochs = create_simulation_epochs(initial_orbit.epoch, final_orbit.epoch, NUM_INTERVALS)\n",
+ "\n",
+ " # Retrieve reference spacecraft data for error analysis.\n",
+ " reference_positions, reference_velocities = retrieve_actual_spacecraft_data(simulation_epochs, NUM_INTERVALS)\n",
+ "\n",
+ " # Define base force flags and simulation cases (each turning off one perturbation).\n",
+ " base_force_flags = {\n",
+ " \"Moon\": True, \"Sun\": True, \"Jupiter\": True, \"Venus\": True, \"Radiation\": True\n",
+ " }\n",
+ " simulation_cases = {\n",
+ " \"Without Moon\": {**base_force_flags, \"Moon\": False},\n",
+ " \"Without Sun\": {**base_force_flags, \"Sun\": False},\n",
+ " \"Without Jupiter\": {**base_force_flags, \"Jupiter\": False},\n",
+ " \"Without Venus\": {**base_force_flags, \"Venus\": False},\n",
+ " \"Without Radiation\": {**base_force_flags, \"Radiation\": False},\n",
+ " \"All Forces\": base_force_flags.copy()\n",
+ " }\n",
+ "\n",
+ " # Dictionaries to hold the fractional errors for each simulation case.\n",
+ " pos_errors_dict = {}\n",
+ " vel_errors_dict = {}\n",
+ " # Convert epochs difference to years for x-axis plotting.\n",
+ " time_axis_years = (simulation_epochs - simulation_epochs[0]).sec / (3600 * 24 * 365.25)\n",
+ "\n",
+ " # Execute simulation for each force configuration.\n",
+ " for case_label, force_config in simulation_cases.items():\n",
+ " # Define the perturbation function using current force configuration.\n",
+ " perturbation_function = lambda t0, state, k, flags=force_config: perturbation_derivatives(\n",
+ " t0, state, k, physical_constants, sun_interp, moon_interp, jupiter_interp, venus_interp, flags\n",
+ " )\n",
+ " simulated_positions, simulated_velocities = simulate_spacecraft_trajectory(initial_orbit, perturbation_function, simulation_epochs)\n",
+ "\n",
+ " # Compute fractional errors between simulation and reference data.\n",
+ " frac_err_pos, frac_err_vel = calculate_fractional_errors(simulated_positions, reference_positions,\n",
+ " simulated_velocities, reference_velocities)\n",
+ " pos_errors_dict[case_label] = frac_err_pos\n",
+ " vel_errors_dict[case_label] = frac_err_vel\n",
+ "\n",
+ " # Plot and save error comparisons.\n",
+ " plot_fractional_errors(time_axis_years, pos_errors_dict, vel_errors_dict)\n",
+ " gc.collect()\n",
+ "\n",
+ "# Execute the simulation only if this script is run directly.\n",
+ "if __name__ == \"__main__\":\n",
+ " main_simulation()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "D3x_mHN-Q6wr"
+ },
+ "source": [
+ "The following figures are obtained from this code:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "g1xBm_YEQkQ5"
+ },
+ "source": [
+ ""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "tQ0tBVlpRG3S"
+ },
+ "source": [
+ "#"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "# Results and Conclusions\n",
+ "\n",
+ "Both graphs show the need to include the third-body effects from the Moon and the Sun to obtain accurate simulations. The fractional error for the magnitude of both the position and the velocity spikes at a very early stage of the simulation when the gravitational forces from these bodies is ignored. The solar radiation force becomes more relevant at around 6 years after the beginning of the simulation, where it starts to clearly diverge from the \"All forces\" case. In contrast, removing the influence of secondary celestial bodies, such as Jupiter and Venus, on the orbit of the simulated spacecraft shows no clear sign of loss of accuracy. Based on these results, it can be concluded that when simulating the orbit of a spacecraft on High-Earth Orbit (HEO), such as TESS, it is crucial to include the gravitational effects of primary celestial bodies like the Moon and the Sun while the influence from planets like Jupiter and Venus can be removed, saving computational time while not compromising a significant loss of accuracy\n",
+ "\n"
+ ],
+ "metadata": {
+ "id": "rTMJYPxyqxQB"
+ }
+ }
+ ],
+ "metadata": {
+ "colab": {
+ "provenance": [],
+ "authorship_tag": "ABX9TyMsVwy+oufGfJ1JcOW2+zMX",
+ "include_colab_link": true
+ },
+ "kernelspec": {
+ "display_name": "Python 3",
+ "name": "python3"
+ },
+ "language_info": {
+ "name": "python"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}
\ No newline at end of file
diff --git a/TESS_Simulation.ipynb b/TESS_Simulation.ipynb
new file mode 100644
index 0000000..e1f0912
--- /dev/null
+++ b/TESS_Simulation.ipynb
@@ -0,0 +1,802 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "view-in-github",
+ "colab_type": "text"
+ },
+ "source": [
+ ""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "ZJ2210Cdwwal"
+ },
+ "source": [
+ "# Analyzing orbital perturbations of TESS orbit\n",
+ "\n",
+ " The goal of this tutorial is to simulate the orbit of the Transiting Exoplanet Survey Satellite (TESS) under the effect of different perturbing forces (solar radiation, gravitational forces, etc.). The simulation will be performed for approximately 8.5 years. Therefore, we aim to understand how we can minimize computational processing time by compromising some degree of accuracy in our simulation. This is achieved by ignoring the effect of each perturbing force and computing the evolution of the fractional error of the position and velocity of TESS over time. The JPL Horizons database is utilized to obtain our reference values for the position and velocity of TESS for the timespan of the simulation. By understanding which perturbing forces can be dismissed while maintaining a low fractional error, we can generalize the results of this simulation to others which have a similar orbit and thus save computational time.\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true,
+ "id": "sLFju77hwW6w"
+ },
+ "outputs": [],
+ "source": [
+ "#Install necessary packages for Google Colab. Note: hapsira is the modern version of poliastro to handle Python 3.12\n",
+ "!pip install hapsira \"astropy<7.0\" jplephem"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "fqRzoMKXwu_r"
+ },
+ "source": [
+ "# Code Imports\n",
+ "This code utilizes several Python libraries for different functionalities:\n",
+ "1. poliastro: defines celestial bodies and physical constants, calculates the evolution of the orbit using Cowell Propagator, and introduces the effects of solar radiation and other forces\n",
+ "\n",
+ "2. astropy: imports the JPL horizons database and handles coordinate system conversions to accurately compare the simulation results\n",
+ "\n",
+ "3. numpy: performs math computations\n",
+ "\n",
+ "4. matplotlib: plots evolution of fractional error for position and velocity over time in an organized graph with axes labels and legend\n",
+ "\n",
+ "5. gc: clean up data\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "NI5Hi5Tj7kF9"
+ },
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt\n",
+ "import astropy.units as u\n",
+ "from astropy.time import Time\n",
+ "from astropy.coordinates import solar_system_ephemeris\n",
+ "from hapsira.bodies import Earth, Moon, Sun, Jupiter, Venus\n",
+ "from hapsira.twobody import Orbit\n",
+ "from hapsira.twobody.propagation import CowellPropagator\n",
+ "from hapsira.core.perturbations import third_body, radiation_pressure\n",
+ "from hapsira.core.propagation import func_twobody\n",
+ "from hapsira.ephem import build_ephem_interpolant, Ephem\n",
+ "from hapsira.constants import rho0_earth, H0_earth\n",
+ "import time\n",
+ "import gc\n",
+ "from astropy.coordinates import (\n",
+ " CartesianRepresentation, CartesianDifferential, SkyCoord, ICRS, GCRS\n",
+ ")\n",
+ "\n",
+ "# ---------------------- Simulation Configuration ----------------------\n",
+ "# Total simulation duration in days.\n",
+ "SIMULATION_DURATION_DAYS = 3050"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "34eJzUCH-WOg"
+ },
+ "source": [
+ "# Environment Initialization\n",
+ "\n",
+ "Ths function defines the starting time of the simulation as June 1st, 2018. It then creates an evenly-spaced position vector for each celestial body with its current position every 1 hour from the beginning to the end of the simulation. This is more effective than looking for the exact position of each celestial body for each small step in the simulation since it would significantly increase the computational time."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "4EjAFUmaAaGw"
+ },
+ "outputs": [],
+ "source": [
+ "def initialize_environment():\n",
+ " \"\"\"\n",
+ " Set up the solar system ephemeris and build ephemeris interpolants.\n",
+ " \"\"\"\n",
+ " # Choose ephemeris data source (JPL in this case)\n",
+ " solar_system_ephemeris.set(\"jpl\")\n",
+ "\n",
+ " # 1. Define simulation epoch as a pure Time object\n",
+ " simulation_start_epoch = Time(\"2018-06-01T00:00:00\", scale=\"utc\", format=\"isot\")\n",
+ "\n",
+ " # 2. Calculate the exact end Time object by adding the days directly to the Time object (NO .jd!)\n",
+ " simulation_end_epoch = simulation_start_epoch + SIMULATION_DURATION_DAYS * u.day\n",
+ "\n",
+ " # Generate an array of epochs at 1-hour intervals for the interpolation\n",
+ " time_span_delta = simulation_end_epoch - simulation_start_epoch\n",
+ " num_steps = int(time_span_delta.to(u.hr).value) # Calculate number of steps to ensure including the end epoch\n",
+ " epochs_for_interpolant = simulation_start_epoch + np.arange(num_steps + 1) * u.hr\n",
+ "\n",
+ " # 3. Build ephemeris interpolants, specifying Earth as the attractor\n",
+ " moon_ephem_interp = build_ephem_interpolant(\n",
+ " Moon, epochs_for_interpolant, attractor=Earth\n",
+ " )\n",
+ " sun_ephem_interp = build_ephem_interpolant(\n",
+ " Sun, epochs_for_interpolant, attractor=Earth\n",
+ " )\n",
+ " jupiter_ephem_interp = build_ephem_interpolant(\n",
+ " Jupiter, epochs_for_interpolant, attractor=Earth\n",
+ " )\n",
+ " venus_ephem_interp = build_ephem_interpolant(\n",
+ " Venus, epochs_for_interpolant, attractor=Earth\n",
+ " )\n",
+ "\n",
+ " # Clean up memory\n",
+ " gc.collect()\n",
+ "\n",
+ " return simulation_start_epoch, moon_ephem_interp, sun_ephem_interp, jupiter_ephem_interp, venus_ephem_interp"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "xPL8KK8kAhek"
+ },
+ "source": [
+ "# Setting up physical constants\n",
+ "\n",
+ "This function mainly computes the required physical constants to perform the orbital simulations and converts them to the correct units"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "bzZ33Af7BMOW"
+ },
+ "outputs": [],
+ "source": [
+ "# ---------------------- Physical Constants ----------------------\n",
+ "def define_physical_constants():\n",
+ " \"\"\"\n",
+ " Retrieve and compute key physical constants required for the simulation.\n",
+ "\n",
+ " Returns:\n",
+ " dict: A dictionary containing various physical parameters.\n",
+ " \"\"\"\n",
+ " earth_radius_km = Earth.R.to(u.km).value\n",
+ " gravitational_param = Earth.k.to(u.km**3 / u.s**2).value\n",
+ "\n",
+ " # Atmospheric drag coefficient (negligible for TESS orbit, hence set as a constant)\n",
+ " drag_coefficient = 2.2\n",
+ "\n",
+ " # Area-to-mass ratio, converting from square meters/kg to square km/kg\n",
+ " area_to_mass = ((3) * (u.m**2) / (350 * u.kg)).to_value(u.km**2 / u.kg)\n",
+ "\n",
+ " # Atmospheric density at Earth's surface and scale height\n",
+ " surface_density = rho0_earth.to(u.kg / u.km**3).value\n",
+ " scale_height = H0_earth.to(u.km).value\n",
+ "\n",
+ " # Radiation pressure coefficient and solar flux-to-speed constant ratio\n",
+ " radiation_coefficient = 1.3\n",
+ " solar_flux_over_c = (3.828e20) / 299792.458\n",
+ "\n",
+ " return {\n",
+ " \"earth_radius\": earth_radius_km,\n",
+ " \"grav_param\": gravitational_param,\n",
+ " \"drag_coefficient\": drag_coefficient,\n",
+ " \"area_to_mass\": area_to_mass,\n",
+ " \"surface_density\": surface_density,\n",
+ " \"scale_height\": scale_height,\n",
+ " \"radiation_coefficient\": radiation_coefficient,\n",
+ " \"solar_flux_over_c\": solar_flux_over_c,\n",
+ " }"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "bMhDPYC3BUr9"
+ },
+ "source": [
+ "# Computing Acceleration Vector\n",
+ "\n",
+ "This function creates a 3D vector representing the acceleration in the x, y, and z directions for our spacecraft. It then checks whether it is desired to calculate the effect of a perturbing force by checking the force_flags dictionary. If it is listed as \"True,\" it calls the corresponding function from poliastro (e.g. \"third_body\" for celestial bodies) to add the acceleration produced by such force. In this way, we can bypass the effect of a force to analyze how accurate the simulation can still be. The function returns an acceleration vector is the addition of all the accelerations produced by the desired perturbing forces"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "_iUBJbRcD4l0"
+ },
+ "outputs": [],
+ "source": [
+ "# ---------------------- Perturbation Acceleration Functions ----------------------\n",
+ "def compute_perturbation_acceleration(t0, state, grav_param, constants, sun_interp, moon_interp,\n",
+ " jupiter_interp, venus_interp, force_flags):\n",
+ " \"\"\"\n",
+ " Computes the total perturbation acceleration due to third-body effects and radiation pressure.\n",
+ "\n",
+ " Args:\n",
+ " t0 (float): Initial time.\n",
+ " state (ndarray): Current state vector of the spacecraft.\n",
+ " grav_param (float): Gravitational parameter of the primary body.\n",
+ " constants (dict): Physical constants (e.g., Earth radius, radiation coefficients).\n",
+ " sun_interp: Ephemeris interpolant for the Sun.\n",
+ " moon_interp: Ephemeris interpolant for the Moon.\n",
+ " jupiter_interp: Ephemeris interpolant for Jupiter.\n",
+ " venus_interp: Ephemeris interpolant for Venus.\n",
+ " force_flags (dict): Dictionary indicating which forces to include in the simulation.\n",
+ "\n",
+ " Returns:\n",
+ " ndarray: Combined acceleration vector due to all selected perturbation forces.\n",
+ " \"\"\"\n",
+ " total_acceleration = np.zeros(3)\n",
+ "\n",
+ " # Moon's third-body effect\n",
+ " if force_flags.get(\"Moon\", True):\n",
+ " total_acceleration += third_body(t0, state, grav_param,\n",
+ " k_third=Moon.k.to(u.km**3/u.s**2).value,\n",
+ " perturbation_body=moon_interp)\n",
+ " # Sun's third-body effect (using sun_interp as the star body)\n",
+ " if force_flags.get(\"Sun\", True):\n",
+ " total_acceleration += third_body(t0, state, grav_param,\n",
+ " k_third=Sun.k.to(u.km**3/u.s**2).value,\n",
+ " perturbation_body=sun_interp)\n",
+ " # Jupiter's third-body effect\n",
+ " if force_flags.get(\"Jupiter\", True):\n",
+ " total_acceleration += third_body(t0, state, grav_param,\n",
+ " k_third=Jupiter.k.to(u.km**3/u.s**2).value,\n",
+ " perturbation_body=jupiter_interp)\n",
+ " # Venus's third-body effect\n",
+ " if force_flags.get(\"Venus\", True):\n",
+ " total_acceleration += third_body(t0, state, grav_param,\n",
+ " k_third=Venus.k.to(u.km**3/u.s**2).value,\n",
+ " perturbation_body=venus_interp)\n",
+ " # Solar radiation pressure effect\n",
+ " if force_flags.get(\"Radiation\", True):\n",
+ " total_acceleration += radiation_pressure(\n",
+ " t0, state, grav_param,\n",
+ " constants[\"earth_radius\"],\n",
+ " constants[\"radiation_coefficient\"],\n",
+ " constants[\"area_to_mass\"],\n",
+ " constants[\"solar_flux_over_c\"],\n",
+ " sun_interp\n",
+ " )\n",
+ " return total_acceleration"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "gQVWwLQnD--r"
+ },
+ "source": [
+ "# Computing the derivative of the state vector\n",
+ "\n",
+ "In order to use the Cowell Propagator from poliastro, this function computes how fast the position and velocity of the spacecraft is changing at a specific point in time. It calculates the velocity and acceleration vector based on a 2-body Keplerian motion model between the Earth and the spacecraft. Then, it adds all the third-body effects to this vector to obtain a final vector describing how both the position and velocity are changing through time at that instant."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "WA3qxCw4KLXE"
+ },
+ "outputs": [],
+ "source": [
+ "def perturbation_derivatives(t0, state, grav_param, constants, sun_interp, moon_interp,\n",
+ " jupiter_interp, venus_interp, force_flags):\n",
+ " \"\"\"\n",
+ " Computes the time derivatives of the state vector including both Keplerian and perturbation accelerations.\n",
+ "\n",
+ " Args:\n",
+ " t0 (float): Initial time.\n",
+ " state (ndarray): Current state vector.\n",
+ " grav_param (float): Gravitational parameter for the primary body.\n",
+ " constants (dict): Dictionary of physical constants.\n",
+ " sun_interp, moon_interp, jupiter_interp, venus_interp: Ephemeris interpolants.\n",
+ " force_flags (dict): Flags indicating which perturbing forces are active.\n",
+ "\n",
+ " Returns:\n",
+ " ndarray: Time derivative of the state vector.\n",
+ " \"\"\"\n",
+ " # Two-body (Keplerian) derivative\n",
+ " keplerian_deriv = func_twobody(t0, state, grav_param)\n",
+ "\n",
+ " # Compute total perturbation acceleration.\n",
+ " accel_perturb = compute_perturbation_acceleration(\n",
+ " t0, state, grav_param, constants, sun_interp, moon_interp, jupiter_interp, venus_interp, force_flags\n",
+ " )\n",
+ "\n",
+ " # Append the perturbation acceleration (affecting velocity components) to the Keplerian derivative.\n",
+ " return keplerian_deriv + np.array([0, 0, 0, accel_perturb[0], accel_perturb[1], accel_perturb[2]])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "eiUzvvn5KQDl"
+ },
+ "source": [
+ "# Initializing the orbit\n",
+ "\n",
+ "This function defines the Orbit object needed to start the simulation. To do so, it extracts the the state vector of TESS at the starting time of our simulation from the JPL Horizons database. It performs some necessary unit conversions and translates the vector from ICRS (Sun-centered) to GCRS (Earth-centered) coordinate system, as poliastro uses GCRS as default for orbit propagations."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "_AaJQNupKg2Y"
+ },
+ "outputs": [],
+ "source": [
+ "# ---------------------- Orbit Initialization ----------------------\n",
+ "def initialize_orbit(simulation_epoch):\n",
+ " \"\"\"\n",
+ " Create the initial orbit state from spacecraft ephemeris data (using TESS as reference).\n",
+ " The state vector is transformed from the ICRS frame to the GCRS frame.\n",
+ "\n",
+ " Args:\n",
+ " simulation_epoch (Time): The epoch of the simulation start.\n",
+ "\n",
+ " Returns:\n",
+ " Orbit: The initial Orbit object.\n",
+ " \"\"\"\n",
+ " # Retrieve TESS ephemeris data from Horizons system.\n",
+ " tess_ephem = Ephem.from_horizons(\"TESS\", epochs=simulation_epoch)\n",
+ "\n",
+ " # Get position and velocity vectors in km and km/s respectively.\n",
+ " position_icrf = np.array(tess_ephem.rv()[0].to(u.km).value).reshape(-1)\n",
+ " velocity_icrf = np.array(tess_ephem.rv()[1].to(u.km / u.s).value).reshape(-1)\n",
+ "\n",
+ " # Create initial orbit in ICRF frame from position and velocity.\n",
+ " orbit_icrf = Orbit.from_vectors(Earth, position_icrf * u.km, velocity_icrf * (u.km / u.s), epoch=simulation_epoch)\n",
+ "\n",
+ " # Extract ICRF vectors\n",
+ " r_icrf, v_icrf = orbit_icrf.rv()\n",
+ "\n",
+ " # Transform state vector from ICRS/ICRF to the GCRS frame.\n",
+ " cart_representation = CartesianRepresentation(*r_icrf)\n",
+ " cart_differentials = CartesianDifferential(*v_icrf)\n",
+ " icrs_coord = SkyCoord(cart_representation.with_differentials(cart_differentials),\n",
+ " frame=ICRS, obstime=simulation_epoch)\n",
+ " gcrs_coord = icrs_coord.transform_to(GCRS(obstime=simulation_epoch))\n",
+ "\n",
+ " # Retrieve transformed position and velocity components.\n",
+ " r_gcrs = gcrs_coord.data.represent_as(CartesianRepresentation).xyz\n",
+ " v_gcrs = gcrs_coord.data.differentials[\"s\"].d_xyz\n",
+ "\n",
+ " # Clean up temporary variables.\n",
+ " del cart_representation, cart_differentials, icrs_coord, gcrs_coord\n",
+ " gc.collect()\n",
+ "\n",
+ " return Orbit.from_vectors(Earth, r_gcrs, v_gcrs, epoch=simulation_epoch)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "PRC-BO0IMUfx"
+ },
+ "source": [
+ "# Orbit propagation\n",
+ "\n",
+ "This function simply propagates the starting orbit using poliastro's Cowell Propagator over the specified time interval. It returns the final orbit object after the simulation, from which we can extract the state vector. In principle, this function is called for each time step during the simulation to obtain a continuous evolution of the position and velocity of the spacecraft"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "vmGWMSUdMXER"
+ },
+ "outputs": [],
+ "source": [
+ "# ---------------------- Orbit Propagation & Data Processing ----------------------\n",
+ "def propagate_orbit(orbit_obj, time_interval, perturbation_function):\n",
+ " \"\"\"\n",
+ " Propagate the orbit over a specified time interval using the Cowell propagator.\n",
+ "\n",
+ " Args:\n",
+ " orbit_obj (Orbit): The initial Orbit object.\n",
+ " time_interval (Quantity): Duration over which to propagate.\n",
+ " perturbation_function (function): Function to compute perturbations.\n",
+ "\n",
+ " Returns:\n",
+ " Orbit: Final orbit after propagation.\n",
+ " \"\"\"\n",
+ " propagated_orbit = orbit_obj.propagate(\n",
+ " time_interval, method=CowellPropagator(rtol=3e-14, f=perturbation_function)\n",
+ " )\n",
+ " gc.collect()\n",
+ " return propagated_orbit"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "uSBQa23XNQiR"
+ },
+ "source": [
+ "# Creating evenly-spaced time intervals\n",
+ "\n",
+ "This function creates epoch objects which are evenly spaced. The number of intervals is user-defined and can affect the how precise our simulation is. A good value should be chosen so that we can have enough precision while still running the program in a reasonable amount of time. Thus, 100 intervals were chosen."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "yTf9vpShNREu"
+ },
+ "outputs": [],
+ "source": [
+ "def create_simulation_epochs(start_epoch, end_epoch, num_intervals):\n",
+ " \"\"\"\n",
+ " Generate an array of epochs spanning from the start to the end epoch.\n",
+ " \"\"\"\n",
+ " # Safely convert the TimeDelta into raw seconds\n",
+ " total_seconds = (end_epoch - start_epoch).to(u.s).value\n",
+ "\n",
+ " # Generate the checkpoints and add them back to the start time\n",
+ " epochs = start_epoch + np.linspace(0, total_seconds, num_intervals + 1) * u.s\n",
+ "\n",
+ " gc.collect()\n",
+ " return epochs"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "NX7MXI5aOFr6"
+ },
+ "source": [
+ "# Extracting data for error analysis\n",
+ "\n",
+ "This function calls the Cowell Propagator for each time interval to obtain the Orbit object for each time step. Then, it extracts the position and velocity vector from each object. These values are converted from GCRS to the ICRS coordinate system for comparison with the JPL Horizons database"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "JcRYLKiZOEcQ"
+ },
+ "outputs": [],
+ "source": [
+ "def simulate_spacecraft_trajectory(orbit_obj, perturbation_function, epoch_array):\n",
+ " \"\"\"\n",
+ " Propagate the spacecraft orbit for each epoch in the provided array and convert the results\n",
+ " to ICRS frame state vectors.\n",
+ "\n",
+ " Args:\n",
+ " orbit_obj (Orbit): Initial orbit.\n",
+ " perturbation_function (function): Function to compute perturbation accelerations.\n",
+ " epoch_array (Quantity): Array of epochs at which to propagate the orbit.\n",
+ "\n",
+ " Returns:\n",
+ " tuple: (r_state, v_state) state arrays in the ICRS frame.\n",
+ " \"\"\"\n",
+ " start_time = time.perf_counter()\n",
+ " propagated_orbits = [\n",
+ " orbit_obj.propagate(ep, method=CowellPropagator(rtol=3e-14, f=perturbation_function))\n",
+ " for ep in epoch_array\n",
+ " ]\n",
+ " end_time = time.perf_counter()\n",
+ " print(\"Propagation time: {:.3f} seconds\".format(end_time - start_time))\n",
+ "\n",
+ " positions, velocities, epochs_list = [], [], []\n",
+ " for orbit_instance in propagated_orbits:\n",
+ " positions.append(orbit_instance.r)\n",
+ " velocities.append(orbit_instance.v)\n",
+ " epochs_list.append(orbit_instance.epoch)\n",
+ "\n",
+ " # Convert lists to numpy arrays\n",
+ " pos_array = np.array(positions)\n",
+ " vel_array = np.array(velocities)\n",
+ " epochs_array = Time(epochs_list)\n",
+ "\n",
+ " # Convert coordinates from GCRS to ICRS for output and analysis.\n",
+ " conversion_start = time.perf_counter()\n",
+ " sky_coordinates = SkyCoord(\n",
+ " CartesianRepresentation(\n",
+ " pos_array[:, 0] * u.km,\n",
+ " pos_array[:, 1] * u.km,\n",
+ " pos_array[:, 2] * u.km,\n",
+ " differentials=CartesianDifferential(\n",
+ " vel_array[:, 0] * (u.km / u.s),\n",
+ " vel_array[:, 1] * (u.km / u.s),\n",
+ " vel_array[:, 2] * (u.km / u.s)\n",
+ " )\n",
+ " ),\n",
+ " frame=GCRS,\n",
+ " obstime=epochs_array\n",
+ " )\n",
+ " icrs_coordinates = sky_coordinates.transform_to(ICRS)\n",
+ " r_state = icrs_coordinates.cartesian.xyz.T.value\n",
+ " v_state = icrs_coordinates.cartesian.differentials[\"s\"].d_xyz.T.value\n",
+ "\n",
+ " # Clean up and free memory.\n",
+ " del positions, velocities, pos_array, vel_array, epochs_list, epochs_array, sky_coordinates, icrs_coordinates\n",
+ " gc.collect()\n",
+ "\n",
+ " conversion_end = time.perf_counter()\n",
+ " print(\"Coordinate conversion time: {:.3f} seconds\".format(conversion_end - conversion_start))\n",
+ " return r_state, v_state"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "WDbOY6Y8PP0E"
+ },
+ "source": [
+ "# Retrieving TESS data from JPL Horizons\n",
+ "\n",
+ "At each time step, this function retrieves the position and velocity vectors from JPL Horizons for TESS"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "bmRkh3r_PQae"
+ },
+ "outputs": [],
+ "source": [
+ "def retrieve_actual_spacecraft_data(epoch_array, num_intervals):\n",
+ " \"\"\"\n",
+ " Obtain reference spacecraft data from Horizons (using TESS as reference) for error comparison.\n",
+ "\n",
+ " Args:\n",
+ " epoch_array (Quantity): Array of simulation epochs.\n",
+ " num_intervals (int): Number of intervals (subdivisions used during simulation).\n",
+ "\n",
+ " Returns:\n",
+ " tuple: (r_actual, v_actual) state arrays from reference data.\n",
+ " \"\"\"\n",
+ " tess_ephem_list = [Ephem.from_horizons(\"TESS\", epochs=ep) for ep in epoch_array]\n",
+ " actual_positions = np.array([ephem.rv()[0].to(u.km).value for ephem in tess_ephem_list])\n",
+ " actual_velocities = np.array([ephem.rv()[1].to(u.km / u.s).value for ephem in tess_ephem_list])\n",
+ "\n",
+ " # Ensure that the arrays have the expected dimensions.\n",
+ " r_actual = actual_positions.reshape(num_intervals + 1, 3)\n",
+ " v_actual = actual_velocities.reshape(num_intervals + 1, 3)\n",
+ " del tess_ephem_list\n",
+ " gc.collect()\n",
+ " return r_actual, v_actual"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "8vTnimfdPo_m"
+ },
+ "source": [
+ "# Computing and plotting fractional error over time\n",
+ "\n",
+ "These functions calculate and plot the fractional error as a function of time for the duration of the simulation. The fractional error for the position and velocity vectors are computed as the difference between the Euclidean norm between the simulation vector and JPL Horizons vector dividided by the JPL Horizons vector. This is done for each time step using vectorized operations with the help of numpy. The results are plotted in two different figures: one for the position and one for the velocity vector"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "FlquX_O1PpeE"
+ },
+ "outputs": [],
+ "source": [
+ "# ---------------------- Error Analysis ----------------------\n",
+ "def calculate_fractional_errors(simulated_positions, reference_positions, simulated_velocities, reference_velocities):\n",
+ " \"\"\"\n",
+ " Calculate the fractional errors between simulated and reference position and velocity vectors.\n",
+ "\n",
+ " Args:\n",
+ " simulated_positions (ndarray): Simulated positions.\n",
+ " reference_positions (ndarray): Reference positions.\n",
+ " simulated_velocities (ndarray): Simulated velocities.\n",
+ " reference_velocities (ndarray): Reference velocities.\n",
+ "\n",
+ " Returns:\n",
+ " tuple: Fractional errors for positions and velocities.\n",
+ " \"\"\"\n",
+ " # Position error calculations\n",
+ " position_error_norm = np.linalg.norm(simulated_positions - reference_positions, axis=1)\n",
+ " reference_position_norm = np.linalg.norm(reference_positions, axis=1)\n",
+ " fractional_error_position = position_error_norm / reference_position_norm\n",
+ "\n",
+ " # Velocity error calculations\n",
+ " velocity_error_norm = np.linalg.norm(simulated_velocities - reference_velocities, axis=1)\n",
+ " reference_velocity_norm = np.linalg.norm(reference_velocities, axis=1)\n",
+ " fractional_error_velocity = velocity_error_norm / reference_velocity_norm\n",
+ "\n",
+ " return fractional_error_position, fractional_error_velocity\n",
+ "\n",
+ "# ---------------------- Plotting Functions ----------------------\n",
+ "def plot_fractional_errors(time_years, position_errors_dict, velocity_errors_dict):\n",
+ " \"\"\"\n",
+ " Plot and save the fractional error comparisons for positions and velocities.\n",
+ "\n",
+ " Args:\n",
+ " time_years (ndarray): Time axis in years.\n",
+ " position_errors_dict (dict): Dictionary of position error arrays keyed by simulation case label.\n",
+ " velocity_errors_dict (dict): Dictionary of velocity error arrays keyed by simulation case label.\n",
+ " \"\"\"\n",
+ " # Plot position fractional errors.\n",
+ " fig_pos, ax_pos = plt.subplots(figsize=(10, 6))\n",
+ " for label, error_array in position_errors_dict.items():\n",
+ " ax_pos.plot(time_years, error_array, 'o-', label=label)\n",
+ " ax_pos.set_title('Fractional Error in Position Magnitude')\n",
+ " ax_pos.set_xlabel('Time Since Start (years)')\n",
+ " ax_pos.set_ylabel('Fractional Error')\n",
+ " ax_pos.legend()\n",
+ " fig_pos.savefig(\"position_magnitude_error.png\", dpi=130)\n",
+ "\n",
+ " # Plot velocity fractional errors.\n",
+ " fig_vel, ax_vel = plt.subplots(figsize=(10, 6))\n",
+ " for label, error_array in velocity_errors_dict.items():\n",
+ " ax_vel.plot(time_years, error_array, 'o-', label=label)\n",
+ " ax_vel.set_title('Fractional Error in Velocity Magnitude')\n",
+ " ax_vel.set_xlabel('Time Since Start (years)')\n",
+ " ax_vel.set_ylabel('Fractional Error')\n",
+ " ax_vel.legend()\n",
+ " fig_vel.savefig(\"velocity_magnitude_error.png\", dpi=130)\n",
+ "\n",
+ " plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "-_FqFJ3-RMho"
+ },
+ "source": [
+ "# Main program function\n",
+ "\n",
+ "This is the main function of the program that handles all other functions. Its main purpose is to update the force_flags dictionary each time with a new force turned \"off\" and plot the results for that corresponding case"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "DyzoIjwGRM_F",
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "# ---------------------- Main Simulation Routine ----------------------\n",
+ "def main_simulation():\n",
+ " # Initialize simulation environment and physical constants.\n",
+ " sim_epoch, moon_interp, sun_interp, jupiter_interp, venus_interp = initialize_environment()\n",
+ " physical_constants = define_physical_constants()\n",
+ "\n",
+ " # Initialize the spacecraft orbit using reference (TESS) data.\n",
+ " initial_orbit = initialize_orbit(sim_epoch)\n",
+ "\n",
+ " # Define the number of intervals for simulation sampling.\n",
+ " NUM_INTERVALS = 100\n",
+ "\n",
+ " # Propagate orbit without any modifications to define final simulation end epoch.\n",
+ " final_orbit = propagate_orbit(\n",
+ " initial_orbit,\n",
+ " SIMULATION_DURATION_DAYS * u.day,\n",
+ " lambda t0, state, k: perturbation_derivatives(\n",
+ " t0, state, k, physical_constants, sun_interp, moon_interp, jupiter_interp,\n",
+ " venus_interp, {\"Moon\": True, \"Sun\": True, \"Jupiter\": True, \"Venus\": True, \"Radiation\": True}\n",
+ " )\n",
+ " )\n",
+ "\n",
+ " # Generate an array of epochs between the start and final orbit epochs.\n",
+ " simulation_epochs = create_simulation_epochs(initial_orbit.epoch, final_orbit.epoch, NUM_INTERVALS)\n",
+ "\n",
+ " # Retrieve reference spacecraft data for error analysis.\n",
+ " reference_positions, reference_velocities = retrieve_actual_spacecraft_data(simulation_epochs, NUM_INTERVALS)\n",
+ "\n",
+ " # Define base force flags and simulation cases (each turning off one perturbation).\n",
+ " base_force_flags = {\n",
+ " \"Moon\": True, \"Sun\": True, \"Jupiter\": True, \"Venus\": True, \"Radiation\": True\n",
+ " }\n",
+ " simulation_cases = {\n",
+ " \"Without Moon\": {**base_force_flags, \"Moon\": False},\n",
+ " \"Without Sun\": {**base_force_flags, \"Sun\": False},\n",
+ " \"Without Jupiter\": {**base_force_flags, \"Jupiter\": False},\n",
+ " \"Without Venus\": {**base_force_flags, \"Venus\": False},\n",
+ " \"Without Radiation\": {**base_force_flags, \"Radiation\": False},\n",
+ " \"All Forces\": base_force_flags.copy()\n",
+ " }\n",
+ "\n",
+ " # Dictionaries to hold the fractional errors for each simulation case.\n",
+ " pos_errors_dict = {}\n",
+ " vel_errors_dict = {}\n",
+ " # Convert epochs difference to years for x-axis plotting.\n",
+ " time_axis_years = (simulation_epochs - simulation_epochs[0]).sec / (3600 * 24 * 365.25)\n",
+ "\n",
+ " # Execute simulation for each force configuration.\n",
+ " for case_label, force_config in simulation_cases.items():\n",
+ " # Define the perturbation function using current force configuration.\n",
+ " perturbation_function = lambda t0, state, k, flags=force_config: perturbation_derivatives(\n",
+ " t0, state, k, physical_constants, sun_interp, moon_interp, jupiter_interp, venus_interp, flags\n",
+ " )\n",
+ " simulated_positions, simulated_velocities = simulate_spacecraft_trajectory(initial_orbit, perturbation_function, simulation_epochs)\n",
+ "\n",
+ " # Compute fractional errors between simulation and reference data.\n",
+ " frac_err_pos, frac_err_vel = calculate_fractional_errors(simulated_positions, reference_positions,\n",
+ " simulated_velocities, reference_velocities)\n",
+ " pos_errors_dict[case_label] = frac_err_pos\n",
+ " vel_errors_dict[case_label] = frac_err_vel\n",
+ "\n",
+ " # Plot and save error comparisons.\n",
+ " plot_fractional_errors(time_axis_years, pos_errors_dict, vel_errors_dict)\n",
+ " gc.collect()\n",
+ "\n",
+ "# Execute the simulation only if this script is run directly.\n",
+ "if __name__ == \"__main__\":\n",
+ " main_simulation()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "D3x_mHN-Q6wr"
+ },
+ "source": [
+ "The following figures are obtained from this code:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "g1xBm_YEQkQ5"
+ },
+ "source": [
+ ""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "tQ0tBVlpRG3S"
+ },
+ "source": [
+ "#"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "# Results and Conclusions\n",
+ "\n",
+ "Both graphs show the need to include the third-body effects from the Moon and the Sun to obtain accurate simulations. The fractional error for the magnitude of both the position and the velocity spikes at a very early stage of the simulation when the gravitational forces from these bodies is ignored. The solar radiation force becomes more relevant at around 6 years after the beginning of the simulation, where it starts to clearly diverge from the \"All forces\" case. In contrast, removing the influence of secondary celestial bodies, such as Jupiter and Venus, on the orbit of the simulated spacecraft shows no clear sign of loss of accuracy. Based on these results, it can be concluded that when simulating the orbit of a spacecraft on High-Earth Orbit (HEO), such as TESS, it is crucial to include the gravitational effects of primary celestial bodies like the Moon and the Sun while the influence from planets like Jupiter and Venus can be removed, saving computational time while not compromising a significant loss of accuracy\n",
+ "\n"
+ ],
+ "metadata": {
+ "id": "rTMJYPxyqxQB"
+ }
+ }
+ ],
+ "metadata": {
+ "colab": {
+ "provenance": [],
+ "authorship_tag": "ABX9TyMUbF1AMZ3NdtaQOd740GBx",
+ "include_colab_link": true
+ },
+ "kernelspec": {
+ "display_name": "Python 3",
+ "name": "python3"
+ },
+ "language_info": {
+ "name": "python"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}
\ No newline at end of file
diff --git a/corner_code.ipynb b/corner_code.ipynb
new file mode 100644
index 0000000..fc6031d
--- /dev/null
+++ b/corner_code.ipynb
@@ -0,0 +1,281 @@
+{
+ "nbformat": 4,
+ "nbformat_minor": 0,
+ "metadata": {
+ "colab": {
+ "provenance": [],
+ "authorship_tag": "ABX9TyNHjjBoncuSIBPpLueRZOte",
+ "include_colab_link": true
+ },
+ "kernelspec": {
+ "name": "python3",
+ "display_name": "Python 3"
+ },
+ "language_info": {
+ "name": "python"
+ }
+ },
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "view-in-github",
+ "colab_type": "text"
+ },
+ "source": [
+ ""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "# Project Overview\n",
+ "\n",
+ "This code analyzes how different launch parameters (like launch latitude, exhaust velocity, etc.) affect the required delta-V and propellant mass to insert a payload in a geosynchronous orbit with a desired inclination. It creates a corner plots to understand the dependency of the delta-V and fuel mass on these launch variables"
+ ],
+ "metadata": {
+ "id": "zCNyd43A3Y-6"
+ }
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true,
+ "id": "FCGMNdHQ2wJz"
+ },
+ "outputs": [],
+ "source": [
+ "#Install necessary packages for Google Colab. Note: hapsira is the modern version of poliastro to handle Python 3.12\n",
+ "!pip install hapsira \"astropy<7.0\" corner"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "# Code Imports\n",
+ "This code utilizes several Python libraries for different functionalities:\n",
+ "\n",
+ "1. poliastro & astropy: for orbital state representation\n",
+ "\n",
+ "2. multiprocessing: to handle multiple orbital simulations simultaneously to reduce computational time\n",
+ "\n",
+ "3. corner & matplotlib: to plot results\n",
+ "\n",
+ "4. random: to randomize some launch parameters\n",
+ "\n",
+ "5. time: to keep track of the computational time\n",
+ "\n",
+ "6. numpy: to make complex vector calculations\n",
+ "\n",
+ "7. simulate_orbit: function that simulates the insertion of a payload into GEO orbit and computes the total delta-V and propellant mass needed given specific launch parameters (from GEO_orbital_simulation tutorial)"
+ ],
+ "metadata": {
+ "id": "WmVwZAuv4Iba"
+ }
+ },
+ {
+ "cell_type": "code",
+ "source": [
+ "import numpy as np\n",
+ "import time\n",
+ "import random\n",
+ "from itertools import product\n",
+ "import multiprocessing as mp\n",
+ "import matplotlib as mpl\n",
+ "from matplotlib import pyplot as plt\n",
+ "import corner\n",
+ "\n",
+ "from astropy import units as u\n",
+ "from hapsira.bodies import Earth\n",
+ "from hapsira.twobody import Orbit\n",
+ "from astropy.coordinates import CartesianRepresentation\n",
+ "\n",
+ "# Import the simulation function from the simulation module.\n",
+ "# Ensure that the function 'simulate' is available in the 'simulation.py' module in your environment.\n",
+ "from simulation import simulate_orbit"
+ ],
+ "metadata": {
+ "collapsed": true,
+ "id": "RCxEH_aK4Onk"
+ },
+ "execution_count": null,
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "# Defining launch parameters\n",
+ "\n",
+ "All the launch parameters that will be analyzed are initialized here. Discrete values are chosen for all quantities and are randomized within a range except the final inclination and exhaust velocities of the two stages of the rocket."
+ ],
+ "metadata": {
+ "id": "SJ0Gu2JM7r9d"
+ }
+ },
+ {
+ "cell_type": "code",
+ "source": [
+ "# ===============================================================\n",
+ "# Define Parameter Lists for Simulation\n",
+ "# ===============================================================\n",
+ "# Fixed list of inclination values (in degrees)\n",
+ "inclination_values = [5, 13, 26, 34, 46, 62]\n",
+ "\n",
+ "# Fixed list of stage 1 exhaust velocities (in m/s)\n",
+ "exhaust_velocity_stage1 = [2800, 2925, 3050, 3175, 3300]\n",
+ "\n",
+ "# Fixed list of stage 2 exhaust velocities (in m/s)\n",
+ "exhaust_velocity_stage2 = [4200, 4300, 4400, 4500, 4600]\n",
+ "\n",
+ "# Generate randomized list for payload exhaust velocity (propulsive exhaust) in m/s\n",
+ "payload_exhaust_lin = np.linspace(20000, 30000, 5)\n",
+ "payload_exhaust_values = [\n",
+ " np.clip(val + np.random.uniform(-1000, 1000), 20000, 30000)\n",
+ " for val in payload_exhaust_lin\n",
+ "]\n",
+ "\n",
+ "# Generate randomized list for LEO altitude values (in km)\n",
+ "leo_altitude_lin = np.linspace(200, 800, 10)\n",
+ "leo_altitude_values = [\n",
+ " np.clip(val + np.random.uniform(-10, 10), 200, 800)\n",
+ " for val in leo_altitude_lin\n",
+ "]\n",
+ "\n",
+ "# Generate randomized list for payload mass values (in kg)\n",
+ "payload_mass_lin = np.linspace(200, 300, 5)\n",
+ "payload_mass_values = [\n",
+ " np.clip(val + np.random.uniform(-10, 10), 200, 300)\n",
+ " for val in payload_mass_lin\n",
+ "]\n",
+ "\n",
+ "# Generate randomized list for inclination change percentage (in %)\n",
+ "inclination_change_lin = np.linspace(0, 100, 10)\n",
+ "inclination_change_values = [\n",
+ " np.clip(val + np.random.uniform(-5, 5), 0, 100)\n",
+ " for val in inclination_change_lin\n",
+ "]"
+ ],
+ "metadata": {
+ "id": "XAR7Ntb97sca"
+ },
+ "execution_count": 7,
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "# Main code\n",
+ "\n",
+ "Here, we combine all the launch parameters into a tuple and pass them over to the simulation_orbit() function. The multiprocessing tool is utilized to run multiple simulations with different parameters at the same time to run the code faster. Lastly, we create a corner plot using the corner function and adjust font size, dpi of the picture, and other variabes to create a polished figure with the help of matplotlib"
+ ],
+ "metadata": {
+ "id": "158lPIJX8QTI"
+ }
+ },
+ {
+ "cell_type": "code",
+ "source": [
+ "# ===============================================================\n",
+ "# Main Entry Point for Simulation and Plotting\n",
+ "# ===============================================================\n",
+ "if __name__ == '__main__':\n",
+ " # Mark start time of simulation\n",
+ " start_time = time.perf_counter()\n",
+ "\n",
+ " # Create a list of all possible parameter combinations as tuples:\n",
+ " # (inclination, LEO altitude, stage 1 exhaust, stage 2 exhaust, payload exhaust, payload mass, inclination change)\n",
+ " parameter_combinations = list(product(\n",
+ " inclination_values,\n",
+ " leo_altitude_values,\n",
+ " exhaust_velocity_stage1,\n",
+ " exhaust_velocity_stage2,\n",
+ " payload_exhaust_values,\n",
+ " payload_mass_values,\n",
+ " inclination_change_values\n",
+ " ))\n",
+ "\n",
+ " # Use a multiprocessing pool to parallelize the simulation.\n",
+ " with mp.Pool(processes=mp.cpu_count()) as pool:\n",
+ " simulation_results = pool.map(simulate_orbit, parameter_combinations)\n",
+ "\n",
+ " # Convert the results to a NumPy array for further processing/plotting.\n",
+ " results_array = np.array(simulation_results)\n",
+ "\n",
+ " # Set global font size for plotting with matplotlib\n",
+ " mpl.rcParams.update({'font.size': 115})\n",
+ "\n",
+ " # Define axis labels for the corner plot\n",
+ " plot_labels = [\n",
+ " \"Initial Inclination (deg)\",\n",
+ " \"Inclination Change (%)\",\n",
+ " \"LEO Altitude (km)\",\n",
+ " \"Stage 1 EV (m/s)\",\n",
+ " \"Stage 2 EV (m/s)\",\n",
+ " \"Payload EV (m/s)\",\n",
+ " \"Payload Mass (kg)\",\n",
+ " \"Propellant Mass (kg)\",\n",
+ " \"Total ΔV (m/s)\"\n",
+ " ]\n",
+ "\n",
+ " # Generate the corner plot with the simulation results\n",
+ " corner_fig = corner.corner(\n",
+ " results_array,\n",
+ " labels=plot_labels,\n",
+ " label_kwargs={\"fontsize\": 120},\n",
+ " show_titles=True,\n",
+ " title_fmt=\".2f\",\n",
+ " bins=60,\n",
+ " smooth=1.0,\n",
+ " fill_contours=True,\n",
+ " plot_contours=True,\n",
+ " max_n_ticks=4,\n",
+ " hist2d_kwargs={\n",
+ " 'plot_datapoints': True,\n",
+ " 'alpha': 0.15,\n",
+ " 'plot_density': False\n",
+ " },\n",
+ " scatter_kwargs={\n",
+ " 'alpha': 0.05,\n",
+ " 's': 3\n",
+ " }\n",
+ " )\n",
+ "\n",
+ " # Set figure size and layout adjustments, then save and display the plot.\n",
+ " corner_fig.set_size_inches(220, 220)\n",
+ " corner_fig.tight_layout()\n",
+ " plt.savefig(\"Corner.pdf\", dpi=130, bbox_inches='tight')\n",
+ "\n",
+ " # Mark end time and output the duration of the simulation\n",
+ " end_time = time.perf_counter()\n",
+ " print(\"Total simulation and plotting time: \", end_time - start_time, \"seconds\")\n",
+ " plt.show()"
+ ],
+ "metadata": {
+ "id": "L0pRoWks8Qrc"
+ },
+ "execution_count": null,
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "Given the size of this figure, google colab crashes out when it tries to display this picture. Thus, a link is provided with direct access to this figure. Click here to view a PDF version of the corner plot: https://drive.google.com/drive/folders/1mz0qIrBFHchnG6LQ4o484kGEIt2xbMU1?usp=sharing"
+ ],
+ "metadata": {
+ "id": "rkm5DPtP9Ugu"
+ }
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "# Results and Conclusions\n",
+ "\n",
+ "As shown in the corner plot above, the three main factors to consider to otpimize the payload mass are the launch site latitude, the exhaust velocity of the stage 2 of the rocket, and the mass of the payload. Similarly, the percent change in final inclination and initial launch site latitude are the factors that most heavily influence the total required delta-V. These results highlight the need for a strategic design of the payload to reduce the fuel mass, as well as selecting high-performance engines, especially for the stage 2 of the rocket. Additionally, the initial and final inclinations of the orbit must be close in value if minimizing the delta-V is a priority."
+ ],
+ "metadata": {
+ "id": "v7KAZ3SkAgGl"
+ }
+ }
+ ]
+}
\ No newline at end of file