{ "cells": [ { "cell_type": "markdown", "id": "0", "metadata": {}, "source": [ "# Propagator Comparisons (IAS15 & DOP853)\n", "---\n", "Last revised by Z. Ellis on 2026 MAR 30\n", "\n", "## Objectives\n", "This tutorial will demonstrate the usage of the two different integrators supported within Scarabaeus' Propagator class, DOP853 and IAS15, across increasingly complex dynamics configurations. The configurations are:\n", "\n", "1) Keplerian orbit\n", "2) Keplerian orbit + third body perturbations\n", "3) Keplerian orbit + third body perturbations + cannonball solar radiation pressure model\n", "4) Keplerian orbit + third body perturbations + N-plate solar radiation pressure model\n", "5) Keplerian orbit + third body perturbations + N-plate solar radiation pressure model + spherical harmonics gravity\n", "\n", "The differences are highlighted by examining an inclined, elliptical geocentric orbit and noting where and why one integrator calculates different results than the other.\n", "Performance metrics like computational efficiency are also examined." ] }, { "cell_type": "markdown", "id": "1", "metadata": {}, "source": [ "## Imports and Set Up\n", "\n", "Here we'll import the necessary libraries, define units and frames, and load a metakernel for SPICE time conversions." ] }, { "cell_type": "code", "execution_count": null, "id": "2", "metadata": {}, "outputs": [], "source": [ "import scarabaeus as scb\n", "\n", "import os, time\n", "import pyrootutils\n", "\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "os.chdir(pyrootutils.setup_root(os.getcwd(), indicator = 'pyproject.toml', pythonpath = True))\n", "\n", "## units, frames, kernels\n", "kg, km, sec = scb.Units.get_units(['kg', 'km', 'sec'])\n", "J2000 = scb.Frame('J2000')\n", "\n", "scb.SpiceManager.load_kernel_from_mkfile(os.getcwd() + '/data/kernels/locked/locked_generic.tm')" ] }, { "cell_type": "markdown", "id": "3", "metadata": {}, "source": [ "## Define Propagation Parameters\n", "\n", "Next, we'll create an object to represent Earth, the body we'll be orbiting in this tutorial, and its spherical harmonics information `sh_info`. We also need to create an orbiting spacecraft `sc` and an N-plate model for it.\n", "\n", "We'll also define the time interval we want to propagate over and a few orbital elements." ] }, { "cell_type": "code", "execution_count": null, "id": "4", "metadata": {}, "outputs": [], "source": [ "## earth and spacecraft\n", "earth = scb.CelestialBody.from_constants(\"EARTH\")\n", "sh_info = {'order' : 3,\n", " 'cs_file' : 'data/dynamic_setup/sph_coefficients/Earth_100.json',\n", " 'body' : scb.EARTH,\n", " 'norm_flag': True}\n", "\n", "n_plate = scb.nPlateModel('data/dynamic_setup/nplate_coefficients/VandV_case5a_nPlate_config.json')\n", "sc = scb.Spacecraft(name = 'ORCCA_SC',\n", " spice_id = -1000,\n", " tot_mass = scb.ArrayWUnits(2000, kg),\n", " area = scb.ArrayWUnits(1e-06, km**2),\n", " ref_coeff = scb.ArrayWUnits(1.5, None),\n", " n_plate_model = n_plate,\n", " attitudeMode = 'nadir_pointing_to_sun',\n", " nadirPointingAxis = np.array([-1, 0, 0]))\n", "\n", "## orbital parameters\n", "t0 = scb.SpiceManager.cal2et('2026 JAN 01 00:00:00.000')\n", "tf = scb.SpiceManager.cal2et('2026 JAN 01 10:00:00.000') # 10 hours\n", "epochs = scb.EpochArray(np.arange(t0, tf, 500), 'TDB')\n", "\n", "# semi major axis, ecccentricity, inclination, gravitational parameter\n", "a, e, i, mu = 25000, 0.5, np.deg2rad(35), earth.grav_param.values" ] }, { "cell_type": "markdown", "id": "5", "metadata": {}, "source": [ "We want to examine how our integrators differ from one another, and the point that we'll likely see the largest difference will be at periapsis where the spacecraft is moving fastest. In order to see this better, we want to examine a few different trajectories with initial conditions at different points in the orbit. We'll create a function `state_at_nu` which will find the state corresponding to a given true anomaly nu so we can test these different initial true anomalies when we start propagating." ] }, { "cell_type": "code", "execution_count": null, "id": "6", "metadata": {}, "outputs": [], "source": [ "def state_at_nu(nu):\n", " \"\"\"\n", " Function to calculate state at given true anomaly nu.\n", " \"\"\"\n", " # perifocal position and velocity\n", " h = np.sqrt(mu*a*(1 - e**2)) # angular momentum \n", " r = h**2 / (mu*(1 + e*np.cos(nu)))\n", " r0 = np.array([r*np.cos(nu), r*np.sin(nu), 0])\n", " v0 = np.array([-(mu/h)*np.sin(nu), (mu/h) * (e + np.cos(nu)), 0])\n", "\n", " # rotate to intertial intertial\n", " peri_to_n = np.array([[1, 0, 0],\n", " [0, np.cos(i), -np.sin(i)],\n", " [0, np.sin(i), np.cos(i)]])\n", " r0 = (peri_to_n @ r0.reshape((3, 1))).flatten()\n", " v0 = (peri_to_n @ v0.reshape((3, 1))).flatten()\n", "\n", " # convert to AWF and return\n", " pos0 = scb.ArrayWFrame(r0, km , J2000)\n", " vel0 = scb.ArrayWFrame(v0, km/sec, J2000)\n", "\n", " return (pos0, vel0)" ] }, { "cell_type": "markdown", "id": "7", "metadata": {}, "source": [ "## Set up Propagation (and Plotting)\n", "Since we want to look at a few different propagation configurations, we'll define another function that can:\n", "1) propagate any arbitrary dynamics configuration\n", "2) plot the resulting propagation\n", "\n", "We'll do this under `prop_and_plot_fm`, which will first initialize a figure, then initialize four different initial conditions at different true anomalies using `state_at_nu`. For each true anomaly, it will then propagate the initial conditions using both the DOP853 and the IAS15 integrator and compare their differences as percents for each state component. We'll also return the average percent differences so that we can look at how each set of dynamics changes the results." ] }, { "cell_type": "code", "execution_count": null, "id": "8", "metadata": {}, "outputs": [], "source": [ "## plotting\n", "def prop_and_plot_fm(model, title):\n", " \"\"\"\n", " Function to propagate the given dynamics using both DOP853 and IAS15 \n", " integrators. Plots the results and returns the average percent errors \n", " for each state component.\n", "\n", " Also returns average computation time for both integrators.\n", " \"\"\"\n", " # initialize plot\n", " fig, axes = plt.subplot_mosaic([['o', 'o', '$X$', '$Y$', '$Z$'],\n", " ['o', 'o', '$\\dot{X}$', '$\\dot{Y}$', '$\\dot{Z}$']],\n", " per_subplot_kw = {'o' : {'projection' : '3d'}},\n", " constrained_layout = True, figsize = (17, 9))\n", " \n", " # examine four different initial true anomaly cases\n", " nus = [-90, 0, 90, 180] # approaching peri, at peri, approaching apo, at apo\n", "\n", " # save computation time for each propagator and differences between them\n", " tot_avg_pos_err, tot_avg_vel_err, times, traj_errs = [], [], [], []\n", " dop_time, ias_time = [], []\n", " for j, nu in enumerate(nus):\n", " # get state that corresponds to this case's true anomaly\n", " pos0, vel0 = state_at_nu(np.rad2deg(nu))\n", " state = [(\"position\", 3, \"estimated\", \"dynamic\", sc, pos0),\n", " (\"velocity\", 3, \"estimated\", \"dynamic\", sc, vel0)]\n", " state_vector = scb.StateArray(epoch = epochs[0],\n", " origin = earth,\n", " state = scb.StateDefinition.from_components(state))\n", " \n", " # propagate with DOP853\n", " dop_start = time.time()\n", " prop_dop = scb.Propagator(primary_body = sc,\n", " state_vector = state_vector,\n", " tspan = epochs,\n", " integrator = 'DOP853',\n", " force_models = model,\n", " propagate_STM = False)\n", " prop_dop.propagate()\n", " dop_time.append(time.time() - dop_start)\n", " ys_dop = prop_dop.ys\n", "\n", " # and IAS15\n", " ias_start = time.time()\n", " prop_ias = scb.Propagator(primary_body = sc,\n", " state_vector = state_vector,\n", " tspan = epochs,\n", " integrator = 'IAS15',\n", " force_models = model,\n", " propagate_STM = False)\n", " prop_ias.propagate()\n", " ias_time.append(time.time() - ias_start)\n", " ts_ias = prop_ias.times.times.values\n", " ys_ias = prop_ias.ys\n", " times.append(ts_ias)\n", "\n", " ## plot 3d orbit\n", " # IAS15\n", " if j == 0: # only need legend labels for one\n", " axes['o'].plot(ys_ias[:, 0], ys_ias[:, 1], ys_ias[:, 2], \n", " 'b', label = 'IAS15')\n", " axes['o'].plot(ys_ias[0, 0], ys_ias[0, 1], ys_ias[0, 2],\n", " 'go', label = 'IAS15 $t_0$')\n", " axes['o'].plot(ys_ias[-1, 0], ys_ias[-1, 1], ys_ias[-1, 2],\n", " 'rs', label = 'IAS15 $t_f$')\n", " else:\n", " axes['o'].plot(ys_ias[:, 0], ys_ias[:, 1], ys_ias[:, 2], \n", " linewidth = 1.4)\n", " axes['o'].plot(ys_ias[0, 0], ys_ias[0, 1], ys_ias[0, 2], 'go')\n", " axes['o'].plot(ys_ias[-1, 0], ys_ias[-1, 1], ys_ias[-1, 2], 'rs')\n", " \n", " # note which true anomaly only for IAS15\n", " match nu:\n", " case -90: offset = (0, 1000, 0, 1000)\n", " case 0 : offset = (1000, -1000, 1000, 1000)\n", " case 90 : offset = (0, 1000, 0, 2000)\n", " case 180: offset = (0, -2000, 0, 1000)\n", " axes['o'].text(ys_ias[0, 0]+offset[0], ys_ias[0, 1], ys_ias[0, 2]+offset[1],\n", " f'\\u03bd$_0 =$ {nu}$^\\circ$')\n", " axes['o'].text(ys_ias[-1, 0]+offset[2], ys_ias[-1, 1], ys_ias[-1, 2]+offset[3],\n", " f'\\u03bd$_0 =$ {nu}$^\\circ$')\n", "\n", " # DOP853\n", " if j == 0: # only need legend labels for one\n", " axes['o'].plot(ys_dop[:, 0], ys_dop[:, 1], ys_dop[:, 2],\n", " 'k--', alpha = 0.6, label = 'DOP853')\n", " axes['o'].plot(ys_dop[0, 0], ys_dop[0, 1], ys_dop[0, 2],\n", " 'm.', label = 'DOP853 $t_0$')\n", " axes['o'].plot(ys_dop[-1, 0], ys_dop[-1, 1], ys_dop[-1, 2],\n", " 'b+', label = 'DOP853 $t_f$')\n", " else:\n", " axes['o'].plot(ys_dop[:, 0], ys_dop[:, 1], ys_dop[:, 2],\n", " '--', linewidth = 0.9)\n", " axes['o'].plot(ys_dop[0, 0], ys_dop[0, 1], ys_dop[0, 2], 'm.')\n", " axes['o'].plot(ys_dop[-1, 0], ys_dop[-1, 1], ys_dop[-1, 2], 'b+')\n", "\n", " # percent errors\n", " err = np.nan_to_num(abs((ys_ias[:] - ys_dop[:, 0:6]) / ys_ias)*100, nan = 0)\n", " traj_errs.append(err)\n", "\n", " # calculate average percent errors and save for this initial true anomaly\n", " avg_pos_err = np.mean(err[:, 0:3])\n", " avg_vel_err = np.mean(err[:, 3:6])\n", " tot_avg_pos_err.append(avg_pos_err), tot_avg_vel_err.append(avg_vel_err)\n", " \n", " # primary body\n", " axes['o'].plot(0, 0, 0, 'ko')\n", "\n", " # formatting\n", " axes['o'].set_title(f'IAS15 Avg Computation Time = {np.mean(ias_time):.3f} sec\\n'\n", " f'DOP853 Avg Computation Time = {np.mean(dop_time):.3f} sec')\n", " axes['o'].legend(bbox_to_anchor = (0.5, -0.15), loc = 'lower center', ncols = 2)\n", " axes['o'].set_aspect('equal'), axes['o'].view_init(elev = 33, azim = -60)\n", "\n", " ## plot percent errors for each initial true case\n", " # plot values with smaller average error in front of ones with larger\n", " zorder_inds = sorted(range(len(traj_errs)), key = lambda err: np.mean(err), reverse = True)\n", " traj_errs = [traj_errs[j] for j in zorder_inds]\n", " nus = [nus[j] for j in zorder_inds]\n", " for t, traj_err, nu in zip(times, traj_errs, nus):\n", " for j, ax in enumerate(axes):\n", " if ax != 'o':\n", " axes[ax].plot((t - t0)/3600, traj_err[:, j - 1], '-o', markersize = 3,\n", " label = f'\\u03bd$_0 =$ {nu}$^\\circ$')\n", " axes[ax].set_xlabel('Time [hrs]')\n", " axes[ax].set_ylabel('IAS15 - DOP853 % Difference')\n", " axes[ax].set_title(ax + f'\\n Mean Difference: {np.mean(err[:, j - 1]):.4e} %')\n", " axes[ax].grid(), axes[ax].legend()\n", "\n", " # finalize plot and return average valeus\n", " fig.suptitle(title)\n", " plt.show()\n", " return (np.mean(tot_avg_pos_err), np.mean(tot_avg_vel_err))" ] }, { "cell_type": "markdown", "id": "9", "metadata": {}, "source": [ "## Dynamics Models and Analysis\n", "Now that we have a function that can propagate and plot any dynamics we give it, we need to define those dynamics. As a reminder, these are the five dynamics cases we will examine:\n", "\n", "1) Keplerian orbit\n", "2) Keplerian orbit + third body perturbations\n", "3) Keplerian orbit + third body perturbations + cannonball solar radiation pressure model\n", "4) Keplerian orbit + third body perturbations + N-plate solar radiation pressure model\n", "5) Keplerian orbit + third body perturbations + N-plate solar radiation pressure model + spherical harmonics gravity\n", "\n", "We'll define each of these cases using `ForceModelTranslation` and then feed it to `prop_and_plot_fm` to compare the two propagators. We're also saving the percent differences computed by `prop_and_plot_fm` to compare the two integrators for every force model later." ] }, { "cell_type": "code", "execution_count": null, "id": "10", "metadata": {}, "outputs": [], "source": [ "## create and plot increasingly complex dynamics\n", "# two body orbit\n", "twobp = scb.ForceModelTranslation(primary_body = sc)\n", "twobp_errs = prop_and_plot_fm(twobp, f'a = {a} km, i = {np.rad2deg(i):.2f}$^\\circ$\\n'\n", " 'Earth Orbit')" ] }, { "cell_type": "code", "execution_count": null, "id": "11", "metadata": {}, "outputs": [], "source": [ "# same as above with third bodies added\n", "nBP = scb.ForceModelTranslation(primary_body = sc, \n", " third_bodies = ['MERCURY', 'VENUS'])\n", "nBP_errs = prop_and_plot_fm(nBP, f'a = {a} km, i = {np.rad2deg(i):.2f}$^\\circ$\\n'\n", " 'Earth Orbit | 3rd Bodies: Mercury & Venus')" ] }, { "cell_type": "code", "execution_count": null, "id": "12", "metadata": {}, "outputs": [], "source": [ "# all above + cannonball SRP\n", "cannon = scb.ForceModelTranslation(primary_body = sc,\n", " third_bodies = ['MERCURY', 'VENUS'],\n", " cannonball_SRP = True)\n", "cannon_errs = prop_and_plot_fm(cannon, f'a = {a} km, i = {np.rad2deg(i):.2f}$^\\circ$\\n'\n", " 'Earth Orbit | 3rd Bodies: Mercury & Venus | Cannonball SRP')" ] }, { "cell_type": "code", "execution_count": null, "id": "13", "metadata": {}, "outputs": [], "source": [ "# same as above with N-plate SRP instead of cannonball\n", "npl = scb.ForceModelTranslation(primary_body = sc,\n", " third_bodies = ['MERCURY', 'VENUS'],\n", " nplate_SRP = True)\n", "npl_errs = prop_and_plot_fm(npl, f'a = {a} km, i = {np.rad2deg(i):.2f}$^\\circ$\\n'\n", " 'Earth Orbit | 3rd Bodies: Mercury & Venus | N-Plate SRP')" ] }, { "cell_type": "code", "execution_count": null, "id": "14", "metadata": {}, "outputs": [], "source": [ "# all above + spherical harmonics\n", "sph = scb.ForceModelTranslation(primary_body = sc, \n", " third_bodies = ['MERCURY', 'VENUS'], \n", " nplate_SRP = True,\n", " sph_harm = True,\n", " sph_harm_order = sh_info['order'],\n", " sph_harm_cs_file = sh_info['cs_file'],\n", " sph_harm_body = sh_info['body'],\n", " sph_harm_norm_flag = sh_info['norm_flag'])\n", "sph_errs = prop_and_plot_fm(sph, f'a = {a} km, i = {np.rad2deg(i):.2f}$^\\circ$\\n'\n", " 'Earth Orbit | 3rd Bodies: Mercury & Venus | N-Plate SRP | 3rd Order SPH')" ] }, { "cell_type": "markdown", "id": "15", "metadata": {}, "source": [ "Now plot the average difference between DOP853 and IAS15 for position and velocity for all force model configurations." ] }, { "cell_type": "code", "execution_count": null, "id": "16", "metadata": {}, "outputs": [], "source": [ "## plot average percent differences\n", "# set up titles and colors\n", "cases = ['2BP + 3rd Bodies\\n& N-Plate & SPH', '2BP + 3rd Bodies\\n& N-Plate', \n", " '2BP + 3rd Bodies\\n& Cannonball', '2BP + 3rd Bodies', '2BP']\n", "ind, width = np.arange(len(cases)), 0.4\n", "cmap = plt.get_cmap('tab10')\n", "c = cmap(np.linspace(0, 1, len(cases)))\n", "\n", "# position differences\n", "fig_bar_pos, ax_bar_pos = plt.subplots(figsize = (15, 8))\n", "# rust vs dop\n", "bars = ax_bar_pos.barh(ind + width/2, [sph_errs[0], npl_errs[0], cannon_errs[0], nBP_errs[0], twobp_errs[0]], width, color = c)\n", "ax_bar_pos.bar_label(bars, fmt='%.2e', padding=5)\n", "\n", "ax_bar_pos.set(yticks = ind, yticklabels = cases)\n", "ax_bar_pos.set_xlabel('Average Position Difference [%]'), ax_bar_pos.set_ylabel('Cases')\n", "ax_bar_pos.set_xlim(0, 1.1*ax_bar_pos.get_xlim()[1]) # padding for bar labels\n", "ax_bar_pos.set_title('Average Position Differences Between IAS15 and DOP853')\n", "\n", "# velocity differences\n", "fig_bar_vel, ax_bar_vel = plt.subplots(figsize = (15, 8))\n", "# rust vs dop\n", "bars = ax_bar_vel.barh(ind + width/2, [sph_errs[1], npl_errs[1], cannon_errs[1], nBP_errs[1], twobp_errs[1]], width, color = c)\n", "ax_bar_vel.bar_label(bars, fmt='%.2e', padding=5)\n", "\n", "ax_bar_vel.set(yticks = ind, yticklabels = cases)\n", "ax_bar_vel.set_xlabel('Average Velocity Difference [%]'), ax_bar_vel.set_ylabel('Cases')\n", "ax_bar_vel.set_xlim(0, 1.1*ax_bar_vel.get_xlim()[1]) # padding for bar labels\n", "ax_bar_vel.set_title('Average Velocity Differences Between IAS15 and DOP853')" ] }, { "cell_type": "markdown", "id": "17", "metadata": {}, "source": [ "## Conclusion\n", "FILL OUT" ] } ], "metadata": { "kernelspec": { "display_name": ".venv", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.9" } }, "nbformat": 4, "nbformat_minor": 5 }