Examples¶
Below are some examples of how the package could be used in astronomical studies. While there are more use cases, the most common are outlined below.
High-order resonances¶
If you want to find high-order resonances, you can set up the custom config and use the code below. Important: make sure that you removed old csv files in cache directory (or custom if you have changed it) before running the code below. If the files are not removed, the app will use them and you will not get correct results.
Note that you need to specify the following variables: resonant order, maximum of the primary coefficient, and the maximum of the coefficient. You can do this for both two-body and three-body resonances.
import resonances
resonances.config.set('MATRIX_3BODY_PRIMARY_MAX', 20)
resonances.config.set('MATRIX_3BODY_COEF_MAX', 20)
resonances.config.set('MATRIX_3BODY_ORDER_MAX', 20)
resonances.config.set('MATRIX_2BODY_PRIMARY_MAX', 100)
resonances.config.set('MATRIX_2BODY_COEF_MAX', 100)
resonances.config.set('MATRIX_2BODY_ORDER_MAX', 100)
sim = resonances.find([847, 1020, 8089], ['Earth', 'Mars', 'Jupiter', 'Saturn'])
sim.config.plot = 'nonzero'
sim.config.save = 'nonzero'
sim.config.plot_type = 'save'
sim.run(progress=True)
Searching NASA Horizons for 'Sun'... Found: Sun (10) Searching NASA Horizons for 'Mercury'... Found: Mercury Barycenter (199) (chosen from query 'Mercury') Searching NASA Horizons for 'Venus'... Found: Venus Barycenter (299) (chosen from query 'Venus') Searching NASA Horizons for 'Earth'... Found: Earth-Moon Barycenter (3) (chosen from query 'Earth') Searching NASA Horizons for 'Mars'... Found: Mars Barycenter (4) (chosen from query 'Mars') Searching NASA Horizons for 'Jupiter'... Found: Jupiter Barycenter (5) (chosen from query 'Jupiter') Searching NASA Horizons for 'Saturn'... Found: Saturn Barycenter (6) (chosen from query 'Saturn') Searching NASA Horizons for 'Uranus'... Found: Uranus Barycenter (7) (chosen from query 'Uranus') Searching NASA Horizons for 'Neptune'... Found: Neptune Barycenter (8) (chosen from query 'Neptune') Searching NASA Horizons for 'Pluto'... Found: Pluto Barycenter (9) (chosen from query 'Pluto')
100%|██████████| 6283/6283 [01:55<00:00, 54.22it/s]
Secular resonances with custom config¶
For secular resonances, it is quite often that you need different configuration values for frequencies, integration time, integration parameters, etc. Furthermore, you might want to specify the concrete resonances you want to check. For this, you can set up the custom config and use the code below.
import resonances
import numpy as np
config = {
'plot': 'all',
'save': 'all',
'integration_years': 1000000, # 1,000,000 years integration
'oscillations_cutoff': 0.00005, # 20,000 yrs
'libration_period_min': 10000,
'libration_period_critical': 40000,
'periodogram_frequency_min': 0.000001, # 1 Myr
'periodogram_frequency_max': 0.0002, # 5,000 yrs
'image_type': 'png',
'plot_type': 'both',
'dt': 1.0,
'integrator': 'saba(10,6,4)',
}
z1_asteroids = [847, 363, 1020]
z2_asteroids = [8089]
z1_resonance = resonances.GeneralSecularResonance(formula='g-g6')
print(f"Testing z1 secular resonance: {z1_resonance.to_s()}")
z2_resonance = resonances.GeneralSecularResonance(formula='2(g-g6)+(s-s6)')
print(f"Testing z2 secular resonance: {z2_resonance.to_s()}")
# Create simulation
sim = resonances.Simulation(
name="test_z1_secular",
tmax=int(2000000 * 2 * np.pi),
**config,
)
sim.create_solar_system()
# Add asteroids to check for z1 resonance
for asteroid in z1_asteroids:
sim.add_body(asteroid, z1_resonance, name=f"{asteroid}")
print(f"Added asteroid {asteroid} to check for z1 resonance")
# Add asteroids to check for z2 resonance
for asteroid in z2_asteroids:
sim.add_body(asteroid, z2_resonance, name=f"{asteroid}")
print(f"Added asteroid {asteroid} to check for z2 resonance")
# Run the simulation
print("\nRunning simulation...")
sim.run(progress=True)
# Get results
summary = sim.data_manager.get_simulation_summary(sim.bodies)
print("\n" + "=" * 60)
print("RESULTS: Secular Resonance Checks")
print("=" * 60)
# Helper function to interpret status
def interpret_status(status):
# Status meanings:
# 0 = circulation (not in resonance)
# 1 = transient
# 2 = libration (in resonance)
# -1, -2 = controversial (transient or pure)
if abs(status) == 2:
return "✓ YES - Strong libration detected (in resonance)"
elif abs(status) == 0:
return "✗ NO - Circulation detected (not in resonance)"
elif status == 1:
return "? TRANSIENT"
else:
return f"? CONTROVERSIAL (status={status})"
# Display z1 results
print("\nz1 resonance (g-g6):")
for asteroid in z1_asteroids:
status = summary.loc[(summary['name'] == str(asteroid)) & (summary['resonance'] == 'g-g6'), 'status'].iloc[0]
result = interpret_status(status)
print(f" Asteroid {asteroid:4d}: {result}")
# Display z2 results
print("\nz2 resonance (2(g-g6)+(s-s6)):")
for asteroid in z2_asteroids:
status = summary.loc[(summary['name'] == str(asteroid)) & (summary['resonance'] == '2(g-g6)+(s-s6)'), 'status'].iloc[0]
result = interpret_status(status)
print(f" Asteroid {asteroid:4d}: {result}")
Testing z1 secular resonance: g-g6 Testing z2 secular resonance: 2(g-g6)+(s-s6) Added asteroid 847 to check for z1 resonance Added asteroid 363 to check for z1 resonance Added asteroid 1020 to check for z1 resonance Added asteroid 8089 to check for z2 resonance Running simulation...
100%|██████████| 125663/125663 [02:47<00:00, 749.55it/s]
============================================================ RESULTS: Secular Resonance Checks ============================================================ z1 resonance (g-g6): Asteroid 847: ? TRANSIENT Asteroid 363: ✓ YES - Strong libration detected (in resonance) Asteroid 1020: ? TRANSIENT z2 resonance (2(g-g6)+(s-s6)): Asteroid 8089: ? TRANSIENT
AstDyS Catalog¶
By default, the application uses NASA Horizon database for the initial data. However, there is another option for numbered asteroids — AstDyS catalog. The app will download the latest version, convert it, and start using it in simulation.
This is especially useful if you need to run some statistical experiments, such as finding all asteroids trapped in a specific resonance.
import resonances
sim = resonances.Simulation(name='test_astdys', source='astdys')
sim.create_solar_system()
sim.add_body(463, '4J-2S-1', '463 Lola')
import warnings
warnings.filterwarnings("ignore")
Now, let's run the simulation with AstDyS data. Note that we did not specify the date of the simulation — when using AstDyS, the app will take the date of the AstDyS catalog. It is possible to change this but it is not recommended.
sim.run()
Configuring the integrator¶
By default, the package uses SABA(10,6,4)timestamp. However, sometimes, one might want to perform a quick integration using another integrator, i.e., whfast or a precision integrator like ias15. The default functions does not support such a customization. However, it is easy to create a custom simulation and proceed with it.
import resonances
sim = resonances.Simulation(name='test_whfast', date='2025-01-01 00:00', integrator='whfast', dt=0.1)
sim.create_solar_system()
sim.add_body(463, '4J-2S-1', '463 Lola')
Searching NASA Horizons for 'Sun'... Found: Sun (10) Searching NASA Horizons for 'Mercury'... Found: Mercury Barycenter (199) (chosen from query 'Mercury') Searching NASA Horizons for 'Venus'... Found: Venus Barycenter (299) (chosen from query 'Venus') Searching NASA Horizons for 'Earth'... Found: Earth-Moon Barycenter (3) (chosen from query 'Earth') Searching NASA Horizons for 'Mars'... Found: Mars Barycenter (4) (chosen from query 'Mars') Searching NASA Horizons for 'Jupiter'... Found: Jupiter Barycenter (5) (chosen from query 'Jupiter') Searching NASA Horizons for 'Saturn'... Found: Saturn Barycenter (6) (chosen from query 'Saturn') Searching NASA Horizons for 'Uranus'... Found: Uranus Barycenter (7) (chosen from query 'Uranus') Searching NASA Horizons for 'Neptune'... Found: Neptune Barycenter (8) (chosen from query 'Neptune') Searching NASA Horizons for 'Pluto'... Found: Pluto Barycenter (9) (chosen from query 'Pluto') Searching NASA Horizons for '463;'... Found: 463 Lola (A900 UK)
sim.run()
Accessing the initial data¶
The package stores the initial and the integration data in memory. One might want to access the initial data.
bodies = sim.bodies
print(bodies)
Each body is an instance of the Body class. One might want to access the initial data of a given body.
body = sim.bodies[0]
print(body.initial_data)
print(body.axis)
print(body.ecc)
print(body.inc)
print(body.Omega)
print(body.omega)
print(body.M)
print(body.longitude)
print(body.varpi)
print(body.angles)
print(body.mmrs)
The full list of available methods and properties is available in the Body class.
Find only two-body MMRs near the semimajor axis¶
In Quick Start, we have seen how to find all MMRs (two-body and three-body) for a given value of semi-major axis through the function resonances.find_resonances. However, one might want to find only two-body or three-body resonances. To perform this, there are special classes ThreeBodyMatrix and TwoBodyMatrix that can be used to find only three-body and two-body resonances, respectively.
from resonances import ThreeBodyMatrix, TwoBodyMatrix
mmrs = TwoBodyMatrix.find_resonances(2.35, sigma=0.1, planets=['Mars', 'Jupiter'])
for mmr in mmrs:
print(mmr.to_short())
mmrs = ThreeBodyMatrix.find_resonances(2.35, sigma=0.05, planets=['Mars', 'Jupiter', 'Saturn'])
for mmr in mmrs:
print(mmr.to_short())
1M-2 1M+2 5M-9 6M-11 7J-2 7J+2 10J-3 1M-3J-1 2M-3J-3 2M+4J-5 3M-2J-5 3M+1J-6 3M+4J-7 4M-2J-7 1M+1S-2 2M+1S-4 3M+1S-6 3M+2S-6 2J+3S-1 4J-2S-1 5J-4S-1 6J+1S-2 7J-1S-2
The code above performs the following:
- The first
find_resonancesfinds two-body MMRs resonances within the range $2.35\pm 0.1 = [2.25, 2.45]$. - The second
find_resonancesfinds three-body MMRs resonances within the range $2.35\pm 0.05 = [2.3, 2.4]$.
Find asteroids in a given resonance¶
One might want to find all asteroids in a given resonance. This can be done by the function find_asteroids_in_resonance. It has the following inputs:
mmr(Union[MMR, str]): The mean motion resonance to search for. Can be either an MMR object or a string representation (e.g., "4J-2S-1")sigma(float, default=0.1): Width parameter for resonance search around the resonant axisper_iteration(int, default=500): Number of asteroids to process in each batch/chunk for memory efficiency
The function returns a list of DataFrames containing simulation results for asteroids found in the specified resonance.
This function is time-consuming because it integrates the orbits of all relevant asteroids. One might want to get a list of candidates first based on the closeness of the resonant axis. To do this, we can use astdys package.
import astdys, resonances
mmr = resonances.create_mmr('4J-2S-1')
df_asteroids = astdys.search_by_axis(mmr.resonant_axis, sigma=0.01)
print('Number of objects found: {}'.format(len(df_asteroids)))
df_asteroids.head(5)
Number of objects found: 3619
| a | e | inc | Omega | omega | M | epoch | |
|---|---|---|---|---|---|---|---|
| num | |||||||
| 20 | 2.407925 | 0.143731 | 0.012380 | 3.595020 | 4.493478 | 4.962143 | 60600.0 |
| 25 | 2.400450 | 0.254243 | 0.377134 | 3.736615 | 1.574323 | 5.989823 | 60600.0 |
| 60 | 2.392187 | 0.184544 | 0.062843 | 3.342778 | 4.725915 | 0.377345 | 60600.0 |
| 63 | 2.394722 | 0.128203 | 0.100759 | 5.893842 | 5.161197 | 5.439140 | 60600.0 |
| 192 | 2.402742 | 0.245596 | 0.118624 | 5.988033 | 0.533790 | 2.924088 | 60600.0 |
The function search_by_axis is powerful and can be used to find asteroids by a given value of the semi-major axis for other purposes as well (i.e., to plot some of them).
If one still wants to classify whether any of the asteroids found are really resonant, one can use the function find_asteroids_in_mmr described above.
data = resonances.find_asteroids_in_mmr('6J-3S-2', sigma=0.01)
By default, this method also saves output data to cache/%current_datetime% subdirectory. You can access the results either through the variable data or by reading summary.csv file in the proper directory.