anise-py/tutorials/Tutorial 01 - Querying SPK files.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c6f0bc55-c0be-46c6-bb54-32df56aeba1f",
   "metadata": {},
   "source": [
    "# ANISE\n",
    "\n",
    "ANISE is a modern rewrite of NAIF SPICE, written in Rust and providing interfaces to other languages including 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 be able to know how to load local files and query the Almanac for the position and velocity of different celestial objects at different times, expressed in different time systems.\n",
    "\n",
    "Let's start by installing ANISE: `pip install anise`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01c54376-a0cb-47b1-bdcd-da9a85242602",
   "metadata": {},
   "source": [
    "## Introduction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2ca6dd47-507d-4e9a-a767-e46d2cc5a836",
   "metadata": {},
   "outputs": [],
   "source": [
    "import anise"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc4ff306-68be-407f-9571-0ce23db65e27",
   "metadata": {},
   "source": [
    "ANISE allows one to query SPICE files, like those created by JPL, for position and rotation information. A contrario to SPICE, ANISE is fully thread safe. To support this, ANISE has an object called the `Almanac` which stores all of the loaded files, whereas SPICE would have a single thread-unsafe shared context.\n",
    "\n",
    "Let's import that class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "1502349d-444c-4512-90ae-01bcae54a573",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[0;31mInit signature:\u001b[0m \u001b[0mAlmanac\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpath\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",
       "An Almanac contains all of the loaded SPICE and ANISE data.\n",
       "\n",
       "# Limitations\n",
       "The stack space required depends on the maximum number of each type that can be loaded.\n",
       "\u001b[0;31mFile:\u001b[0m           ~/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages/anise/__init__.py\n",
       "\u001b[0;31mType:\u001b[0m           type\n",
       "\u001b[0;31mSubclasses:\u001b[0m     "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from anise import Almanac\n",
    "Almanac?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f215f0f3-9bd7-4c45-b097-aefa55ff523a",
   "metadata": {},
   "source": [
    "As much as possible, the functions and classes in ANISE are documented. In typical Python, you can do `help(Almanac)` and in Jupyter notebook, you can simply enter the name of what you're seeking help on followed by a question mark, and execute the cell (Ctrl+Enter).\n",
    "\n",
    "Let's load on the files provided in the ANISE repo. Note that ANISE supports automatically downloading remote files as well, but that is the topic of another tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "4c91ffd3-ebe3-4a52-8617-5a8b7293ecfe",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Almanac: #SPK = 1\t#BPC = 0\n"
     ]
    }
   ],
   "source": [
    "almanac = Almanac(\"../../data/de440s.bsp\")\n",
    "print(almanac)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3759ef17-cd5f-435e-869d-a68b21eedeee",
   "metadata": {},
   "source": [
    "You may also use the `describe` function to see the details of what is loaded, optionally showing only some of the loaded data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "67ee44ce-d5b6-4c86-ba49-7ea89254a9b3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[0;31mSignature:\u001b[0m\n",
       "\u001b[0malmanac\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdescribe\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mspk\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mbpc\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mplanetary\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mtime_scale\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mround_time\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\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",
       "Pretty prints the description of this Almanac, showing everything by default. Default time scale is TDB.\n",
       "If any parameter is set to true, then nothing other than that will be printed.\n",
       "\u001b[0;31mType:\u001b[0m      builtin_function_or_method"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "almanac.describe?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d120774f-09a5-4e38-bc96-f632e6fd746b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=== SPK #0 ===\n",
      "┌────────────────┬─────────────────────────────┬───────────────────────────────┬───────────────────────────────────┬───────────────────────────────────┬─────────────┬────────────────────┐\n",
      "│ Name           │ Target                      │ Center                        │ Start epoch                       │ End epoch                         │ Duration    │ Interpolation kind │\n",
      "├────────────────┼─────────────────────────────┼───────────────────────────────┼───────────────────────────────────┼───────────────────────────────────┼─────────────┼────────────────────┤\n",
      "│ DE-0440LE-0440 │ Mercury J2000               │ Solar System Barycenter J2000 │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet  │\n",
      "│ DE-0440LE-0440 │ Venus J2000                 │ Solar System Barycenter J2000 │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet  │\n",
      "│ DE-0440LE-0440 │ Earth-Moon Barycenter J2000 │ Solar System Barycenter J2000 │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet  │\n",
      "│ DE-0440LE-0440 │ Mars Barycenter J2000       │ Solar System Barycenter J2000 │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet  │\n",
      "│ DE-0440LE-0440 │ Jupiter Barycenter J2000    │ Solar System Barycenter J2000 │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet  │\n",
      "│ DE-0440LE-0440 │ Saturn Barycenter J2000     │ Solar System Barycenter J2000 │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet  │\n",
      "│ DE-0440LE-0440 │ Uranus Barycenter J2000     │ Solar System Barycenter J2000 │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet  │\n",
      "│ DE-0440LE-0440 │ Neptune Barycenter J2000    │ Solar System Barycenter J2000 │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet  │\n",
      "│ DE-0440LE-0440 │ Pluto Barycenter J2000      │ Solar System Barycenter J2000 │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet  │\n",
      "│ DE-0440LE-0440 │ Sun J2000                   │ Solar System Barycenter J2000 │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet  │\n",
      "│ DE-0440LE-0440 │ Moon J2000                  │ Earth-Moon Barycenter J2000   │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet  │\n",
      "│ DE-0440LE-0440 │ Earth J2000                 │ Earth-Moon Barycenter J2000   │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet  │\n",
      "│ DE-0440LE-0440 │ body 199 J2000              │ Mercury J2000                 │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet  │\n",
      "│ DE-0440LE-0440 │ body 299 J2000              │ Venus J2000                   │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet  │\n",
      "└────────────────┴─────────────────────────────┴───────────────────────────────┴───────────────────────────────────┴───────────────────────────────────┴─────────────┴────────────────────┘\n"
     ]
    }
   ],
   "source": [
    "almanac.describe(bpc=True) # This will print nothing because no BPC is loaded\n",
    "almanac.describe(spk=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "256e5365-1f6d-4c98-b6cb-f63a8ad18c51",
   "metadata": {},
   "source": [
    "As of version 0.2.2, in Python, one can only perform transformations, and simple translations, but not extract the rotation information like the DCM (it is only available in Rust for now). What this means in practice is that one can query for the state of the Moon in the Earth body fixed frame (IAU Earth for example), and ANISE will perform all of the required transformations, but you cannot only query for the rotation parameters.\n",
    "\n",
    "**To query something, we need several pieces of information:**\n",
    "\n",
    "1. The epoch (or \"datetime\") at which to compute the state;\n",
    "2. The frames we want the observer and target to be in;\n",
    "3. The aberration corrections desired.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f9ee8b4-1d39-4d5c-a713-d8dd6d7f2eb1",
   "metadata": {},
   "source": [
    "## Time management\n",
    "For time management, ANISE uses [Hifitime](https://nyxspace.com/hifitime/intro/), arguably the most precise time management library in the world. Let's create an epoch by importing the hifitime classes, which are re-exported from ANISE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "c4c1d46e-4701-4ccb-9a27-0832265f973e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2010-12-21T03:04:05 UTC\n",
      "2010-12-21T03:03:31 UTC\n",
      "-34 s\n"
     ]
    }
   ],
   "source": [
    "from anise.time import Epoch\n",
    "\n",
    "epoch_utc = Epoch(\"2010-12-21T03:04:05\")\n",
    "print(epoch_utc)\n",
    "\n",
    "epoch_tai = Epoch(\"2010-12-21T03:04:05 TAI\")\n",
    "print(epoch_tai)\n",
    "print(epoch_tai.timedelta(epoch_utc)) # This will print the number of leap seconds between UTC and TAI."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a31ba37-1173-4e07-a5af-785799e1c08e",
   "metadata": {},
   "source": [
    "We created an Epoch in the TAI and in the UTC time scales. UTC is a glitchy time scale, which is never used in astrodynamics or astronomy. If the notion of a time scale is new to you, that's fine, just please read the Hifitime link above which explains the concept of time scales and why they are so important.\n",
    "\n",
    "So if UTC is not used in astrodynamics, what is? Well that would be the GPST time scale for low Earth orbiters and the TDB time scale for most things out of LEO."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "5526b332-d48a-469c-a24d-eb59c85dbd90",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2010-12-21T03:04:05 UTC is 11307 days 3 h 4 min 20 s after the GPS reference time of 1980-01-06T00:00:00 UTC\n",
      "2010-12-21T03:04:05 UTC is 4006 days 15 h 5 min 11 s 183 ms 602 μs 262 ns after the Dynamical Barycentric Time reference time of 2000-01-01T11:58:55.816072704 UTC\n"
     ]
    }
   ],
   "source": [
    "epoch = epoch_utc\n",
    "\n",
    "print(f\"{epoch} is {epoch.to_gpst_duration()} after the GPS reference time of {Epoch.init_from_gpst_nanoseconds(0)}\")\n",
    "\n",
    "print(f\"{epoch} is {epoch.to_tdb_duration()} after the Dynamical Barycentric Time reference time of {Epoch.init_from_tdb_seconds(0.0)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "725f779a-9a42-4dcd-b65e-353d824bde87",
   "metadata": {},
   "source": [
    "^^^ That is one of the most useful features of Hifitime: seamless conversion between time scales, all with nanosecond precision, and for ~65,000 centuries around 01 January 1900.\n",
    "\n",
    "At this point, we've got one epoch. More analysis requires us to check stuff at different times.\n",
    "\n",
    "To do so, Hifitime provides arithmetics on epochs. Let's try that out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "2ea6b1b1-ce75-45c3-9324-133e05190762",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2010-12-22T16:40:33.800000020 UTC\n",
      "2010-12-22T16:40:33.800000020 UTC\n"
     ]
    }
   ],
   "source": [
    "from anise.time import Duration, Unit\n",
    "\n",
    "print(epoch + Duration(\"1.567 days 20 ns\"))\n",
    "print(epoch + Unit.Day*1.567 + Unit.Nanosecond*20)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0cc14e3f-05bb-47b6-a565-661fbf86f5d3",
   "metadata": {},
   "source": [
    "Durations can be initialized using either their string representation using standard unit abbreviations.\n",
    "\n",
    "If we know exactly that we want to sample a certain number of times between two epoch, we can either write a for loop and add a duration, or we can use the `TimeSeries` class of Hifitime."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "5e234050-3153-48b3-bae1-005c215df7e9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2010-12-21T03:04:05 UTC\n",
      "2010-12-21T06:46:05 UTC\n",
      "2010-12-21T10:28:05 UTC\n",
      "2010-12-21T14:10:05 UTC\n",
      "2010-12-21T17:52:05 UTC\n",
      "2010-12-21T21:34:05 UTC\n",
      "2010-12-22T01:16:05 UTC\n",
      "2010-12-22T04:58:05 UTC\n",
      "2010-12-22T08:40:05 UTC\n",
      "2010-12-22T12:22:05 UTC\n",
      "2010-12-22T16:04:05 UTC\n",
      "2010-12-22T19:46:05 UTC\n",
      "2010-12-22T23:28:05 UTC\n",
      "2010-12-23T03:10:05 UTC\n",
      "2010-12-23T06:52:05 UTC\n",
      "2010-12-23T10:34:05 UTC\n",
      "2010-12-23T14:16:05 UTC\n",
      "2010-12-23T17:58:05 UTC\n",
      "2010-12-23T21:40:05 UTC\n",
      "2010-12-24T01:22:05 UTC\n",
      "2010-12-24T05:04:05 UTC\n",
      "2010-12-24T08:46:05 UTC\n",
      "2010-12-24T12:28:05 UTC\n",
      "2010-12-24T16:10:05 UTC\n",
      "2010-12-24T19:52:05 UTC\n",
      "2010-12-24T23:34:05 UTC\n"
     ]
    }
   ],
   "source": [
    "from anise.time import TimeSeries\n",
    "\n",
    "start = Epoch(\"2010-12-21T03:04:05\")\n",
    "\n",
    "for epoch_from_series in TimeSeries(start=start, end=start + Unit.Day*4, step=Unit.Hour*3.7, inclusive=True):\n",
    "    print(epoch_from_series)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc3f816b-9393-4075-88b0-c1048e5cbb34",
   "metadata": {},
   "source": [
    "Finally, as a last bit of this introduction on time, let's demonstrate how durations can be rounded in hifitime.\n",
    "Most humans don't usually care about nanosecond precision, so we can round those away."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "bd49354c-f5cc-46b8-bb56-4ff485da6b35",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "13 days 5 h 36 min 23 ns\n",
      "13 days\n",
      "13 days 5 h 36 min\n",
      "9131 days 6 h\n"
     ]
    }
   ],
   "source": [
    "orig = Duration(\"13 d 5.6 h 23 ns\")\n",
    "print(orig)\n",
    "\n",
    "print(orig.round(Unit.Day*1)) # Rounded to the nearest day\n",
    "\n",
    "print(orig.floor(Unit.Second*1)) # Floored to the nearest second\n",
    "\n",
    "print(orig.ceil(Unit.Century*0.25)) # Ceiled to the nearest quarter of a century"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0150cf41-df82-43e9-afdb-7281c39ea1fa",
   "metadata": {},
   "source": [
    "### Learning exercise\n",
    "\n",
    "Your exercise is simple: create two dates in different time systems, run a time series over them at your desired time step, and print the difference between the start date and \"current\" date of the iterator as the printing happens. Don't hesitate to look at the [Python user guide of Hifitime](https://nyxspace.com/hifitime/python/).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68eb5caf-252c-4495-a314-31bd1e3e7795",
   "metadata": {},
   "source": [
    "## Frame management\n",
    "\n",
    "Frames are absolutely crucial in astrodynamics. They defined the point of origin and the orientation of that frame for each state.\n",
    "\n",
    "In terms of orientation, one will typically use the J2000 orientation. This is an inertial frame (well, quasi-inertial) which is tilted by about 23 degrees compared to the Ecliptic frame. The Ecliptic frame can be thought of as the average plane in which all of the planets rotate around the Sun. As such, the J2000 frame can be thought of the frame where the Z axis points _on average_ in the direction of the axis of rotation of the Earth.\n",
    "\n",
    "If a state is communicated in the \"J2000 frame\" without any centeral object provided, it's safe to assume that's the \"EME2000\" frame, which is the Earth Mean Equator J2000 orientation frame.\n",
    "\n",
    "There are _lots_ of frames, since one can define any frame with any center and any orientation. ANISE tries to help out: you can create a frame simply specifying its central object (as it NAIF ID) and its orientation ID (also matches the NAIF IDs). ANISE also provides constants for most common frames."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "286bfaf7-8ffd-4f77-bc61-c41b03f52983",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['EARTH_ECLIPJ2000',\n",
       " 'EARTH_ITRF93',\n",
       " 'EARTH_J2000',\n",
       " 'EARTH_MOON_BARYCENTER_J2000',\n",
       " 'EME2000',\n",
       " 'IAU_EARTH_FRAME',\n",
       " 'IAU_JUPITER_FRAME',\n",
       " 'IAU_MARS_FRAME',\n",
       " 'IAU_MERCURY_FRAME',\n",
       " 'IAU_MOON_FRAME',\n",
       " 'IAU_NEPTUNE_FRAME',\n",
       " 'IAU_SATURN_FRAME',\n",
       " 'IAU_URANUS_FRAME',\n",
       " 'IAU_VENUS_FRAME',\n",
       " 'JUPITER_BARYCENTER_J2000',\n",
       " 'MARS_BARYCENTER_J2000',\n",
       " 'MERCURY_J2000',\n",
       " 'MOON_J2000',\n",
       " 'NEPTUNE_BARYCENTER_J2000',\n",
       " 'PLUTO_BARYCENTER_J2000',\n",
       " 'SATURN_BARYCENTER_J2000',\n",
       " 'SSB_J2000',\n",
       " 'SUN_J2000',\n",
       " 'URANUS_BARYCENTER_J2000',\n",
       " 'VENUS_J2000',\n",
       " '__class__',\n",
       " '__delattr__',\n",
       " '__dir__',\n",
       " '__doc__',\n",
       " '__eq__',\n",
       " '__format__',\n",
       " '__ge__',\n",
       " '__getattribute__',\n",
       " '__getstate__',\n",
       " '__gt__',\n",
       " '__hash__',\n",
       " '__init__',\n",
       " '__init_subclass__',\n",
       " '__le__',\n",
       " '__lt__',\n",
       " '__module__',\n",
       " '__ne__',\n",
       " '__new__',\n",
       " '__reduce__',\n",
       " '__reduce_ex__',\n",
       " '__repr__',\n",
       " '__setattr__',\n",
       " '__sizeof__',\n",
       " '__str__',\n",
       " '__subclasshook__']"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from anise.astro.constants import Frames\n",
    "dir(Frames) # List of built-in frames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "e87a0926-8018-4ef8-a581-2d8274e19f9e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Earth-Moon Barycenter J2000\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "# One can build a frame from its ephemeris and orientation IDs, which are the same as the NAIF codes.\n",
    "# For example, ephemeris center 3 is the Earth Moon Barycenter, and orientation 1 is J2000\n",
    "from anise.astro import Frame\n",
    "custom_frame = Frame(3, 1)\n",
    "print(custom_frame)\n",
    "print(custom_frame == Frames.EARTH_MOON_BARYCENTER_J2000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "f3e821b2-78f4-4acb-8592-10a741d39674",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['__class__',\n",
       " '__delattr__',\n",
       " '__dir__',\n",
       " '__doc__',\n",
       " '__eq__',\n",
       " '__format__',\n",
       " '__ge__',\n",
       " '__getattribute__',\n",
       " '__getnewargs__',\n",
       " '__getstate__',\n",
       " '__gt__',\n",
       " '__hash__',\n",
       " '__init__',\n",
       " '__init_subclass__',\n",
       " '__le__',\n",
       " '__lt__',\n",
       " '__module__',\n",
       " '__ne__',\n",
       " '__new__',\n",
       " '__reduce__',\n",
       " '__reduce_ex__',\n",
       " '__repr__',\n",
       " '__setattr__',\n",
       " '__sizeof__',\n",
       " '__str__',\n",
       " '__subclasshook__',\n",
       " 'ephem_origin_id_match',\n",
       " 'ephem_origin_match',\n",
       " 'ephemeris_id',\n",
       " 'flattening',\n",
       " 'is_celestial',\n",
       " 'is_geodetic',\n",
       " 'mean_equatorial_radius_km',\n",
       " 'mu_km3_s2',\n",
       " 'orient_origin_id_match',\n",
       " 'orient_origin_match',\n",
       " 'orientation_id',\n",
       " 'polar_radius_km',\n",
       " 'semi_major_radius_km',\n",
       " 'shape',\n",
       " 'with_ephem',\n",
       " 'with_orient']"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eme2k = Frames.EARTH_J2000\n",
    "dir(eme2k)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f460a3de-30ba-4935-9560-98fa1097de93",
   "metadata": {},
   "source": [
    "Frames include _potentially_ a lot of information. This is because the frame data is used to compute orbital elements, which we'll cover later on.\n",
    "\n",
    "However, ANISE does not provide any defaults for this data. You must load a Planetary Constants ANISE (\"PCA\") file for ANISE to know that these exist. Otherwise, you'll get an exception:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "c9a2eb6c-0783-49d8-ab5b-8cf914f5803e",
   "metadata": {},
   "outputs": [
    {
     "ename": "Exception",
     "evalue": "retrieving gravitational parameter requires the frame Earth J2000 to have mu_km3_s2 defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mException\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[14], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[43meme2k\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mmu_km3_s2\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n",
      "\u001b[0;31mException\u001b[0m: retrieving gravitational parameter requires the frame Earth J2000 to have mu_km3_s2 defined"
     ]
    }
   ],
   "source": [
    "eme2k.mu_km3_s2()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "ad5af23d-3a04-4122-ac51-a818bf4451e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's load the PCA\n",
    "almanac = almanac.load(\"../../data/pck08.pca\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca9cf9b1-bb8b-40b8-8364-de9fbf4f0aae",
   "metadata": {},
   "source": [
    "**VERY IMPORTANT:** Note how we've _reassigned_ the `almanac` variable in the `load` call. This is _also_ part of the thread safety: the almanac structure is _read only_ and the reassignment is what allows us to grab a reference to the almanac that has the newly loaded file.\n",
    "\n",
    "(Note that this also means you can have _multiple_ Almanac objects and compare the data loaded into them, similar to a \"BOA\" for those familiar with JPL MONTE.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "f4374ba5-2c0c-4768-869c-cdb6edb7651b",
   "metadata": {},
   "outputs": [
    {
     "ename": "Exception",
     "evalue": "retrieving gravitational parameter requires the frame Earth J2000 to have mu_km3_s2 defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mException\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[16], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[43meme2k\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mmu_km3_s2\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m)\n",
      "\u001b[0;31mException\u001b[0m: retrieving gravitational parameter requires the frame Earth J2000 to have mu_km3_s2 defined"
     ]
    }
   ],
   "source": [
    "print(eme2k.mu_km3_s2())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e550acba-dbe5-45cb-bb41-cd627b1d07b3",
   "metadata": {},
   "source": [
    "Uh! One way for ANISE to be thread safe without mutexes is for the data to be stored solely in the almanac structure. So loading the data in the almanac is not sufficient: we need to retrieve it.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "43ac7a62-2c1e-4c34-9b62-4646079c11d0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Earth J2000 (μ = 398600.435436096 km^3/s^2, eq. radius = 6378.14 km, polar radius = 6356.75 km, f = 0.0033536422844278)\n"
     ]
    }
   ],
   "source": [
    "loaded_eme2k = almanac.frame_info(eme2k)\n",
    "print(loaded_eme2k)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12a7a0d9-f4c4-460e-b059-999fe8dd9ab0",
   "metadata": {},
   "source": [
    "All of this data was created from the KPL/TPC files that NAIF provides. Another tutorial will guide you through modifying these."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc8ce818-393f-4781-8343-4f04818747e9",
   "metadata": {},
   "source": [
    "### Exercise #1\n",
    "\n",
    "1. Investigate what are some of the available fields in this loaded EME2000 frame.\n",
    "2. Load the Venus J2000 frame and repeat point 1."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d09ea2e7-b41c-4c5e-88fa-8a64c098fa26",
   "metadata": {},
   "source": [
    "## Querying the almanac\n",
    "\n",
    "So we've already seen how to build an epoch, and how to find a frame. Let's now quickly look into aberration corrections. The best explanation you'll find is that from NAIF itself: <https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/abcorr.html>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "6ecc0641-ec60-4c62-98b0-efecfbe6ee2c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[0;31mInit signature:\u001b[0m \u001b[0mAberration\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mname\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",
       "Represents the aberration correction options in ANISE.\n",
       "\n",
       "In space science and engineering, accurately pointing instruments (like optical cameras or radio antennas) at a target is crucial. This task is complicated by the finite speed of light, necessitating corrections for the apparent position of the target.\n",
       "\n",
       "This structure holds parameters for aberration corrections applied to a target's position or state vector. These corrections account for the difference between the target's geometric (true) position and its apparent position as observed.\n",
       "\n",
       "# Rule of tumb\n",
       "In most Earth orbits, one does _not_ need to provide any aberration corrections. Light time to the target is less than one second (the Moon is about one second away).\n",
       "In near Earth orbits, e.g. inner solar system, preliminary analysis can benefit from enabling unconverged light time correction. Stellar aberration is probably not required.\n",
       "For deep space missions, preliminary analysis would likely require both light time correction and stellar aberration. Mission planning and operations will definitely need converged light-time calculations.\n",
       "\n",
       "For more details, <https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/abcorr.html>.\n",
       "\n",
       "# SPICE Validation\n",
       "\n",
       "The validation test `validate_jplde_de440s_aberration_lt` checks 101,000 pairs of ephemeris computations and shows that the unconverged Light Time computation matches the SPICE computations almost all the time.\n",
       "More specifically, the 99th percentile of error is less than 5 meters, the 75th percentile is less than one meter, and the median error is less than 2 millimeters.\n",
       "\u001b[0;31mFile:\u001b[0m           ~/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages/anise/__init__.py\n",
       "\u001b[0;31mType:\u001b[0m           type\n",
       "\u001b[0;31mSubclasses:\u001b[0m     "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from anise import Aberration\n",
    "Aberration?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "4410e354-d9dc-42a8-94b9-09cffd57dc96",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LT (@0x7fd4c0478ee0)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lt_aber = Aberration(\"LT\") # This is the same as the LT enum in SPICE.\n",
    "lt_aber"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "935bf989-bf59-49ba-b31f-8b64267f26d1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "unconverged light-time aberration\n"
     ]
    }
   ],
   "source": [
    "# If you want the details of what will be computed, using the string representation instead of the `__repr__`:\n",
    "print(lt_aber)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "fba7ff18-4669-42d1-92c9-00e3568c736c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[0;31mSignature:\u001b[0m \u001b[0malmanac\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtranslate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtarget_frame\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mobserver_frame\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mepoch\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mab_corr\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\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",
       "Returns the Cartesian state of the target frame as seen from the observer frame at the provided epoch, and optionally given the aberration correction.\n",
       "\n",
       "# SPICE Compatibility\n",
       "This function is the SPICE equivalent of spkezr: `spkezr(TARGET_ID, EPOCH_TDB_S, ORIENTATION_ID, ABERRATION, OBSERVER_ID)`\n",
       "In ANISE, the TARGET_ID and ORIENTATION are provided in the first argument (TARGET_FRAME), as that frame includes BOTH\n",
       "the target ID and the orientation of that target. The EPOCH_TDB_S is the epoch in the TDB time system, which is computed\n",
       "in ANISE using Hifitime. THe ABERRATION is computed by providing the optional Aberration flag. Finally, the OBSERVER\n",
       "argument is replaced by OBSERVER_FRAME: if the OBSERVER_FRAME argument has the same orientation as the TARGET_FRAME, then this call\n",
       "will return exactly the same data as the spkerz SPICE call.\n",
       "\n",
       "# Warning\n",
       "This function only performs the translation and no rotation whatsoever. Use the `transform` function instead to include rotations.\n",
       "\n",
       "# Note\n",
       "This function performs a recursion of no more than twice the [MAX_TREE_DEPTH].\n",
       "\u001b[0;31mType:\u001b[0m      builtin_function_or_method"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Looking at the documentation here can be quite useful\n",
    "almanac.translate?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "a2c9b478-e8b0-411b-828a-0bf50cc82a6b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[0;31mSignature:\u001b[0m \u001b[0malmanac\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtransform\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtarget_frame\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mobserver_frame\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mepoch\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mab_corr\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\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",
       "Returns the Cartesian state needed to transform the `from_frame` to the `to_frame`.\n",
       "\n",
       "# SPICE Compatibility\n",
       "This function is the SPICE equivalent of spkezr: `spkezr(TARGET_ID, EPOCH_TDB_S, ORIENTATION_ID, ABERRATION, OBSERVER_ID)`\n",
       "In ANISE, the TARGET_ID and ORIENTATION are provided in the first argument (TARGET_FRAME), as that frame includes BOTH\n",
       "the target ID and the orientation of that target. The EPOCH_TDB_S is the epoch in the TDB time system, which is computed\n",
       "in ANISE using Hifitime. THe ABERRATION is computed by providing the optional Aberration flag. Finally, the OBSERVER\n",
       "argument is replaced by OBSERVER_FRAME: if the OBSERVER_FRAME argument has the same orientation as the TARGET_FRAME, then this call\n",
       "will return exactly the same data as the spkerz SPICE call.\n",
       "\n",
       "# Note\n",
       "The units will be those of the underlying ephemeris data (typically km and km/s)\n",
       "\u001b[0;31mType:\u001b[0m      builtin_function_or_method"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# The transform function will also compute the rotations where needed.\n",
    "almanac.transform?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "aa3e6682-c0d0-4335-a488-69e9eef0dccf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Solar System Barycenter J2000] 2010-12-21T03:04:05 UTC\tposition = [2006351.284308, 135110969.586735, 58576035.963861] km\tvelocity = [-30.254436, 0.375013, 0.163924] km/s\n"
     ]
    }
   ],
   "source": [
    "state = almanac.translate(Frames.EARTH_J2000, Frames.SSB_J2000, epoch, None)\n",
    "print(state)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8463e96-1b38-49f9-bdf2-36137e3b5212",
   "metadata": {},
   "source": [
    "**TA DA!** We've just queried the almanac (i.e. the DE440s.bsp file) to get the position of the Earth in the J2000 frame as seen from the solar system barycenter (\"SSB\") in the J2000 frame, at the epoch we defined earlier, and when no aberration correction is performed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "0ab91d31-1792-46ca-921f-e3320c226835",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8 min 11 s 259 ms 84 μs 389 ns\n"
     ]
    }
   ],
   "source": [
    "state = almanac.translate(Frames.EARTH_J2000, Frames.SSB_J2000, epoch, Aberration(\"LT\"))\n",
    "print(state.light_time())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9644e04-1f62-4694-9d27-a70e41ad9f43",
   "metadata": {},
   "source": [
    "Calling \"light_time\" will return a Hifitime duration instance. You can grab the exact seconds by converting that high precision duration into the unit of your choice:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "7e50475e-d2a1-44f5-82bb-b721165c2625",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.005685869032280092\n",
      "491.259084389\n",
      "491.259084389\n"
     ]
    }
   ],
   "source": [
    "print(state.light_time().to_unit(Unit.Day))\n",
    "print(state.light_time().to_unit(Unit.Second))\n",
    "print(state.light_time().to_seconds())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84d60190-cc88-4ab8-86b5-cad83297f787",
   "metadata": {},
   "source": [
    "## Exerices\n",
    "\n",
    "1. Query the position of the Earth at different times using a time series. Note that the DE440s file is only valid for 200 years centered on 1950.\n",
    "2. Query the position of Venus as seen from the Earth, both in the J2000 frame, and make a plot of the light time between these two planets throughout a given year."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6f54988-411b-456d-b3b6-0e1cb5902a95",
   "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
}