{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Improved state estimation with MPC\n", "Copyright (c) 2010-2024, Delft University of Technology. All rights reserved. This file is part of the Tudat. Redistribution and use in source and binary forms, with or without modification, are permitted exclusively under the terms of the Modified BSD license. You should have received a copy of the license with this file. If not, please visit: http://tudat.tudelft.nl/LICENSE.\n", "\n", "## Objectives\n", "This example extends the previous [Initial state estimation with Minor Planet Center Observations](estimation_with_mpc.ipynb). In an attempt to improve the results from the previous example, we introduce and compare the effects of including satellite data, star catalog corrections, observation weighting and more expansive acceleration models. It essential to be familiar with the previous example as many concepts will be reused here without explanation. \n", "\n", "As in the previous example we will estimate the initial state of [433 Eros](https://en.wikipedia.org/wiki/433_Eros). In addition to observation data from MPC and metadata from SBDB, we now also use ephemeris data from JPL Horizons to retrieve position data for observing space telescopes, additional perturbing bodies and as a method of comparison. This is accomplished using Tudat's HorizonsQuery Interface." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Import statements" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "# Tudat imports for propagation and estimation\n", "from tudatpy.interface import spice\n", "from tudatpy.dynamics import environment_setup, parameters_setup, parameters, propagation, propagation_setup\n", "from tudatpy.estimation import observable_models_setup,observable_models, observations_setup, observations, estimation_analysis\n", "from tudatpy.constants import GRAVITATIONAL_CONSTANT\n", "from tudatpy.astro.frame_conversion import inertial_to_rsw_rotation_matrix\n", "from tudatpy.astro.time_representation import DateTime\n", "from tudatpy.astro import element_conversion\n", "\n", "# import MPC, SBDB and Horizons interface\n", "from tudatpy.data.mpc import BatchMPC\n", "from tudatpy.data.horizons import HorizonsQuery\n", "from tudatpy.data.sbdb import SBDBquery\n", "\n", "\n", "# other useful modules\n", "import numpy as np\n", "import datetime\n", "import pandas as pd\n", "\n", "import matplotlib.pyplot as plt\n", "import matplotlib.cm as cm\n", "from tudatpy.astro import time_representation\n", "from tudatpy.astro.time_representation import DateTime\n", "from astropy.table import Table\n", "\n", "# SPICE KERNELS\n", "spice.load_standard_kernels()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Preparing the environment and observations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Setting the input constants\n", "Let's setup some constants that are used throughout the tutorial. The MPC code for Eros is 433. We also set a start and end date for our observations, the number of iterations for our estimation, a timestep for our integrator and a 1 month buffer to avoid interpolation errors in our analysis.\n", "\n", "We use a spice kernel to get a guess for our initial state and to check our estimation afterwards. The default spice kernel `codes_300ast_20100725.bsp` contains many popular asteroids, however they are not all identified by name (433 Eros is `\"Eros\"` but 16 Psyche is `\"2000016\"` etc.). To ensure this example works dynamically, for any single MPC code as input we use the SDBD to retrieve the name and SPK-ID used for the spice kernel.\n", "\n", "For our frame origin we use the Solar System Barycentre. The data from MPC is presented in the J2000 reference frame, currently BatchMPC does not support conversion to other reference frames and as such we match it in our environment. \n", "\n", "For this extended example, a longer observation period of 9 years is used." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Direct inputs:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "target_mpc_code = \"433\"\n", "\n", "observations_start = datetime.datetime(2015, 1, 1)\n", "observations_end = datetime.datetime(2024, 1, 1)\n", "\n", "# number of iterations for our estimation\n", "number_of_pod_iterations = 6\n", "\n", "# timestep of 24 hours for our estimation\n", "timestep_global = 24 * 3600.0\n", "\n", "# 2 month time buffer used to avoid interpolation errors:\n", "time_buffer = 2 * 31 * 86400.0\n", "\n", "# define the frame origin and orientation.\n", "global_frame_origin = \"SSB\"\n", "global_frame_orientation = \"J2000\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Derived inputs:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "SPK ID for 433 Eros is: Eros\n" ] } ], "source": [ "target_sbdb = SBDBquery(target_mpc_code)\n", "\n", "mpc_codes = [target_mpc_code] # the BatchMPC interface requires a list.\n", "target_spkid = target_sbdb.codes_300_spkid # the ID used by the\n", "target_name = target_sbdb.shortname # the ID used by the\n", "\n", "print(f\"SPK ID for {target_name} is: {target_spkid}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Eros Ephemeris Uncertainty\n", "\n", "Additionally, we will retrieve the published ephemeris uncertainty from JPL Horizons.\n", "At the moment, this is not directly supported through Tudat interfaces, but will be added in a future release.\n", "In this case, the ephemeris uncertainty of Eros has been downloaded manually and is provided in the [data/Eros-Ephemeris-Uncertainty.ecsv](data/Eros-Ephemeris-Uncertainty.ecsv) file, using the code in the following `astroquery` PR: https://github.com/astropy/astroquery/pull/3273" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "targetname, datetime_jd, datetime_str, H, G, x, y, z, vx, vy, vz, x_s, y_s, z_s, vx_s, vy_s, vz_s, r_s, t_s, n_s, vr_s, vt_s, vn_s\n" ] } ], "source": [ "uncertainty_file_path = \"data/Eros-Ephemeris-Uncertainty.ecsv\"\n", "ephemeris_uncertainty_table = Table.read(uncertainty_file_path)\n", "\n", "print(\", \".join(ephemeris_uncertainty_table.colnames))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As the data provided by JPL is not in SI units, we will first convert it to units compatible with Tudat." ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "ephemeris_uncertainty_table[\"ephemeris_time\"] = [\n", " julian_day_to_seconds_since_epoch(\n", " jd,\n", " )\n", " for jd in ephemeris_uncertainty_table[\"datetime_jd\"]\n", "]\n", "\n", "for col in ephemeris_uncertainty_table.colnames:\n", " if col in [\"ephemeris_time\", \"datetime_jd\", \"datetime_str\", \"targetname\", \"H\", \"G\"]:\n", " continue\n", " ephemeris_uncertainty_table[col] = ephemeris_uncertainty_table[col].quantity.si" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Combinations and additional body setup\n", "There are various ways to change our estimation. We can create a system of setups to compare those various options and to facilitate comparison. Throughout the example, the following options are considered:\n", "\n", "- [`accel_levels`] Different acceleration settings, for this example, 3 options are created in increasing order of realism\n", "\n", " - LVL 1 - Only point-mass gravity for the sun and the 8 mayor planets as well as Schwarzschild relativistic correction for the sun.\n", " - LVL 2 - LVL 1 + point-mass gravity for the mayor moons of Jupiter, Saturn, Earth and Mars.\n", " - LVL 3 - LVL 2 + SHG for the Earth and point-mass gravity for Triton, Titania, Pluto, the mayor Near Earth Asteroids (NEA) and largest Main Body Asteroids (MBA). These additional bodies are retrieved through the JPL Horizons interface.\n", "\n", "- [`use_sat_data`] Observations by space telescope WISE\n", "- [`use_catalog_cor`] Star catalog corrections as described in \"Star catalog position and proper motion corrections in asteroid astrometry II: The Gaia era\" by Eggl et al.\n", "- [`use_weighting`] Estimation weights as described in \"Statistical analysis of astrometric errors for the most productive asteroid surveys\" by Veres et al.\n", "\n", "A function, `perform_estimation`, is be created below which will perform the estimation based on these settings. The settings are described by a series of lists below, with a list of setup names to describe them.\n", "\n", "The acceleration model is expected to have the most effect on the simulation. For the first round of comparison, only the acceleration models will be changed with the remainder all set to False. Three setups are constructed below. We also define constants to later set up satellite data." ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "setup_names = [\"LVL1 Accelerations\", \"LVL2 Accelerations\", \"LVL3 Accelerations\"]\n", "\n", "accel_levels = [1, 2, 3]\n", "use_sat_data = [False, False, False]\n", "use_catalog_cor = [False, False, False]\n", "use_weighting = [False, False, False]\n", "\n", "satellites_names = [\"WISE\"]\n", "satellites_MPC_codes = [\"C51\"] # C51 is the observatory code MPC uses for WISE\n", "satellites_Horizons_codes = [\n", " \"-163\"\n", "] # -163 is the query ID for WISE in Horizons see explanation below.\n", "\n", "\n", "# Consider trying out different combinations of satellites.\n", "# Note that you must change the dates to use TESS as it launched in April 2018\n", "# satellites_names = [\"WISE\", \"TESS\"]\n", "# satellites_MPC_codes = [\"C51\", \"C57\"]\n", "# satellites_Horizons_codes = [\"-163\", \"-95\"]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For LVL3 accelerations, the point-mass gravitational acceleration of Pluto, Triton and Titania are added using JPL Horizons. Horizons only provides an ephemeris, the masses are retrieved and added manually. Note that JPL Horizons has a unique querying scheme in which Pluto is best accessed using the ID 999. The API documentation for the `HorizonsQuery()` class provides an extensive but not exhaustive explanation of these IDs. For now it is sufficient to understand that mayor bodies such as Earth are denoted `399` (3rd mayor body), Asteroids/Minor bodies are denoted with a semicolon like `433;` for Eros (MPC code + ;), and satellites are denote with a minus sign like `-163` for WISE.\n", "\n", "JPL Horizons will also be used to retrieve the ephemeris for mayor NEA and MBA. Again their masses will be added through other means, in this case we use [SiMDA](https://astro.kretlow.de/simda/), which is an archive of published mass and diameter estimates for minor bodies. \n", "\n", "All NEAs from the archive are retrieved, as well as all MBA with a mass greater than 1e20 kg. Consider altering this filter to see the effects." ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "lvl3_extra_bodies = [\"999\", \"Triton\", \"Titania\"] # here 999 is Pluto in JPL Horizons\n", "lvl3_extra_bodies_masses = [1.3025e22, 2.1389e22, 3.4550e21]" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of additional bodies from SiMDA: 16\n" ] }, { "data": { "text/html": [ "
| \n", " | NUM | \n", "DESIGNATION | \n", "DIAM | \n", "DYN | \n", "MASS | \n", "
|---|---|---|---|---|---|
| 19 | \n", "2 | \n", "Pallas | \n", "520.8 | \n", "MBA | \n", "2.050000e+20 | \n", "
| 384 | \n", "1036 | \n", "Ganymed | \n", "35.7 | \n", "NEA | \n", "6.580000e+16 | \n", "
| 395 | \n", "3671 | \n", "Dionysus | \n", "0.9 | \n", "NEA | \n", "8.380000e+11 | \n", "
| 398 | \n", "5381 | \n", "Sekhmet | \n", "1.0 | \n", "NEA | \n", "1.040000e+12 | \n", "
| 399 | \n", "25143 | \n", "Itokawa | \n", "0.3 | \n", "NEA | \n", "3.500000e+10 | \n", "
| 401 | \n", "35107 | \n", "1991 VH | \n", "1.1 | \n", "NEA | \n", "1.400000e+12 | \n", "
| 407 | \n", "65803 | \n", "Didymos | \n", "0.8 | \n", "NEA | \n", "5.240000e+11 | \n", "
| 408 | \n", "66063 | \n", "1998 RO1 | \n", "0.7 | \n", "NEA | \n", "3.600000e+11 | \n", "
| 409 | \n", "66391 | \n", "1999 KW4 | \n", "1.3 | \n", "NEA | \n", "2.350000e+12 | \n", "
| 418 | \n", "136617 | \n", "1994 CC | \n", "0.6 | \n", "NEA | \n", "2.590000e+11 | \n", "
| 421 | \n", "164121 | \n", "2003 YT1 | \n", "1.1 | \n", "NEA | \n", "1.270000e+12 | \n", "
| 422 | \n", "175706 | \n", "1996 FG3 | \n", "1.8 | \n", "NEA | \n", "4.270000e+12 | \n", "
| 423 | \n", "185851 | \n", "2000 DP107 | \n", "1.6 | \n", "NEA | \n", "4.600000e+11 | \n", "
| 424 | \n", "276049 | \n", "2002 CE26 | \n", "3.5 | \n", "NEA | \n", "1.950000e+13 | \n", "
| 425 | \n", "311066 | \n", "2004 DC | \n", "0.3 | \n", "NEA | \n", "3.570000e+10 | \n", "
| 426 | \n", "494658 | \n", "2000 UG11 | \n", "0.3 | \n", "NEA | \n", "9.350000e+09 | \n", "