main.py

import json
from math import radians

import numpy as np
from rebound import Particle, Simulation
from reboundx import Extras

data = [
    [5.202887e0, 0.04838624e0, 1.30439695e0, 274.2545707e0, 100.4739091e0, 19.66796068e0, 0.000954792e0, "JUPITER"],
    [9.53667594e0, 0.05386179e0, 2.48599187e0, 338.9364538e0, 113.6624245e0, 317.3553659e0, 0.000285886e0, "SATURN"],
    [19.18916464e0, 0.04725744e0, 0.77263783e0, 96.93735127e0, 74.01692503e0, 142.2838282e0, 4.36624e-05, "URANUS"],
    [30.06992276e0, 0.00859048e0, 1.77004347e0, 273.1805365e0, 131.7842257e0, 259.915208e0, 5.15139e-05, "NEPTUNE"]
]


def main():
    sim = Simulation()
    rebx = Extras(sim)
    rebx.register_param("water", "REBX_TYPE_DOUBLE")
    sim.units = ('yr', 'AU', 'Msun')
    sim.add(m=1.0)  # Sun
    for planet in data:
        part = Particle(a=planet[0], e=planet[1], inc=radians(planet[2]), omega=radians(planet[3]),
                        Omega=radians(planet[4]), M=radians(planet[5]), m=planet[6], simulation=sim)
        sim.add(part)
    max_n = sim.N
    sim.particles[1].params["water"] = 0.3
    print("start")
    tmax = 1e5
    savesteps = 1000
    times = np.linspace(0., tmax, savesteps)
    sim.move_to_com()
    sim.automateSimulationArchive("sa.bin", interval=savesteps)
    for i, t in enumerate(times):
        sim.integrate(t)
        print(f"done: {i / 10}% ({sim.N}")
    with open("sa.meta.json", "w") as f:
        meta = {"tmax": tmax, "savesteps": savesteps, "max_n": max_n}
        json.dump(meta, f)


if __name__ == '__main__':
    main()