From MPC Orbit to N-Body Simulation with Rebound¶
This tutorial demonstrates how to take an orbit from the MPC's Orbits API and use it in an N-body simulation.¶
We will:
- Fetch the orbit of the TNO Sedna from the MPC Orbits API
- Parse it with the
mpc_orbPython package to access Cartesian elements and the covariance matrix - Sample "orbital clones" from the covariance matrix to represent orbital uncertainty
- Set up a Rebound N-body simulation with the Sun, planets, and our asteroid clones
- Integrate forward for 10 million years and plot the orbital evolution
This workflow connects MPC data products directly to the widely-used Rebound dynamical simulation package.
Further information:
Environment Setup¶
This notebook requires rebound, mpc-orb, astropy, matplotlib, numpy, tqdm, and requests. We recommend creating a dedicated conda environment and registering it as a Jupyter kernel so that all dependencies are isolated and reproducible.
Run the following commands once in a terminal before opening this notebook:
# Create and populate the environment
conda create -n orbit-tutorial python=3.11 -y
conda activate orbit-tutorial
pip install rebound mpc-orb astropy matplotlib numpy tqdm requests ipykernel
# Register the environment as a Jupyter kernel
python -m ipykernel install --user --name orbit-tutorial --display-name "Python (orbit-tutorial)"
Then, in Jupyter, select Kernel > Change Kernel > Python (orbit-tutorial) before running the cells below.
Install and Import Packages¶
We need rebound for the N-body simulation, mpc-orb for parsing MPC orbit data, and astropy for time conversions.
In [2]:
Copied!
import requests
import numpy as np
import matplotlib.pyplot as plt
from pprint import pprint
import rebound
from mpc_orb.parse import MPCORB
from astropy.time import Time
from tqdm import tqdm
import requests
import numpy as np
import matplotlib.pyplot as plt
from pprint import pprint
import rebound
from mpc_orb.parse import MPCORB
from astropy.time import Time
from tqdm import tqdm
Step 1: Fetch and Parse the Orbit¶
We query the MPC Orbits API for Sedna, a well-known TNO, and parse the result with mpc_orb.
The MPCORB class provides structured access to orbital elements, uncertainties, covariance matrices, and metadata.
In [4]:
Copied!
# Fetch the orbit of Sedna from the MPC
object_name = "Sedna"
response = requests.get(
"https://data.minorplanetcenter.net/api/get-orb",
json={"desig": object_name}
)
response.raise_for_status()
mpc_orb_dict = response.json()[0]['mpc_orb'][0]
# Parse with mpc_orb
M = MPCORB(mpc_orb_dict)
# Display basic information
print(f"Object Designation Data:")
pprint(M.designation_data, indent=2, depth=2 )
print(f"Epoch (MJD): {M.epoch_data['epoch']}")
print(f"Time system: {M.epoch_data['timesystem']}")
print()
# Display Cartesian elements with uncertainties
print("Cartesian orbital elements (AU, AU/day):")
print(f"{'Name':<6} {'Value':>20} {'Uncertainty':>16}")
print("-" * 44)
for name, val, unc in zip(
M.CAR.coefficient_names,
M.CAR.coefficient_values,
M.CAR.coefficient_uncertainties
):
print(f"{name:<6} {val:>20.15f} {unc:>16.6e}")
# Fetch the orbit of Sedna from the MPC
object_name = "Sedna"
response = requests.get(
"https://data.minorplanetcenter.net/api/get-orb",
json={"desig": object_name}
)
response.raise_for_status()
mpc_orb_dict = response.json()[0]['mpc_orb'][0]
# Parse with mpc_orb
M = MPCORB(mpc_orb_dict)
# Display basic information
print(f"Object Designation Data:")
pprint(M.designation_data, indent=2, depth=2 )
print(f"Epoch (MJD): {M.epoch_data['epoch']}")
print(f"Time system: {M.epoch_data['timesystem']}")
print()
# Display Cartesian elements with uncertainties
print("Cartesian orbital elements (AU, AU/day):")
print(f"{'Name':<6} {'Value':>20} {'Uncertainty':>16}")
print("-" * 44)
for name, val, unc in zip(
M.CAR.coefficient_names,
M.CAR.coefficient_values,
M.CAR.coefficient_uncertainties
):
print(f"{name:<6} {val:>20.15f} {unc:>16.6e}")
Object Designation Data:
{ 'citation': '(90377) Sedna = 2003 VB12<br><br>Sedna is the Inuit goddess of '
'the sea and the mother of all sea creatures. She rewards the '
'people of the land with food from the sea. Without her '
'blessing, hunts fail and the people starve. She is thus one of '
'the most important figures in Inuit legend. ',
'designation_count': 1,
'iau_name': '',
'name': 'Sedna',
'orbfit_name': '90377',
'packed_primary_provisional_designation': 'K03V12B',
'packed_secondary_provisional_designations': [],
'permid': '90377',
'unpacked_primary_provisional_designation': '2003 VB12',
'unpacked_secondary_provisional_designations': []}
Epoch (MJD): 61000.0
Time system: TDT
Cartesian orbital elements (AU, AU/day):
Name Value Uncertainty
--------------------------------------------
x 39.500123133868897 1.351020e-04
y 71.045083701809901 2.399870e-04
z -17.059777137478701 5.869980e-05
vx -0.002486517552993 1.578040e-08
vy 0.000602970985653 2.286030e-08
vz 0.000201474545227 5.740050e-09
Step 2: Examine the Covariance Matrix¶
The mpc_orb package reconstructs the full symmetric covariance matrix from the stored lower-triangular elements. This matrix encodes the correlations between all orbital parameters.
We verify the matrix is positive-definite by checking that all eigenvalues are positive.
In [5]:
Copied!
# Full covariance matrix (already assembled by mpc_orb)
cov_full = M.CAR.covariance_array
print(f"Covariance matrix shape: {cov_full.shape}")
print()
print("Covariance matrix:")
np.set_printoptions(precision=4, linewidth=120)
print(cov_full)
print()
# Compute eigenvalues to verify positive-definiteness
eigenvalues = np.linalg.eigvalsh(cov_full)
print("Eigenvalues of covariance matrix:")
for i, ev in enumerate(eigenvalues):
print(f" lambda_{i} = {ev:.6e}")
# Sanity check: all eigenvalues must be positive for a valid covariance matrix
assert np.all(eigenvalues > 0), "Covariance matrix is not positive-definite!"
print("\nAll eigenvalues positive: covariance matrix is valid.")
# Full covariance matrix (already assembled by mpc_orb)
cov_full = M.CAR.covariance_array
print(f"Covariance matrix shape: {cov_full.shape}")
print()
print("Covariance matrix:")
np.set_printoptions(precision=4, linewidth=120)
print(cov_full)
print()
# Compute eigenvalues to verify positive-definiteness
eigenvalues = np.linalg.eigvalsh(cov_full)
print("Eigenvalues of covariance matrix:")
for i, ev in enumerate(eigenvalues):
print(f" lambda_{i} = {ev:.6e}")
# Sanity check: all eigenvalues must be positive for a valid covariance matrix
assert np.all(eigenvalues > 0), "Covariance matrix is not positive-definite!"
print("\nAll eigenvalues positive: covariance matrix is valid.")
Covariance matrix shape: (6, 6) Covariance matrix: [[ 1.8253e-08 3.2373e-08 -7.8984e-09 -2.9356e-14 1.7614e-12 -3.1995e-13] [ 3.2373e-08 5.7594e-08 -1.4042e-08 -2.0778e-13 2.9579e-12 -5.2034e-13] [-7.8984e-09 -1.4042e-08 3.4457e-09 4.3914e-14 -7.2821e-13 1.3413e-13] [-2.9356e-14 -2.0778e-13 4.3914e-14 2.4902e-16 2.9151e-16 -7.9405e-17] [ 1.7614e-12 2.9579e-12 -7.2821e-13 2.9151e-16 5.2260e-16 -1.2586e-16] [-3.1995e-13 -5.2034e-13 1.3413e-13 -7.9405e-17 -1.2586e-16 3.2948e-17]] Eigenvalues of covariance matrix: lambda_0 = 1.491403e-19 lambda_1 = 3.993476e-19 lambda_2 = 3.197831e-16 lambda_3 = 2.044994e-11 lambda_4 = 4.240028e-11 lambda_5 = 7.922912e-08 All eigenvalues positive: covariance matrix is valid.
Step 3: Sample Orbital Clones¶
To represent the orbital uncertainty, we draw samples from a multivariate normal distribution where:
- The mean is the nominal (best-fit) orbit
- The covariance defines the uncertainty ellipsoid in 6D phase space
Though MPC orbit fits currently do not explicitly reference a prior (and therefore implicitly use flat priors), the mean and covariance matrix can be interpreted as the Gaussian approximation of the posterior distribution for the orbital elements. Thus samples drawn from this distribution ("clones") each represent a statistically plausible orbit, and by propagating a distribution of samples, we can quantify how the orbital uncertainty grows over time. Since more clones will be drawn from near the mean than towards the edges, the on-sky density of clones at any time is proportional to the likelihood of finding Sedna in that region
In [6]:
Copied!
# Use only the 6 position/velocity components (x, y, z, vx, vy, vz)
mean = np.array(M.CAR.coefficient_values[:6])
cov = cov_full[:6, :6]
# Sample 100 orbital clones
n_clones = 100
np.random.seed(42) # For reproducibility
clones = np.random.multivariate_normal(mean, cov, n_clones)
# Display the nominal orbit and clones
names = M.CAR.coefficient_names[:6]
print(f"{'':>10}", end="")
for name in names:
print(f"{name:>20}", end="")
print()
print("-" * 130)
print(f"{'Nominal':>10}", end="")
for v in mean:
print(f"{v:>20.15f}", end="")
print()
for i, clone in enumerate(clones):
print(f"{'Clone ' + str(i):>10}", end="")
for v in clone:
print(f"{v:>20.15f}", end="")
print()
# Use only the 6 position/velocity components (x, y, z, vx, vy, vz)
mean = np.array(M.CAR.coefficient_values[:6])
cov = cov_full[:6, :6]
# Sample 100 orbital clones
n_clones = 100
np.random.seed(42) # For reproducibility
clones = np.random.multivariate_normal(mean, cov, n_clones)
# Display the nominal orbit and clones
names = M.CAR.coefficient_names[:6]
print(f"{'':>10}", end="")
for name in names:
print(f"{name:>20}", end="")
print()
print("-" * 130)
print(f"{'Nominal':>10}", end="")
for v in mean:
print(f"{v:>20.15f}", end="")
print()
for i, clone in enumerate(clones):
print(f"{'Clone ' + str(i):>10}", end="")
for v in clone:
print(f"{v:>20.15f}", end="")
print()
x y z vx vy vz ---------------------------------------------------------------------------------------------------------------------------------- Nominal 39.500123133868897 71.045083701809901 -17.059777137478701 -0.002486517552993 0.000602970985653 0.000201474545227 Clone 0 39.500055061295633 71.044964381095539 -17.059750948143776 -0.002486535047046 0.000602941585839 0.000201481116801 Clone 1 39.499914508397133 71.044702727682989 -17.059682511127708 -0.002486513307365 0.000602954044219 0.000201477927934 Clone 2 39.500080224057399 71.045033195739634 -17.059755632371782 -0.002486533744880 0.000602951551896 0.000201482429683 Clone 3 39.500237125194928 71.045304786723136 -17.059836978601542 -0.002486532953466 0.000602966821794 0.000201474994355 Clone 4 39.500197670872666 71.045214986285615 -17.059803895405572 -0.002486520646722 0.000602974549439 0.000201475819537 Clone 5 39.500214896216356 71.045222241502870 -17.059812021116610 -0.002486485363607 0.000603017510058 0.000201461971409 Clone 6 39.500084281640135 71.045040939328331 -17.059759331034005 -0.002486542131393 0.000602940195772 0.000201483879961 Clone 7 39.500137576166182 71.045113694210755 -17.059777413454778 -0.002486513166770 0.000602978940388 0.000201474759916 Clone 8 39.500066601849277 71.045006541248910 -17.059758720619904 -0.002486533872166 0.000602948997984 0.000201480434995 Clone 9 39.499989575520566 71.044834083240062 -17.059712940137508 -0.002486502209702 0.000602974589341 0.000201473196067 Clone 10 39.500187165067075 71.045200238028698 -17.059800320290570 -0.002486506811970 0.000602991151487 0.000201470738523 Clone 11 39.500138429940790 71.045097491565869 -17.059782801941431 -0.002486498573141 0.000602993535967 0.000201468053277 Clone 12 39.500137837084210 71.045089631048967 -17.059767416347974 -0.002486507678467 0.000602981229089 0.000201474484080 Clone 13 39.500099521043197 71.045068162265039 -17.059771091507724 -0.002486544448809 0.000602938646669 0.000201482599109 Clone 14 39.500229073597787 71.045278501511490 -17.059828583743936 -0.002486527505402 0.000602969700806 0.000201475177158 Clone 15 39.500115799826979 71.045057953687817 -17.059768206230952 -0.002486503094408 0.000602987594625 0.000201470768272 Clone 16 39.500084646394612 71.045011813044155 -17.059759788390139 -0.002486512152783 0.000602974418853 0.000201474107541 Clone 17 39.500164889129316 71.045168619964159 -17.059796603979169 -0.002486531098794 0.000602958662235 0.000201477328375 Clone 18 39.500088653033920 71.045023806215880 -17.059753588442707 -0.002486517001881 0.000602967292190 0.000201477714252 Clone 19 39.500150827386747 71.045128937841682 -17.059788198438920 -0.002486501439159 0.000602993244292 0.000201468279530 Clone 20 39.500010668222203 71.044895537124091 -17.059737175053872 -0.002486512294862 0.000602967830673 0.000201473023063 Clone 21 39.500253579484379 71.045323100298489 -17.059835641526170 -0.002486519817852 0.000602983362092 0.000201473328440 Clone 22 39.500269558695479 71.045337919736042 -17.059835180720942 -0.002486529200845 0.000602969330108 0.000201477493787 Clone 23 39.500006237513418 71.044892182576604 -17.059730709437464 -0.002486545389504 0.000602926775169 0.000201486279248 Clone 24 39.500092958358202 71.045019943360487 -17.059756339276944 -0.002486493982552 0.000602996603813 0.000201468536035 Clone 25 39.500091543648026 71.045023088791652 -17.059759416611730 -0.002486515852033 0.000602969398087 0.000201475135706 Clone 26 39.499874423538657 71.044635520878089 -17.059662589653527 -0.002486517368534 0.000602945160194 0.000201481309579 Clone 27 39.499961720610585 71.044807437052469 -17.059713704126747 -0.002486530115622 0.000602939157370 0.000201480438524 Clone 28 39.500152296396557 71.045145746037321 -17.059787679655265 -0.002486517774494 0.000602975836446 0.000201474804243 Clone 29 39.500090487031052 71.045014677366382 -17.059760876092298 -0.002486522208980 0.000602957593479 0.000201478054797 Clone 30 39.500034198452788 71.044937208811433 -17.059735907651973 -0.002486531937135 0.000602945477066 0.000201482195694 Clone 31 39.500059153220676 71.044971109822399 -17.059745703798260 -0.002486501824853 0.000602985515898 0.000201471617564 Clone 32 39.500087084000043 71.045036117343400 -17.059765557593529 -0.002486536169860 0.000602946791660 0.000201481284212 Clone 33 39.500108642247959 71.045073038227372 -17.059775483005588 -0.002486536377840 0.000602946898819 0.000201480067316 Clone 34 39.500303568538776 71.045416812021614 -17.059860193805608 -0.002486533771593 0.000602967775318 0.000201476576911 Clone 35 39.500052183054599 71.044942284302408 -17.059747770936859 -0.002486510631948 0.000602969238953 0.000201473353817 Clone 36 39.500226282647695 71.045270324633364 -17.059820262139912 -0.002486521888545 0.000602976288608 0.000201472909981 Clone 37 39.500021500210387 71.044924537268258 -17.059735122334594 -0.002486547697349 0.000602926432665 0.000201486193263 Clone 38 39.500223822349518 71.045253835563571 -17.059815675127126 -0.002486512481844 0.000602986111627 0.000201471925558 Clone 39 39.499838093872846 71.044568990085835 -17.059642596405666 -0.002486510178805 0.000602950509887 0.000201480354217 Clone 40 39.500229270721206 71.045273801693710 -17.059825776187381 -0.002486528924532 0.000602967506715 0.000201476438678 Clone 41 39.500182875047287 71.045198200639035 -17.059812841828894 -0.002486529935602 0.000602962341984 0.000201476313386 Clone 42 39.499843115810158 71.044572516795370 -17.059646048840843 -0.002486498683823 0.000602964991117 0.000201474536880 Clone 43 39.500067975175597 71.044975557679393 -17.059746944950881 -0.002486508225577 0.000602977708636 0.000201475117602 Clone 44 39.500149525179168 71.045146834736158 -17.059799325801226 -0.002486517501864 0.000602977415206 0.000201471813463 Clone 45 39.499919990740445 71.044741376697885 -17.059698133433216 -0.002486533608207 0.000602933879039 0.000201482377885 Clone 46 39.500092823935056 71.045037772731632 -17.059765884145705 -0.002486520179434 0.000602965868658 0.000201475731358 Clone 47 39.499901222179240 71.044705171781601 -17.059693920497374 -0.002486510453412 0.000602961552267 0.000201473333748 Clone 48 39.500081739129662 71.045018425205228 -17.059759862769095 -0.002486519269837 0.000602966290507 0.000201476319259 Clone 49 39.500070669615639 71.044999439183641 -17.059760321810057 -0.002486528448636 0.000602952570694 0.000201477671298 Clone 50 39.500231570499743 71.045283752435694 -17.059829049072324 -0.002486531193156 0.000602965247585 0.000201476203888 Clone 51 39.499947105331742 71.044778502480355 -17.059704859665299 -0.002486521214790 0.000602949802765 0.000201478511485 Clone 52 39.500015864903212 71.044881919833344 -17.059734475115739 -0.002486507789087 0.000602970486632 0.000201471677567 Clone 53 39.500078601228232 71.045006240601211 -17.059758611434969 -0.002486524631658 0.000602956718887 0.000201478612387 Clone 54 39.500251593021765 71.045327942770584 -17.059841324916558 -0.002486540690586 0.000602955991735 0.000201478137275 Clone 55 39.500119654272908 71.045089413358582 -17.059778323363869 -0.002486520846032 0.000602968627114 0.000201474654526 Clone 56 39.500232580125015 71.045282448927551 -17.059827330139061 -0.002486516067388 0.000602985025182 0.000201471875871 Clone 57 39.500087357394463 71.045026932825564 -17.059760791526664 -0.002486526351797 0.000602958830938 0.000201478849934 Clone 58 39.500218781728520 71.045256507350572 -17.059820596761799 -0.002486536063629 0.000602956261861 0.000201478283143 Clone 59 39.500120003868382 71.045091453029201 -17.059778318218665 -0.002486526469715 0.000602962226914 0.000201476791196 Clone 60 39.500061791899768 71.044954322428609 -17.059746029082344 -0.002486503453801 0.000602978838397 0.000201471266516 Clone 61 39.500092921109513 71.045029802627411 -17.059764450762998 -0.002486508959911 0.000602979028681 0.000201472136928 Clone 62 39.499931920001558 71.044730567707319 -17.059701576400325 -0.002486497249038 0.000602975192594 0.000201468728469 Clone 63 39.499823261633672 71.044561515722066 -17.059645369313284 -0.002486519537404 0.000602942338518 0.000201482138652 Clone 64 39.500226332551726 71.045265093103552 -17.059823060516838 -0.002486536165120 0.000602955964661 0.000201478182616 Clone 65 39.500247688708235 71.045298883411007 -17.059823812680282 -0.002486531534219 0.000602963224863 0.000201478133322 Clone 66 39.500362120326187 71.045490609689125 -17.059876703675915 -0.002486516356342 0.000602993346055 0.000201471257097 Clone 67 39.500122859275308 71.045082686997404 -17.059774833311078 -0.002486523715750 0.000602963269422 0.000201477711834 Clone 68 39.500109560592264 71.045052587499555 -17.059773166328977 -0.002486499522719 0.000602991876406 0.000201469139393 Clone 69 39.500073447493151 71.045004974803646 -17.059764657170668 -0.002486527284006 0.000602954637989 0.000201476125777 Clone 70 39.499855049205962 71.044583858579983 -17.059655153269681 -0.002486505216186 0.000602954252835 0.000201476432520 Clone 71 39.500256893693894 71.045312156502050 -17.059838195320651 -0.002486492575880 0.000603016947705 0.000201462671667 Clone 72 39.500163027343575 71.045144770684075 -17.059799443162319 -0.002486510448877 0.000602982059503 0.000201468875507 Clone 73 39.500352601571137 71.045492294338018 -17.059878547962430 -0.002486519457938 0.000602992441202 0.000201471569054 Clone 74 39.500302888027839 71.045394290112100 -17.059855013040178 -0.002486501164982 0.000603009391566 0.000201465499816 Clone 75 39.500137391815564 71.045096579488472 -17.059776304827597 -0.002486511919108 0.000602977790447 0.000201473993536 Clone 76 39.500131895263486 71.045113143094085 -17.059781108396795 -0.002486516766982 0.000602975189619 0.000201473194335 Clone 77 39.500218842724962 71.045251032008679 -17.059817565764362 -0.002486512882709 0.000602985353560 0.000201471312152 Clone 78 39.500191567073045 71.045213064802923 -17.059807058808218 -0.002486500217422 0.000603002498767 0.000201468407536 Clone 79 39.499899464143041 71.044689253405636 -17.059683439565791 -0.002486521842533 0.000602941641169 0.000201478101418 Clone 80 39.500135556230241 71.045118820398187 -17.059777663531378 -0.002486531343201 0.000602958198704 0.000201480371158 Clone 81 39.500203634251946 71.045240808223284 -17.059822605244786 -0.002486539731426 0.000602951823563 0.000201478711103 Clone 82 39.500107229088343 71.045066534025096 -17.059779472751121 -0.002486533659159 0.000602951065985 0.000201478733850 Clone 83 39.500233117806211 71.045297387122105 -17.059832691801972 -0.002486554653527 0.000602937143418 0.000201484483917 Clone 84 39.500208405588914 71.045241890342879 -17.059812674480021 -0.002486514794149 0.000602985147223 0.000201472178221 Clone 85 39.500086429393228 71.045019168920518 -17.059760257714892 -0.002486508222412 0.000602979722183 0.000201472712041 Clone 86 39.500049921113728 71.044966505721575 -17.059748402433478 -0.002486536675328 0.000602941751713 0.000201482864964 Clone 87 39.500216103438206 71.045241567806912 -17.059812462426706 -0.002486492466728 0.000603012602937 0.000201466131106 Clone 88 39.500082695753846 71.045023682377646 -17.059764893982170 -0.002486510903024 0.000602979186402 0.000201471486663 Clone 89 39.500211706555859 71.045240719104370 -17.059811479406257 -0.002486513329595 0.000602985233400 0.000201471783289 Clone 90 39.500003708381087 71.044886277902023 -17.059730733383322 -0.002486546063186 0.000602926250733 0.000201486446156 Clone 91 39.500043947734625 71.044945543990494 -17.059745031330973 -0.002486513233234 0.000602969454438 0.000201473886305 Clone 92 39.499966827775665 71.044802376401179 -17.059710238830803 -0.002486509641723 0.000602965575377 0.000201474547119 Clone 93 39.500067388271049 71.044990126681384 -17.059755393271995 -0.002486544071445 0.000602931286280 0.000201483688084 Clone 94 39.499959421661927 71.044798519983999 -17.059697875101726 -0.002486511017121 0.000602965907469 0.000201478163523 Clone 95 39.500130064897327 71.045073693184435 -17.059777255050147 -0.002486496151960 0.000602995812221 0.000201466519799 Clone 96 39.500096260957861 71.045026012026696 -17.059756692862532 -0.002486520305872 0.000602962568781 0.000201477769355 Clone 97 39.500050571492054 71.044924513957824 -17.059740564494565 -0.002486493448990 0.000602989502242 0.000201468531038 Clone 98 39.500006156922083 71.044887096827125 -17.059728979721935 -0.002486521798086 0.000602955889728 0.000201477464969 Clone 99 39.499980282879953 71.044836577780387 -17.059715289972672 -0.002486512403706 0.000602965165048 0.000201475720873
Step 4: Set Up the Rebound Simulation¶
We create a Rebound simulation at the MPC orbit epoch containing:
- The Sun and 4 planets (loaded from JPL Horizons)
- The nominal orbit of Sedna as a massless test particle
- 100 orbital clones as additional test particles
Unit handling: MPC Cartesian elements use AU for positions and AU/day for velocities. We set sim.units = ('AU', 'day', 'Msun') so that Rebound and Horizons use the same units. This ensures we can directly insert MPC state vectors without any conversion.
We use the WHFast symplectic integrator, which is efficient for long-term integrations.
The timestep is set to 1/20th of Jupiters's orbital period for accuracy.
NB: It is possible that you may experience the following error when adding planets to rebound (which uses the NASA Horizons system to extract the planetary positions): RuntimeError: An error occured while accessing NASA HORIZONS. If this is a SSL certificate issue, you can try disabling the certificate verification by setting rebound.horizons.SSL_CONTEXT = 'unverified'.
- In this scenario, the simplest approach is to just retry the evaluation, and this typically will succeed.
In [7]:
Copied!
# Convert MPC epoch (MJD) to a date string for Horizons
epoch_mjd = M.epoch_data['epoch']
epoch_date = Time(epoch_mjd, format='mjd').iso[:10] # YYYY-MM-DD
print(f"Epoch: MJD {epoch_mjd} = {epoch_date}")
# Create the simulation with AU/day/Msun units to match MPC data
sim = rebound.Simulation()
sim.units = ('AU', 'day', 'Msun')
# Add the Sun and planets from Horizons at the MPC epoch
bodies = ["Sun", "Jupiter", "Saturn", "Uranus", "Neptune"] # "Mercury", "Venus", "Earth", "Mars",
for body in bodies:
sim.add(body, date=epoch_date)
print(f"Added {body}")
# Add the nominal orbit as a test particle (mass=0)
sim.add(x=mean[0], y=mean[1], z=mean[2],
vx=mean[3], vy=mean[4], vz=mean[5])
print(f"Added nominal {object_name} orbit")
# Add the clones as test particles
for i, clone in enumerate(clones):
sim.add(x=clone[0], y=clone[1], z=clone[2],
vx=clone[3], vy=clone[4], vz=clone[5])
print(f"Added {n_clones} orbital clones")
# Configure the integrator
sim.integrator = "whfast"
sim.dt = sim.particles[1].P / 20.0 # 1/20 of inner-most planet's orbital period
sim.move_to_com()
# Verify the nominal orbit looks correct
orb = sim.particles[9].orbit(primary=sim.particles[0])
print(f"\nNominal {object_name} orbit check: a={orb.a:.4f} AU, e={orb.e:.4f}, inc={np.degrees(orb.inc):.2f} deg")
print(f"Total particles: {sim.N} (Sun + 8 planets + 1 nominal + {n_clones} clones)")
print(f"Integrator: {sim.integrator}")
print(f"Timestep: {sim.dt:.4f} days ({sim.dt/365.25:.6f} years)")
# Convert MPC epoch (MJD) to a date string for Horizons
epoch_mjd = M.epoch_data['epoch']
epoch_date = Time(epoch_mjd, format='mjd').iso[:10] # YYYY-MM-DD
print(f"Epoch: MJD {epoch_mjd} = {epoch_date}")
# Create the simulation with AU/day/Msun units to match MPC data
sim = rebound.Simulation()
sim.units = ('AU', 'day', 'Msun')
# Add the Sun and planets from Horizons at the MPC epoch
bodies = ["Sun", "Jupiter", "Saturn", "Uranus", "Neptune"] # "Mercury", "Venus", "Earth", "Mars",
for body in bodies:
sim.add(body, date=epoch_date)
print(f"Added {body}")
# Add the nominal orbit as a test particle (mass=0)
sim.add(x=mean[0], y=mean[1], z=mean[2],
vx=mean[3], vy=mean[4], vz=mean[5])
print(f"Added nominal {object_name} orbit")
# Add the clones as test particles
for i, clone in enumerate(clones):
sim.add(x=clone[0], y=clone[1], z=clone[2],
vx=clone[3], vy=clone[4], vz=clone[5])
print(f"Added {n_clones} orbital clones")
# Configure the integrator
sim.integrator = "whfast"
sim.dt = sim.particles[1].P / 20.0 # 1/20 of inner-most planet's orbital period
sim.move_to_com()
# Verify the nominal orbit looks correct
orb = sim.particles[9].orbit(primary=sim.particles[0])
print(f"\nNominal {object_name} orbit check: a={orb.a:.4f} AU, e={orb.e:.4f}, inc={np.degrees(orb.inc):.2f} deg")
print(f"Total particles: {sim.N} (Sun + 8 planets + 1 nominal + {n_clones} clones)")
print(f"Integrator: {sim.integrator}")
print(f"Timestep: {sim.dt:.4f} days ({sim.dt/365.25:.6f} years)")
Epoch: MJD 61000.0 = 2025-11-21 Searching NASA Horizons for 'Sun'... Found: Sun (10) Added Sun Searching NASA Horizons for 'Jupiter'... Found: Jupiter Barycenter (5) (chosen from query 'Jupiter') Added Jupiter Searching NASA Horizons for 'Saturn'... Found: Saturn Barycenter (6) (chosen from query 'Saturn') Added Saturn Searching NASA Horizons for 'Uranus'... Found: Uranus Barycenter (7) (chosen from query 'Uranus') Added Uranus Searching NASA Horizons for 'Neptune'... Found: Neptune Barycenter (8) (chosen from query 'Neptune') Added Neptune Added nominal Sedna orbit Added 100 orbital clones Nominal Sedna orbit check: a=590.4016 AU, e=0.8709, inc=11.92 deg Total particles: 106 (Sun + 8 planets + 1 nominal + 100 clones) Integrator: whfast Timestep: 216.6078 days (0.593040 years)
Step 5: Run the Simulation¶
We integrate the system forward for 10 million years, recording the semi-major axis, eccentricity, and inclination of the nominal orbit and each clone at 1,000 evenly-spaced output times.
Note: This integration may take a few minutes depending on your machine.
In [8]:
Copied!
# Integration parameters
t_final = 10.0e6 * 365.25 # 10 million years in days
n_outputs = 1000
times = np.linspace(0, t_final, n_outputs)
# Indices of test particles: nominal is after Sun + planets, clones after that ...
i_nominal = len(bodies)
i_first_clone = i_nominal + 1
i_clones = list(range(i_first_clone, i_first_clone + n_clones))
all_test = [i_nominal] + i_clones
n_test = len(all_test)
# Storage arrays
a_arr = np.zeros((n_outputs, n_test))
e_arr = np.zeros((n_outputs, n_test))
inc_arr = np.zeros((n_outputs, n_test))
# Integrate
for k, t in enumerate(tqdm(times)):
sim.integrate(t)
for j, idx in enumerate(all_test):
orb = sim.particles[idx].orbit(primary=sim.particles[0])
a_arr[k, j] = orb.a
e_arr[k, j] = orb.e
inc_arr[k, j] = np.degrees(orb.inc)
print("Integration complete.")
# Convert times to years for plotting
times_yr = times / 365.25
# Integration parameters
t_final = 10.0e6 * 365.25 # 10 million years in days
n_outputs = 1000
times = np.linspace(0, t_final, n_outputs)
# Indices of test particles: nominal is after Sun + planets, clones after that ...
i_nominal = len(bodies)
i_first_clone = i_nominal + 1
i_clones = list(range(i_first_clone, i_first_clone + n_clones))
all_test = [i_nominal] + i_clones
n_test = len(all_test)
# Storage arrays
a_arr = np.zeros((n_outputs, n_test))
e_arr = np.zeros((n_outputs, n_test))
inc_arr = np.zeros((n_outputs, n_test))
# Integrate
for k, t in enumerate(tqdm(times)):
sim.integrate(t)
for j, idx in enumerate(all_test):
orb = sim.particles[idx].orbit(primary=sim.particles[0])
a_arr[k, j] = orb.a
e_arr[k, j] = orb.e
inc_arr[k, j] = np.degrees(orb.inc)
print("Integration complete.")
# Convert times to years for plotting
times_yr = times / 365.25
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1000/1000 [10:27<00:00, 1.59it/s]
Integration complete.
Step 6: Plot the Results¶
We have saved the data into arrays, a_arr, e_arr & inc_arr that have shape (1,000, n_clones+1), such that we have 1,000 data points for the nominal orbit in, e.g., a_arr[:,0], and then 1,000 data points for each of the clones in , e.g., a_arr[:,1:].
We use these arrays to plot the orbital evolution of the nominal orbit (red, bold) and the clones (gray, transparent), showing how the orbital uncertainty spreads over 10 million years.
In [9]:
Copied!
# Semi-major axis vs time
fig, ax = plt.subplots(figsize=(12, 5))
for j in range(1, n_test): # clones
ax.plot(times_yr, a_arr[:, j], color='gray', linewidth=0.5, alpha=0.3, label=f'Clones' if j == 1 else None )
ax.plot(times_yr, a_arr[:, 0], color='red', linewidth=1.0, label='Nominal')
ax.set_xlabel('Time (years)')
ax.set_ylabel('Semi-major axis (AU)')
ax.set_title(f'{object_name}: Semi-major axis evolution')
ax.legend()
plt.tight_layout()
plt.show()
# Semi-major axis vs time
fig, ax = plt.subplots(figsize=(12, 5))
for j in range(1, n_test): # clones
ax.plot(times_yr, a_arr[:, j], color='gray', linewidth=0.5, alpha=0.3, label=f'Clones' if j == 1 else None )
ax.plot(times_yr, a_arr[:, 0], color='red', linewidth=1.0, label='Nominal')
ax.set_xlabel('Time (years)')
ax.set_ylabel('Semi-major axis (AU)')
ax.set_title(f'{object_name}: Semi-major axis evolution')
ax.legend()
plt.tight_layout()
plt.show()
In [10]:
Copied!
# Eccentricity vs time
fig, ax = plt.subplots(figsize=(12, 5))
for j in range(1, n_test): # clones
ax.plot(times_yr, e_arr[:, j], color='gray', linewidth=0.5, alpha=0.3, label=f'Clones' if j == 1 else None )
ax.plot(times_yr, e_arr[:, 0], color='red', linewidth=1.0, label='Nominal')
ax.set_xlabel('Time (years)')
ax.set_ylabel('Eccentricity')
ax.set_title(f'{object_name}: Eccentricity evolution')
ax.legend()
plt.tight_layout()
plt.show()
# Eccentricity vs time
fig, ax = plt.subplots(figsize=(12, 5))
for j in range(1, n_test): # clones
ax.plot(times_yr, e_arr[:, j], color='gray', linewidth=0.5, alpha=0.3, label=f'Clones' if j == 1 else None )
ax.plot(times_yr, e_arr[:, 0], color='red', linewidth=1.0, label='Nominal')
ax.set_xlabel('Time (years)')
ax.set_ylabel('Eccentricity')
ax.set_title(f'{object_name}: Eccentricity evolution')
ax.legend()
plt.tight_layout()
plt.show()
In [11]:
Copied!
# Inclination vs time
fig, ax = plt.subplots(figsize=(12, 5))
for j in range(1, n_test): # clones
ax.plot(times_yr, inc_arr[:, j], color='gray', linewidth=0.5, alpha=0.3, label=f'Clones' if j == 1 else None )
ax.plot(times_yr, inc_arr[:, 0], color='red', linewidth=1.0, label='Nominal')
ax.set_xlabel('Time (years)')
ax.set_ylabel('Inclination (degrees)')
ax.set_title(f'{object_name}: Inclination evolution')
ax.legend()
plt.tight_layout()
plt.show()
# Inclination vs time
fig, ax = plt.subplots(figsize=(12, 5))
for j in range(1, n_test): # clones
ax.plot(times_yr, inc_arr[:, j], color='gray', linewidth=0.5, alpha=0.3, label=f'Clones' if j == 1 else None )
ax.plot(times_yr, inc_arr[:, 0], color='red', linewidth=1.0, label='Nominal')
ax.set_xlabel('Time (years)')
ax.set_ylabel('Inclination (degrees)')
ax.set_title(f'{object_name}: Inclination evolution')
ax.legend()
plt.tight_layout()
plt.show()
In [12]:
Copied!
# (a, e) scatter: initial vs final epoch
fig, ax = plt.subplots(figsize=(8, 6))
# All epochs
ax.scatter(a_arr[:,1:], e_arr[:,1:], color='gray', marker='.', s=20, label='Clones: (All epochs)', zorder=0)
ax.scatter(a_arr[:,0], e_arr[:,0], color='black', marker='.', s=20, label='Nominal: (All epochs)', zorder=1)
# Initial epoch (first output)
ax.scatter(a_arr[0, 1:], e_arr[0, 1:], color='steelblue', alpha=0.7, s=40, label='Clones (initial)', zorder=2)
ax.scatter(a_arr[0, 0], e_arr[0, 0], color='darkblue', marker='*',s=200, label='Nominal (initial)', zorder=3)
# Final epoch (last output)
ax.scatter(a_arr[-1, 1:], e_arr[-1, 1:], color='coral', alpha=0.7,s=40, label='Clones (final)', zorder=2)
ax.scatter(a_arr[-1, 0], e_arr[-1, 0], color='darkred', marker='*',s=200, label='Nominal (final)', zorder=3)
ax.set_xlabel('Semi-major axis (AU)')
ax.set_ylabel('Eccentricity')
ax.set_title(f'{object_name}: (a, e) evolution')
ax.legend()
plt.tight_layout()
plt.show()
# (a, e) scatter: initial vs final epoch
fig, ax = plt.subplots(figsize=(8, 6))
# All epochs
ax.scatter(a_arr[:,1:], e_arr[:,1:], color='gray', marker='.', s=20, label='Clones: (All epochs)', zorder=0)
ax.scatter(a_arr[:,0], e_arr[:,0], color='black', marker='.', s=20, label='Nominal: (All epochs)', zorder=1)
# Initial epoch (first output)
ax.scatter(a_arr[0, 1:], e_arr[0, 1:], color='steelblue', alpha=0.7, s=40, label='Clones (initial)', zorder=2)
ax.scatter(a_arr[0, 0], e_arr[0, 0], color='darkblue', marker='*',s=200, label='Nominal (initial)', zorder=3)
# Final epoch (last output)
ax.scatter(a_arr[-1, 1:], e_arr[-1, 1:], color='coral', alpha=0.7,s=40, label='Clones (final)', zorder=2)
ax.scatter(a_arr[-1, 0], e_arr[-1, 0], color='darkred', marker='*',s=200, label='Nominal (final)', zorder=3)
ax.set_xlabel('Semi-major axis (AU)')
ax.set_ylabel('Eccentricity')
ax.set_title(f'{object_name}: (a, e) evolution')
ax.legend()
plt.tight_layout()
plt.show()
Summary¶
In this tutorial we:
- Fetched the orbit of Sedna from the MPC Orbits API
- Parsed it with the
mpc_orbpackage to access Cartesian state vectors and the full covariance matrix - Sampled 100 orbital clones from the covariance matrix using multivariate normal sampling
- Built a Rebound N-body simulation with the Sun, 4 planets, and all test particles
- Integrated the system forward for 10 million years using the WHFast symplectic integrator
- Plotted the orbital evolution to visualize how uncertainty propagates over time
The spread of the clones over time shows how even small uncertainties in the initial orbit can lead to divergent trajectories over long timescales.
Further Resources¶
For questions or feedback, contact the MPC via the Jira Helpdesk.