How to plot multiple "alternative" orbit evolutions?

[Francyrad]

Dear @hannorein

Dear @hannorein,

I hope this message finds you well.

I've developed a script that visualizes the orbit of my system around the Sun post-collision simulation (illustrated in dark green). Following this, I made some modifications to explore various "alternative scenarios." This included altering the collision direction and adjusting the preferred orbit distance around the sun, starting from AU = Y = 1.

You can view the current representation here: My goal is to simultaneously display the orbit evolution of each alternative scenario on a single graph. Importantly, I don't want the alternative bodies in the simulation to influence each other. The desired output is something similar to what you showcased in the TESLA example:

However, I've encountered some challenges:

When adding multiple alternative scenarios, attempting to display all of them results in a blank (white) screen, while i'd love to see the evolution of each body influenced only from other bodies of the solar system (and not from their conterparts). Is there a method within REBOUND to address this?

For reasons I can't discern, all of my alternative orbits appear in varying colors instead of the desired green. Additionally, adding them seems to alter the coloration of other planets.

Attached, please find a simplified code version that doesn't rely on the original simulation for easier reference.

I deeply appreciate your assistance with these matters and look forward to any insights you can provide.

Thank you for your support and best regards Francesco

import os
import rebound
import numpy as np
import matplotlib.pyplot as plt

def manipulate_velocity(vx_sim, vy_sim, direction):
    """
    Manipulate the velocity based on the desired impact direction.
    :param vx_sim: velocity in x from simulation (before modification)
    :param vy_sim: velocity in y from simulation (before modification)
    :param direction: desired direction for the impact ("+x", "+y", "-x", "-y", "+x+y", "+x-y", "-x+y", "-x-y")
    :return: manipulated velocities (vx, vy)
    """

    if direction == "+x":
        return vx_sim, vy_sim
    elif direction == "+y":
        return -vy_sim, vx_sim
    elif direction == "-x":
        return -vx_sim, -vy_sim
    elif direction == "-y":
        return vy_sim, -vx_sim
    elif direction == "+x+y":
        # Assuming a 45 degree impact for simplicity.
        return np.sqrt((vx_sim**2 + vy_sim**2) / 2), np.sqrt((vx_sim**2 + vy_sim**2) / 2)
    elif direction == "+x-y":
        return np.sqrt((vx_sim**2 + vy_sim**2) / 2), -np.sqrt((vx_sim**2 + vy_sim**2) / 2)
    elif direction == "-x+y":
        return -np.sqrt((vx_sim**2 + vy_sim**2) / 2), np.sqrt((vx_sim**2 + vy_sim**2) / 2)
    elif direction == "-x-y":
        return -np.sqrt((vx_sim**2 + vy_sim**2) / 2), -np.sqrt((vx_sim**2 + vy_sim**2) / 2)
    else:
        raise ValueError("Unsupported direction!")

# Imposta le velocità iniziali manualmente
vx_cm_s_original = -560943.494757  # in cm/s
vy_cm_s_original = 183172.861813   # in cm/s

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

# Modifica la velocità in base alla direzione desiderata e converte in AU/yr per rebound
directions = ["+x", "+y", "-x", "-y", "+x+y", "+x-y", "-x+y", "-x-y"]

for direction in directions:
    # Modifica la velocità per la direzione corrente

    vx_cm_s, vy_cm_s = manipulate_velocity(vx_cm_s_original, vy_cm_s_original, direction)
    print(f"Velocità modificata per direzione {direction}: vx = {vx_cm_s}, vy = {vy_cm_s}")
    vx_cm_s += -29785 * 100  # convert m/s to cm/s
    vx_AU_yr_simulation = vx_cm_s * (3.154e7 / 1.496e13)
    vy_AU_yr_simulation = vy_cm_s * (3.154e7 / 1.496e13)
    print(f"Velocità in AU/yr per direzione {direction}: vx = {vx_AU_yr_simulation}, vy = {vy_AU_yr_simulation}")

    # Aggiungi il pianeta con la velocità modificata alla simulazione
    sim.add(m=3.00273e-6, x=0, y=1, z=0, vx=vx_AU_yr_simulation, vy=vy_AU_yr_simulation, vz=0)

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

# 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")

# Mercurio
sim.add(a=0.3871, e=0.2056, inc=0.020, m=1.65e-7, hash="Mercury")

# Venere
sim.add(a=0.723, e=0.007, inc=0.059, m=4.8675e-6, hash="Venus") # Venus

# Marte
sim.add(a=1.524, e=0.093, inc=0.032, m=3.214e-7, hash="Mars") # Mars

# Giove
sim.add(a=5.2044, e=0.0489, inc=0.0228, m=0.0009543, hash="Jupiter")  # Giove

# Saturno
sim.add(a=9.5826, e=0.0565, inc=0.0433, m=0.0002857, hash="Saturn")  # Saturno

# Urano
sim.add(a=19.2184, e=0.0464, inc=0.0135, m=0.00004365, hash="Uranus")  # Urano

# Nettuno
sim.add(a=30.1104, e=0.0095, inc=0.0309, m=0.0000515, hash="Neptune")  # Nettuno

# Plutone (nota: Plutone ha un'orbita molto inclinata e eccentrica)
sim.add(a=39.482, e=0.2488, inc=17.14175, m=6.55e-9, hash="Pluto")  # Plutone

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")

op2 = rebound.OrbitPlot(sim, particles=[2], ax=op.ax, fig=op.fig, lw=2, color="gray")        # Mercury
op3 = rebound.OrbitPlot(sim, particles=[3], ax=op.ax, fig=op.fig, lw=2, color="beige")     # Venus

op4 = rebound.OrbitPlot(sim, particles=[4], ax=op.ax, fig=op.fig, lw=2, color="lightblue")       # Earth
op5 = rebound.OrbitPlot(sim, particles=[5], ax=op.ax, fig=op.fig, lw=2, color="red")        # Mars

op6 = rebound.OrbitPlot(sim, particles=[6], ax=op.ax, fig=op.fig, lw=2, color="burlywood")     # Jupiter
op7 = rebound.OrbitPlot(sim, particles=[7], ax=op.ax, fig=op.fig, lw=2, color="blanchedalmond")  # Saturn
op8 = rebound.OrbitPlot(sim, particles=[8], ax=op.ax, fig=op.fig, lw=2, color="dodgerblue") # Uranus
op9 = rebound.OrbitPlot(sim, particles=[9], ax=op.ax, fig=op.fig, lw=2, color="blue")       # Neptune
op10 = rebound.OrbitPlot(sim, particles=[10], ax=op.ax, fig=op.fig, lw=2, color="brown")  # Pluto
# Explicitly setting colors based on the planets

# 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

output_name = "orbits_test.png"
op.fig.savefig(output_name)

plt.tight_layout()
plt.show()

# Salvare la simulazione iniziale in un archivio binario
sim.save("initial_simulation.bin")

# Impostare l'integratore a "ias15"
sim.integrator = "ias15"

directions = ["+x", "+y", "-x", "-y", "+x+y", "+x-y", "-x+y", "-x-y"]

# Liste per conservare i dati delle traiettorie
a_values = {direction: [] for direction in directions}
e_values = {direction: [] for direction in directions}

# Tempo di integrazione (ad esempio, 100 anni - puoi modificarlo come desideri)
integration_time = 100000 * 2.0 * np.pi

for direction in directions:
    # Caricare la simulazione iniziale dall'archivio binario
    sim = rebound.Simulation.from_archive("initial_simulation.bin")

    # Manipolare la velocità della particella 1
    planet = sim.particles[1]
    planet.vx, planet.vy = manipulate_velocity(planet.vx, planet.vy, direction)

    # Integrazione
    N_points = 1000
    times = np.linspace(0, integration_time, N_points)
    a_temp = []
    e_temp = []
    for t in times:
        sim.integrate(t)
        orbit = planet.calculate_orbit(primary=sim.particles[0])  # Specifico la particella primaria
        a_temp.append(orbit.a)
        e_temp.append(orbit.e)

    a_values[direction] = a_temp
    e_values[direction] = e_temp

        # Stampare i primi valori per la direzione corrente come debug
    print(f"Per la direzione {direction}:")
    print(f"Primi 10 valori di a: {a_temp[:10]}")
    print(f"Primi 10 valori di e: {e_temp[:10]}")
    print("------")

# Plot delle traiettorie
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 10))
for direction in directions:
    label = f"Direzione {direction}"
    ax1.plot(times, a_values[direction], label=label)
    ax2.plot(times, e_values[direction], label=label)

ax1.set_title("Semi-major Axis Evolution")
ax1.set_xlabel("Time [orbits]")
ax1.set_ylabel("Semi-major Axis [AU]")

ax2.set_title("Eccentricity Evolution")
ax2.set_xlabel("Time [orbits]")
ax2.set_ylabel("Eccentricity")

ax1.legend()
ax2.legend()

plt.tight_layout()
plt.show()

[hannorein]

I'm not sure if your question is more about the simulation itself or the visualization but I think the solution is to run separate simulations, one for each of your alternatives. You can add multiple test particles with zero mass that don't affect each other. But your particle has mass and will therefore influence the planets as well as other the alternative particles.

[Francyrad]

Thank you for the suggestion! I solved saving each simulations in a txt and then loading them all in once, i didn't think about that. Still thank you for the help. When my computer will be free from another sim i'll test the new mpi fix for the moon origin sim.