anise-py/tutorials/Tutorial 03 - Defining and working with the Orbit structure.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e27f64b8-2625-40a7-a6b9-c4b9a9281530",
   "metadata": {},
   "source": [
    "# ANISE\n",
    "\n",
    "ANISE is a modern rewrite of NAIF SPICE, written in Rust and providing interfaces to other languages include Python.\n",
    "\n",
    "Evidently, this tutorial applies to the Python usage of ANISE.\n",
    "\n",
    "## Goal\n",
    "By the end of this tutorial, you should know how to initialize an `Orbit` structure, how to retrieve any of the available orbital elements, and compute the future state of that orbit in the future assuming two body propagation (i.e. _no perturbation of any kind_).\n",
    "\n",
    "For comprehensive and high-fidelity astrodynamics, orbit estimation, and other spaceflight dynamics functionalities essential in real-world missions, consider using [Nyx](https://github.com/nyx-space/nyx/). Nyx has been the choice for several cislunar missions and lunar landers, and, like ANISE, is also open-source and free to use.\n",
    "\n",
    "Let's start by installing ANISE: `pip install anise`\n",
    "\n",
    "## Load the latest Almanac\n",
    "\n",
    "As seen the tutorial #2, we'll just download the latest data using the MetaAlmanac."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "1787d97b-d706-4a29-8666-4addf8a4d95d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Almanac: #SPK = 1\t#BPC = 2\tPlanetaryData with 49 ID mappings and 0 name mappings (@0x5587ec685530)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from anise import MetaAlmanac\n",
    "\n",
    "almanac = MetaAlmanac(\"../../data/latest.dhall\").process()\n",
    "\n",
    "almanac"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d95d9a6-c231-4326-b5b7-c7a19a5ff617",
   "metadata": {},
   "source": [
    "## Orbit structure\n",
    "\n",
    "The `Orbit` structure is an alias for `CartesianState`, which is only available under that name in Rust. In fact, `Orbit` will _always_ store the information in its Cartesian form because it's a non-singular representation of an orbit, regardless of whether it's elliptical, hyperbolic, equatorial and circular, etc.\n",
    "\n",
    "An orbit is defined by its position, velocity, an epoch, and a frame.\n",
    "\n",
    "As we saw in tutorial #1, we can import the frames as they're defined with their ephemeris and orientation IDs. However, recall that ANISE does not provide _any_ default values for the gravitational parameter or shape of the geoid in these frames and they need to be loaded from the Almanac."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f4de9b9-c2fb-4774-837d-e70da5afc84a",
   "metadata": {},
   "source": [
    "Let's grab the frame information from the almanac."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "cbe97591-9aec-48cb-b63f-1f608b10b405",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Earth J2000\n",
      "Earth J2000 (μ = 398600.435436096 km^3/s^2, eq. radius = 6378.14 km, polar radius = 6356.75 km, f = 0.0033536422844278)\n"
     ]
    }
   ],
   "source": [
    "from anise.astro.constants import Frames\n",
    "eme2k = almanac.frame_info(Frames.EME2000)\n",
    "print(Frames.EME2000)\n",
    "# Only the loaded frame has the gravitational parameter and shape information\n",
    "print(eme2k)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c20bfa85-da31-4dda-a9b1-881bf5c0186c",
   "metadata": {},
   "source": [
    "Orbits can be initialized in several ways. The simplest is to initialize it from its cartesian data, provided in kilometers and kilometers per second."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3c457886-95ed-40e4-8b49-9ca8e08b1e7f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[0;31mInit signature:\u001b[0m \u001b[0mOrbit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx_km\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my_km\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mz_km\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvx_km_s\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvy_km_s\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvz_km_s\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mepoch\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mframe\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;31mDocstring:\u001b[0m     \n",
       "Defines a Cartesian state in a given frame at a given epoch in a given time scale. Radius data is expressed in kilometers. Velocity data is expressed in kilometers per second.\n",
       "Regardless of the constructor used, this struct stores all the state information in Cartesian coordinates as these are always non singular.\n",
       "\n",
       "Unless noted otherwise, algorithms are from GMAT 2016a [StateConversionUtil.cpp](https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/src/base/util/StateConversionUtil.cpp).\n",
       "\u001b[0;31mType:\u001b[0m           type\n",
       "\u001b[0;31mSubclasses:\u001b[0m     "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from anise.astro import Orbit\n",
    "Orbit?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b612912a-9f93-472a-957a-ecd4d0e70bb9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Earth J2000] 2010-12-21T02:01:57.816398952 UTC\tposition = [5946.673548, 1656.154606, 2259.012130] km\tvelocity = [-3.098683, 4.579534, 6.246542] km/s\n"
     ]
    }
   ],
   "source": [
    "from anise.time import Epoch\n",
    "\n",
    "epoch = Epoch(\"2010-12-21 02:03:04 TDB\") # Let's create an epoch in the TDB time system\n",
    "\n",
    "# And initialize the orbit from an arbitrary position and velocity.\n",
    "\n",
    "from_cartesian = Orbit.from_cartesian(\n",
    "            5_946.673548288958,\n",
    "            1_656.154606023661,\n",
    "            2_259.012129598249,\n",
    "            -3.098683050943824,\n",
    "            4.579534132135011,\n",
    "            6.246541551539432,\n",
    "            epoch,\n",
    "            eme2k\n",
    "        )\n",
    "print(from_cartesian)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5431a376-8e85-4573-9b1d-9514ee24c78e",
   "metadata": {},
   "source": [
    "Oftentimes, we know the scientific mission of a spacecraft, so we initialize it from its Keplerian parameters and not its Cartesian state.\n",
    "\n",
    "All of the orbit initializers start with `from_...`, making it easy to spot if you have an autocompletion enabled (you should!)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "fc94db3b-9ce3-4113-a96c-5e3c4318bdfb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[0;31mSignature:\u001b[0m \u001b[0mOrbit\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfrom_keplerian\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mecc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mraan\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maop\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mta\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mepoch\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mframe\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;31mDocstring:\u001b[0m\n",
       "Creates a new Orbit around the provided Celestial or Geoid frame from the Keplerian orbital elements.\n",
       "\n",
       "**Units:** km, none, degrees, degrees, degrees, degrees\n",
       "\n",
       "NOTE: The state is defined in Cartesian coordinates as they are non-singular. This causes rounding\n",
       "errors when creating a state from its Keplerian orbital elements (cf. the state tests).\n",
       "One should expect these errors to be on the order of 1e-12.\n",
       "\u001b[0;31mType:\u001b[0m      builtin_function_or_method"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Orbit.from_keplerian?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4542fa57-7473-4f3f-ba67-d381a94f7c7e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Earth J2000] 2010-12-21T02:01:57.816398952 UTC\tposition = [5946.673549, 1656.154605, 2259.012130] km\tvelocity = [-3.098683, 4.579534, 6.246542] km/s\n"
     ]
    }
   ],
   "source": [
    "from_keplerian = Orbit.from_keplerian(\n",
    "    7_712.186_117_895_041,\n",
    "    0.158_999_999_999_999_95,\n",
    "    53.75369,\n",
    "    1.998_632_864_211_17e-5,\n",
    "    359.787_880_000_004,\n",
    "    25.434_003_407_751_188,\n",
    "    epoch,\n",
    "    eme2k\n",
    ")\n",
    "print(from_keplerian)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef0d35eb-d1a3-4a4b-b35f-eebd4919ae7f",
   "metadata": {},
   "source": [
    "### Exercise #1\n",
    "\n",
    "We won't go through all of the initializers, so use your autocompletion to try the following ones out:\n",
    "\n",
    "+ `from_keplerian_mean_anomaly`\n",
    "+ `from_keplerian_apsis_radii`\n",
    "+ `from_keplerian_altitude`\n",
    "+ `from_latlongalt`\n",
    "\n",
    "**Objective:** Create and initialize an Orbit structure with given orbital elements.\n",
    "\n",
    "**Tasks:**\n",
    " 1. Input specific orbital elements to define an orbit, including hyperbolic orbits (eccentricity greater than 1)\n",
    " 2. Display the initialized orbit’s parameters and notice how initializing from the Keplerian parameters won't return _exactly_ the same Keplerian data, as it's converted back and forth from its Keplerian representation to its Cartesian representation.\n",
    "\n",
    "**Learning Goal:** Understand how to define an orbit in ANISE and interpret its parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f0bbad6-ea5c-464f-8907-22acbf569d74",
   "metadata": {},
   "source": [
    "## Fields of the orbit class\n",
    "\n",
    "The orbit class provides direct access to a large range of orbital parameters, including the semi major axis, the true longitude, the orbital period, and so on. All of these accessors include the parameter along with the unit in which it is returned. Let's look at a few more common ones."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "202b81e4-9a6e-45a6-a4f7-aa99c90cd2c1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1334.0462758738768,\n",
       " 25.221883407752955,\n",
       " -25.842246360350657,\n",
       " 107.80853340089652)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from_cartesian.sma_altitude_km(), from_cartesian.aol_deg(), from_cartesian.energy_km2_s2(), from_cartesian.periapsis_altitude_km()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "1377c6ce-ad67-4a60-9f2a-80b2d50b717d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 h 52 min 20 s 269 ms 114 μs 912 ns\n",
      "1 h 52 min 20 s\n",
      "6740.269114912\n",
      "21 ms 926 μs 330 ns\n"
     ]
    }
   ],
   "source": [
    "# The orbital period (and the light time) is returned as a Hifitime Duration, enabling trivial computations\n",
    "from anise.time import Unit\n",
    "\n",
    "print(from_keplerian.period())\n",
    "print(from_keplerian.period().round(Unit.Second*1))\n",
    "print(from_keplerian.period().to_seconds())\n",
    "\n",
    "print(from_keplerian.light_time())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da8c8062-47f8-45d1-8b7e-3a112f7f0ee7",
   "metadata": {},
   "source": [
    "**Importantly**, all of the orbital parameter computations require the frame data to be initialized, otherwise an exception will be raised specifying what data was not available. In almost all cases, this happens because you use the constant frame definition without fetching the frame information from the Almanac, or your Almanac does not include any planetary constants.\n",
    "\n",
    "The Orbit structure supports the equality operation, and it'll check that the epoch, frame, and Cartesian state match to within the centimeter level in position and centimeter per second level in velocity.\n",
    "\n",
    "If that precision is different than the one you wish, you can use the `eq_within` function providing a tolerance in kilometers and kilometers per second."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "c915cdc7-1932-438b-87a9-c7c66591822d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from_keplerian == from_cartesian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ca0931d5-4df0-46e7-97f3-29d5f180d2a2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[0;31mSignature:\u001b[0m \u001b[0mfrom_keplerian\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meq_within\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mother\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mradial_tol_km\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvelocity_tol_km_s\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;31mDocstring:\u001b[0m Returns whether this orbit and another are equal within the specified radial and velocity absolute tolerances\n",
       "\u001b[0;31mType:\u001b[0m      builtin_function_or_method"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from_keplerian.eq_within?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1546978b-d1cd-41bd-bab3-23bb8fe27f0e",
   "metadata": {},
   "source": [
    "### Exercise #2\n",
    "\n",
    "**Objective:** Convert between different sets of orbital elements.\n",
    "\n",
    "**Tasks:**\n",
    " 1. Using the `add_...` methods (e.g. `add_sma`), see how the Cartesian elements change (accessible as getters like `x_km`).\n",
    " 2. Analyze how changes in one element affect others.\n",
    "\n",
    "**Learning Goal:** Gain proficiency in working with various orbital element formats."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "536338b1-392a-4884-8c14-c62c14bbfdf4",
   "metadata": {},
   "source": [
    "## Two body propagation\n",
    "\n",
    "### Important disclaimer\n",
    "\n",
    "Two body propagation is a _tool_ and should not be mistaken for a complete representation of real-world scenarios. In actual mission design, the exclusive use of two-body propagation is _physically incorrect and the wrong approach_. This is because a spacecraft or celestial body is invariably influenced by a multitude of factors. These include gravitational forces from other celestial bodies, as well as subtler forces like solar radiation pressure, atmospheric drag, and gravity fields (as, e.g., spherical harmonic modeling), which account for the uneven distribution of a celestial object's mass. **ANISE will only ever support two-body propagation. For high-fidelity astrodynamics, including orbit estimation, the recommended resource is Nyx. Nyx has been integral to several high-profile lunar missions, offering the necessary models for accurate space mission analysis, planning, and operations.**\n",
    "\n",
    "### Usage\n",
    "\n",
    "Using the two body propagator in ANISE is super straight-forward: simply call `at_epoch` on an `Orbit` instance. This uses the current mean anomaly of this instance to compute its future mean anomaly using only the time difference between the current epoch and the future epoch.\n",
    "\n",
    "Of course, this converts from the internal Cartesian representation to its Keplerian representation, and back again, so some level of error in the Keplerian orbital elements will be seen, but it should be minimal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "98d91960-34cb-4856-986b-0bb9c0116581",
   "metadata": {},
   "outputs": [],
   "source": [
    "current_sma_km = from_keplerian.sma_km()\n",
    "current_true_anomaly_deg = from_keplerian.ta_deg()\n",
    "current_mean_anomaly_deg = from_keplerian.ma_deg()\n",
    "future_state = from_keplerian.at_epoch(from_keplerian.epoch + Unit.Day*1.7565)\n",
    "# Let's see how much these two Keplerian orbital elements have changed\n",
    "future_sma_km = future_state.sma_km()\n",
    "future_true_anomaly_deg = future_state.ta_deg()\n",
    "future_mean_anomaly_deg = future_state.ma_deg()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "0e5f6e9f-2349-4b11-86f2-66e0aa1f0349",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error in SMA:\t9.094947017729282e-13 km\n",
      "Propagated true anomaly:\t172.3281610410407 deg\n",
      "Propagated mean anomaly:\t185.63718609798167 deg\n"
     ]
    }
   ],
   "source": [
    "print(f\"Error in SMA:\\t{abs(current_sma_km - future_sma_km)} km\")\n",
    "print(f\"Propagated true anomaly:\\t{abs(current_true_anomaly_deg - future_true_anomaly_deg)} deg\")\n",
    "print(f\"Propagated mean anomaly:\\t{abs(current_mean_anomaly_deg - future_mean_anomaly_deg)} deg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b25ddb5-1d23-4b8a-94f7-fd81f667b3f1",
   "metadata": {},
   "source": [
    "### Exercise #3\n",
    "\n",
    "All Almanac translation/transform queries return an instance of `Orbit`, populated with the relevant frame information from the Almanac, so gaining familiarity with this structure is important.\n",
    "\n",
    "**Objective:** Simulate the orbit of a body over time using two-body propagation.\n",
    "\n",
    "**Tasks:**\n",
    " 1. Fetch the state of a celestial object, like the Moon, at the epoch of your choice.\n",
    " 3. Use the Orbit instance's `at_epoch` to propagate the position of that object forward in time, assuming it was subjected only to two body dynamics.\n",
    " 4. Fetch the state of that same object from the Almanac at that other time, and look at the difference in Cartesian state and orbital elements.\n",
    "\n",
    "**Learning Goal:** Understand the basics of orbital propagation and its limitations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "002aa86f-9ae2-484a-a6cb-441a9bd26b40",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "language": "python",
   "name": ".venv"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}