Getting too many unbound orbits

[ssah1378]

Hi,

I am defining a system with Sun, Jupiter, 2,000 comets and a passing star and when I plot the histogram of the eccentricities I am getting almost all the comets being unbound and in hyperbolic orbits, while the theories say for the given distance and speed of the passing star, the perturbations should not be strong enough to dump almost all the comets into hyperbolic orbits. Would you be kind to give some suggestions with regards to what could be wrong with the code that I have written?

Here is the code: import rebound import numpy as np import random import matplotlib.pyplot as plt import pandas as pd import scipy as sp import scipy.stats as st from itertools import repeat import time

sec_beg = time.time()

creating values for eccentricity

eccentricities = [] for i in range(0,1000): eccentricities.extend(repeat(i/1000, 2*(1000-i)-1)) #the '-1' works but why?

creating values for the semi-major axis

semimajor = [] for i in range(0,1000): semimajor.extend(repeat((104 + i*10*2), int(91011/((10*4 + i102)2)-0.15)+178))

creating a list for both eccentricity and semimajor, randomly assigning them to eachother

eanda = [] for i in range(0,1000000): eanda.append([semimajor[random.randint(0,999999)], eccentricities[random.randint(0,999999)]])

sim = rebound.Simulation() sim.units = ('yr','AU','Msun') sim.integrator = "whfast" sim.dt=10

sim.add(m=1)

sim.add(m=1e-3, a=5.2, e=0.05)

p = rebound.Particle() p.m = 0.5 p.x = 20000 p.y = 17500000 p.vy = -7 sim.add(p)

comets_init = [] count = 0 while count < 2000: count += 1 sim.add(a=float(eanda[count-1][0]), #N_a = Constant*(1 - a_inner / a) e=float(eanda[count-1][1]), #F_e = 1 - e^2 inc=float(random.uniform(-np.pi/2, np.pi/2)), #isotropy Omega=float(random.uniform(0,np.pi))) #isotropy comets_init.append([sim.particles[count+2].a, sim.particles[count+2].e, sim.particles[count+2].inc, sim.particles[count+2].Omega])

sim.N_active = 3 sim.testparticle_type = 0

sim.integrate(5000000)

sec_end = time.time() print("Runtime =", sec_end - sec_beg)

comets_fin = [] for orbits in sim.calculate_orbits(): comets_fin.append([orbits.a, orbits.e, orbits.inc, orbits.Omega])

pd.DataFrame(np.array(comets_init), columns=['a', 'e', 'i', 'Omega']).to_csv('comets_init_0.csv')

pd.DataFrame(np.array(comets_fin), columns=['a', 'e', 'i', 'Omega']).to_csv('comets_fin_0.csv')

Best, Shasha

[ssah1378]

I case it might help:

as you can see then comets are getting unbound way before the passing particle reaches its perihelion, as if the test particles are not bound to the central mass at all!

[ssah1378]

I needed to calculate everything in heliocentric coordinates! my bad :D

[hannorein]

Glad you've figured it out!