import rebound
import reboundx
from reboundx import constants

sim = rebound.Simulation()
sim.add(m=1.)
sim.add(m=1e-3, a=1., e=0.2)
sim.add(m=1e-3, a=1.9)
sim.move_to_com()
sim.dt = sim.particles[1].P*0.05  # timestep is 5% of orbital period
sim.integrator = "whfast"
rebx = reboundx.Extras(sim)
gr = rebx.load_force("gr_potential")
rebx.add_force(gr)
gr.params["c"] = constants.C
rebx.save("rebxarchive.bin")

"""
Explanation: SimulationArchives with REBOUNDx
We can reproduce REBOUNDx simulations with extra effects bit by bit with the SimulationArchive (see Rein & Tamayo 2017 for details), and make sure to read SimulationArchive.ipynb and SimulationArchiveRestart.ipynb first.
However, this will only work under some conditions. In particular, for bit-by-bit reproducibility, currently one requires that:
- All REBOUNDx effects in the simulation are machine independent (can't have, e.g., trig functions, pow or exp in implementation)
- The effect and particle parameters cannot change throughout the simulation
- Effects must remain on for the entire integration
To use the SimulationArchive with REBOUNDx, we need to save a REBOUNDx binary SavingAndLoadingSimulations.ipynb. Since we can't change effects or particle parameters, it doesn't matter at what point we save this binary (we will just need it to load the SimulationArchive). Here we do a WHFAST integration with the symplectic gr_potential implementation for general relativity corrections. We set up the simulation like we usually would:
End of explanation
"""


sim.automateSimulationArchive("archive.bin", interval=1e3, deletefile=True)
sim.integrate(1.e6)

"""
Explanation: We now set up the SimulationArchive and integrate like we normally would (SimulationArchive.ipynb):
End of explanation
"""


sa = reboundx.SimulationArchive("archive.bin", rebxfilename = "rebxarchive.bin")

"""
Explanation: Once we're ready to inspect our simulation, we use the reboundx.SimulationArchive wrapper that additionally takes a REBOUNDx binary:
End of explanation
"""


sim2, rebx = sa[500]
sim2.t

"""
Explanation: We now create a different simulation from a snapshot in the SimulationArchive halfway through:
End of explanation
"""


sim2.integrate(1.e6)

"""
Explanation: We now integrate our loaded simulation to the same time as above (1.e6):
End of explanation
"""


sim.status()

sim2.status()

"""
Explanation: and see that we obtain exactly the same particle positions in the original and reloaded simulations:
End of explanation
"""


