Bug? Rebound doesn't plot the same orbit if I add speed to my planetary system in 2 differents method

[Francyrad]

Dear developers I wanted to study the trajectory of an asteroidal system following a SPH collision. So, i did the simulation and extract the velocities and masses to plot the system following the collision around the Sun. for example, at a distance of 1 AU.

So, i wanted to add a speed of vx of -29785 m/s in my system supponing that my target was in an orbit of 1 AU with e=0 befor collision

It is not important if I add this speed upstream, so when I extract my data, then plotting them on Rebound, "case A" (simulation velocities + orbital velocities --> REBOUND = orbit) or if I decide to load the raw speeds to REBOUND and then add the orbital speed with REBOUND "case B" (simulation velocities --> REBOUND + REBOUND_orbital_velocity = orbit). The result should always be the same orbit.

And yet I still get a difference that I can’t figure out. I can’t tell if there is any problem in REBOUND summing up the speeds

I also tried to run my simulation by implementing orbital speed added to the impact speed (so my system is comoving with a speed of 29785 m/s). So I’m sure solution A is the correct one. I just can’t find an error of why solution B is wrong if I add the same speed with REBOUND.

Here is part of the script:

import h5py
import os
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
import swiftsimio as sw
import unyt
import woma
from woma.misc import utils
import rebound

CASE A

def load_data(snapshot):
    """
    Load data from the snapshot and calculate quantities of interest.
    """
    # Load data
    data = sw.load(snapshot)
    box_mid = 0.5 * data.metadata.boxsize[0].to(unyt.m)
    pos = data.gas.coordinates.to(unyt.m) - box_mid
    pos = np.array(pos)
    vel = np.array(data.gas.velocities)
    data.gas.velocities.convert_to_mks()
    data.gas.masses.convert_to_mks()
    m = np.array(data.gas.masses)
    data.gas.masses.convert_to_mks()
    vel[:, 0] += -29785/(1e6)      #I ADD THE SPEED HERE

    print("Posizioni iniziali:", pos[:10])  # Stampa le prime 10 posizioni per avere un'idea
    print("Velocità iniziali:", vel[:10])  # Stampa le prime 10 velocità per avere un'idea

    pos_centerM = np.sum(pos * m[:, np.newaxis], axis=0) / np.sum(m)
    #vel_centerM = np.sum(vel * m[:, np.newaxis], axis=0) / np.sum(m)
    print("Centro di massa (posizione):", pos_centerM)

    pos -= pos_centerM

    print("Posizioni corrette:", pos[:10])
    #vel[:, 0] += -29785/(1e6)
    #vel -= vel_centerM
    return pos, vel, m

snapshot_path = find_latest_hdf5_file('output')
if not snapshot_path:
    raise ValueError("No HDF5 file found in 'output' folder.")

# Load data from the latest snapshot
pos, vel, m = load_data(snapshot_path)

print("Massa totale:", np.sum(m))

# Get velocities within 7 Earth radii
within_7_radii = np.sqrt(np.sum(pos**2, axis=1)) <= (6371000 * 7)
print("Velocità entro 7 raggi terrestri:", vel[within_7_radii])
vx_cm_s = np.sum(vel[within_7_radii, 0]) * 100  # convert m/s to cm/s
vy_cm_s = np.sum(vel[within_7_radii, 1]) * 100  # convert m/s to cm/s
vz_cm_s = np.sum(vel[within_7_radii, 2]) * 100  # convert m/s to cm/s

# Convert velocities to AU/yr for rebound
vx_AU_yr_simulation = vx_cm_s * (3.154e7 / 1.496e13)
vy_AU_yr_simulation = vy_cm_s * (3.154e7 / 1.496e13)
vz_AU_yr_simulation = vz_cm_s * (3.154e7 / 1.496e13)
print(f"Velocità convertite in AU/yr: vx = {vx_AU_yr_simulation}, vy = {vy_AU_yr_simulation}, vz = {vz_AU_yr_simulation}")

CASE B

def load_data(snapshot):
    """
    Load data from the snapshot and calculate quantities of interest.
    """
    # Load data
    data = sw.load(snapshot)
    box_mid = 0.5 * data.metadata.boxsize[0].to(unyt.m)
    pos = data.gas.coordinates.to(unyt.m) - box_mid
    pos = np.array(pos)
    vel = np.array(data.gas.velocities)
    data.gas.velocities.convert_to_mks()
    data.gas.masses.convert_to_mks()
    m = np.array(data.gas.masses)
    data.gas.masses.convert_to_mks()

    print("Posizioni iniziali:", pos[:10])  # Stampa le prime 10 posizioni per avere un'idea
    print("Velocità iniziali:", vel[:10])  # Stampa le prime 10 velocità per avere un'idea

    pos_centerM = np.sum(pos * m[:, np.newaxis], axis=0) / np.sum(m)
    #vel_centerM = np.sum(vel * m[:, np.newaxis], axis=0) / np.sum(m)
    print("Centro di massa (posizione):", pos_centerM)

    pos -= pos_centerM

    print("Posizioni corrette:", pos[:10])
    #vel[:, 0] += -29785/(1e6)
    #vel -= vel_centerM
    return pos, vel, m

snapshot_path = find_latest_hdf5_file('output')
if not snapshot_path:
    raise ValueError("No HDF5 file found in 'output' folder.")

# Load data from the latest snapshot
pos, vel, m = load_data(snapshot_path)

print("Massa totale:", np.sum(m))

# Get velocities within 7 Earth radii
within_7_radii = np.sqrt(np.sum(pos**2, axis=1)) <= (6371000 * 7)
print("Velocità entro 7 raggi terrestri:", vel[within_7_radii])
vx_cm_s = np.sum(vel[within_7_radii, 0]) * 100  # convert m/s to cm/s
vy_cm_s = np.sum(vel[within_7_radii, 1]) * 100  # convert m/s to cm/s
vz_cm_s = np.sum(vel[within_7_radii, 2]) * 100  # convert m/s to cm/s

# Convert velocities to AU/yr for rebound
vx_AU_yr_simulation = ((-29785 * 100) * (3.154e7 / 1.496e13)) + (vx_cm_s * (3.154e7 / 1.496e13)) #I ADD THE SPEED HERE
vy_AU_yr_simulation = vy_cm_s * (3.154e7 / 1.496e13)
vz_AU_yr_simulation = vz_cm_s * (3.154e7 / 1.496e13)
print(f"Velocità convertite in AU/yr: vx = {vx_AU_yr_simulation}, vy = {vy_AU_yr_simulation}, vz = {vz_AU_yr_simulation}")

# Corrected conversion for the provided vx
vx_AU_yr = 0
vy_AU_yr = +29785 * (3.154e7 / 1.496e11)
vz_AU_yr = 0

# Convert mass to Msun for rebound
total_mass_msun = np.sum(m[within_7_radii]) / 1.9891e30

# Set up the simulation
sim = rebound.Simulation()
sim.units = ('yr', 'AU', 'Msun')
sim.add(m=1.0)  # Sun

# Add the planet with the velocity and mass from the simulation
sim.add(m=total_mass_msun, x=0, y=1, z=0, vx=vx_AU_yr_simulation, vy=vy_AU_yr_simulation, vz=vz_AU_yr_simulation)

# Add the additional planet with given vx
sim.add(m=1e-99, x=1, y=0, z=0, vx=vx_AU_yr, vy=vy_AU_yr, vz=vz_AU_yr, hash="Terra")

planet = sim.particles[1]
# After adding the bodies to the simulation, extract the planet's orbit.
planet_orbit = planet.calculate_orbit(primary=sim.particles[0])

# Plot the orbit using rebound
op = rebound.OrbitPlot(sim, unitlabel="[AU]", xlim=(-2, 2), ylim=(-2, 2), color=True)
op.primary.set_edgecolor("yellow")

# Plot the orbit of particle 1 with green color on the same plot
op1 = rebound.OrbitPlot(sim, particles=[1], ax=op.ax, fig=op.fig, lw=2, color="green")

# Calculate and display the eccentricity
eccentricity = planet_orbit.e
print("Eccentricità:", eccentricity)
op.ax.text(1.5, 1.8, f"e={eccentricity:.2f}", color="black", ha="center", va="center", bbox=dict(facecolor='white', edgecolor='none'))

# Saving the plot with the desired name
basename = os.path.basename(snapshot_path)
output_name = os.path.splitext(basename)[0] + "_orbit_post_impact.png"
op.fig.savefig(output_name)

plt.tight_layout()
plt.show()

CASE A result (correct)

CASE B result (wrong)

I'm I doing something wrong?

Best regards Francesco

[hannorein]

Hi Francesco,

I don't fully understand the issue you're having. The code is quite long and difficult to read. If you can, please provide a shorter example. But here are a few tips that might help you:

1) Instead of plotting orbits, you might just want to print out some cartesian coordinates of the particles to find where the problem is coming from. 2) If you have particles with a finite mass and you're working with orbital parameters, then the order in which you add particles matters because REBOUND assume Jacobi coordinates by default. I see that you use the primary argument in your code to calculate heliocentric orbital elements, so I assume you are aware of the subtleties. 3) Similarly, maybe there is an issue with calculating the orbital speed because more than one particle is massive? 4) There are a lot of magic numbers in your code. It's impossible for me to check where they are coming from. And units in general can often lead to some confusion. So double checking those might help as well.

I can look into it more, but I'd need to have a shorter example that doesn't rely on other files and I can run myself. Hanno

[Francyrad]

It's not a rebound issue, the problem came from Swiftsim and Swiftsimio. There is a conversion error apparently... So i close for the moment

thank you for the help and availability Francesco

[hannorein]

Hi Francesco, I'm glad you found where the problem was coming from! Feel free to reopen the issue if something else comes up.

Hanno