<jupyter_start><jupyter_text>Integrating eccentric Comets with HERMESIn this example, we study highly eccentric comets which interact with a Neptune mass planet.HERMES, a hybrid integration scheme which combines the WHFAST and IAS15 algorithms, is well suited to study this problem, since the algorithm will automatically switch to IAS15 during the perihelion passage of a comet or a close encounter with the planet but otherwise will use a large timestep with the symplectic WHFast integrator.<jupyter_code>import rebound
import numpy as np<jupyter_output><empty_output><jupyter_text>First let's choose the basic properties required for the HERMES integrator to run correctly. In particular, we are:  * Adding comets as *semi-active* bodies, which means they can influence/be influenced by other active bodies, but are invisible to each other. This is done by setting testparticle_type to 1. Setting testparticle_type to 0 would indicate that we are adding comets as *test* bodies.* Colliding bodies merge perfectly conserving momentum and mass.* Particles that leave our pre-defined box get removed from the simulation.* We track the energy lost due to ejections or collisions.<jupyter_code>sim = rebound.Simulation()
np.random.seed(42)

#integrator options
sim.integrator = "hermes"
sim.ri_hermes.solar_switch_factor = 50.    #boundary for close encounters with sun. Units of physical radius.
sim.ri_hermes.hill_switch_factor = 4.      #Setting HSF "floor", adaptive HSF routine is active by default.
sim.dt = 1
sim.testparticle_type = 1

#collision and boundary options
sim.collision = "direct"
sim.collision_resolve = "merge"
sim.collision_resolve_keep_sorted = 1
sim.boundary = "open"
boxsize = 200.
sim.configure_box(boxsize)
sim.track_energy_offset = 1

#simulation time
tmax = 1e4<jupyter_output><empty_output><jupyter_text>Now that the preliminary setup is complete, it's time to add some particles to the system! When using the HERMES integrator it's always important to add active bodies first and semi-active bodies later. The `sim.N_active` variable distinguishes massive bodies from semi-active/test bodies.<jupyter_code>#massive bodies
sim.add(m=1., r=0.005)                     #sun
a_neptune = 30.05
sim.add(m=5e-5,r=2e-4,a=a_neptune,e=0.01)  #neptune

sim.N_active = sim.N<jupyter_output><empty_output><jupyter_text>Now, let's create some comets! For this simple example we are assuming that all comets have the same mass and radius.<jupyter_code>#semi-active bodies
n_comets = 100
a = np.random.random(n_comets)*10 + a_neptune
e = np.random.random(n_comets)*0.009 + 0.99
inc = np.random.random(n_comets)*np.pi/2.
m=1e-10
r=1e-7

for i in xrange(0,n_comets):
    rand = np.random.random()*2*np.pi
    sim.add(m=m, r=r, a=a[i], e=e[i], inc=inc[i], Omega=0, omega=rand, f=rand)<jupyter_output><empty_output><jupyter_text>We need to move to the COM frame to avoid drifting out of our simulation box. Also it is always good practice to monitor the change in energy over the course of a simulation, which requires us to calculate it before and after the simulation.<jupyter_code>sim.move_to_com()
E0 = sim.calculate_energy()<jupyter_output><empty_output><jupyter_text>We can visualize our setup using `rebound.OrbitPlot`<jupyter_code>%matplotlib inline
fig = rebound.OrbitPlot(sim, trails=True)<jupyter_output><empty_output><jupyter_text>Finally, let's simulate our system for and check that our final relative energy error is small. The energy error is a key measure of whether the integration was performed accurately or not.<jupyter_code>sim.integrate(tmax)
dE = abs((sim.calculate_energy() - E0)/E0)
print dE<jupyter_output>6.01928861878e-07
{'fraction_non_alphanumeric': 0.08034982235583493, 'fraction_numerical': 0.02022410494670675, 'mean_word_length': 5.50088809946714, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 28, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '238805', 'license': 'MIT', 'n_tokens_mistral': 1215, 'n_tokens_neox': 1136, 'n_words': 410, 'path': 'rebound/ipython_examples/EccentricComets.ipynb', 'repo_name': 'rodluger/ttv-devil', 'size': '3659'}