hannorein / rebound

💫 An open-source multi-purpose N-body code.
https://rebound.readthedocs.io/
GNU General Public License v3.0
824 stars 218 forks source link

Close encounter distance of a fly-by differs with planet system inclination #622

Closed Anit5577 closed 1 year ago

Anit5577 commented 1 year ago

I am currently running a very simple simulation where a star with a planetary system of 5 planets has a single close encounter with another star. The fly-by should be governed by the two stars in the simulation. However, the closest encounter distance between the star (with the planetary system) and the fly-by star appears to depend on the inclination of the planet system (very different orbits).

import rebound
import random
import numpy as np

sim = rebound.Simulation()
sim.G = 6.6743e-11
sim.units = ('days', 'AU', 'Msun')
sim.integrator = "ias15"

# using different values for the inclination of the planet system leads to massive differences in the orbit)
inc_val = np.pi 

sim.add(m=1.2, x=0.0, y=0.0, z=0.0, vx=0.0, vy=0.0 ,vz=0.0, r = 0.00511551, hash='Star')
sim.add(primary=sim.particles[0],m=5*3e-6, e=0.0, inc = inc_val, a=0.1,r = 7e-5, hash='P1') # Super-Earth
sim.add(primary=sim.particles[0],m=5*3e-6, e=0.0, inc = inc_val,a=0.1232, r = 7e-5,hash='P2') # Super-Earth
sim.add(primary=sim.particles[0],m=5*3e-6, e=0.0, inc = inc_val,a=0.1518, r = 7e-5,hash='P3') # Super-Earth
sim.add(primary=sim.particles[0],m=5*3e-6, e=0.0, inc = inc_val, a=0.1870, r = 7e-5,hash='P4') # Super-Earth
sim.add(primary=sim.particles[0],m=5e-3,  e=0.0, inc = inc_val, a=10, r = 7e-4,hash='G1') # Cold Jupiter

I am then adding a perturber with its cartesian coordinates. Despite the position and velocity of the perturber being given, its actual orbit past the planetary system depends greatly on the inclination of the planetary system (inc_val above).

# Perturbing star - Fly-by
sim.add(m=0.5, x=7000, y=4554, z=1614,vx=-0.002, vy=-0.0012, vz=-0.0004,r = 0.002, hash='Pert') 

image with inc_val = 1

image with inc_val = 2

When I add the perturber with orbital elements instead, the fly-by remains exactly the same regardless of inclination.

sim.add(m=0.5,a=-98.79795974663652, e=1.9952267131093615, inc=0.8549416902934751,\
      Omega=0.4078757723646595, omega=2.3300573913127174, f=4.2073496026327515,\
    r = 0.002, hash='Pert') 

sim.move_to_com()

I am then running the simulation for 10,000 years, with the fly-by happening right towards the end.

sim_time_init = 1.0e+4*365.25
sim.automateSimulationArchive('integration_test.bin', interval = (sim_time_init/100))

sim.integrate(sim_time_init)

There appears to be no physical reason for this behaviour and it does not appear to happen with all set-ups. Is the mixing of orbital and cartesian elements to blame (despite a mixed example being present on the website "Convenience functions")?

hannorein commented 1 year ago

Is it possible that this occurs because the centre of mass of the star + planets depends on the inclination, so the planetary system is drifting in slightly different directions depending on inc_val. Even if they drift slowly, because it takes a while for the perturber to get to pericenter, this might be a big effect. Try calling sim.move_to_com() after adding all the planets (before adding the perturber).

Anit5577 commented 1 year ago

Hi Hanno, changing the positioning of sim.move_to_com() to earlier seems to have the desired effect. It moves the planetary system's COM to the origin so that the encounter always happens in the same way. The original issue seemed to have been caused by the difference in the drift of the COM of the planetary system depending on its inclination.

hannorein commented 1 year ago

I'm glad the suggestion worked. Closing this issue for now but feel free to reopen it later if you have other questions!