search

hannorein / rebound

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

Problems with reproducibility / inconsistent Keplerian orbital elements from #567

Closed SuperdoerTrav closed 3 years ago

SuperdoerTrav commented 3 years ago

tldr; I find that I receive different Keplerian orbital elements when calculating them in post versus printing them as a REBOUND simulation is running.

I have run a simulation of the full solar system and some test particle asteroids, outputting the X Y Z and VX VY VZ of all particles.

The I've initialized my simulations as such:

sim = rebound.Simulation();
sim.integrator = "IAS15"#"whfast"#
sim.testparticle_type = 0 #test particles do not affect gravity
#Initialize planets

#values obtained from NASA ephemeris for date = "JD2458475.30035" 
sim.add(m=1.0, x=-0.006159568320344337, y=0.006384245192586575, z=9.048779754365662e-05, vx=-0.00042050521613923867, vy=-0.00029933766519877264, vz=1.265154869865768e-05, hash=0)
sim.add(m=1.6601141530543488e-07, x=0.31273969799010204, y=0.10980298144789544, z=-0.020711686152730602, vx=-0.8221345868837158, vy=1.628278370608211, vz=0.20846987935406125, hash=1)
sim.add(m=2.4478382877847715e-06, x=-0.3535614442371774, y=0.6347005998008065, z=0.02876008064261855, vx=-1.0334025942935463, vy=-0.5758197743923686, vz=0.051724162934945674, hash=2)
sim.add(m=3.040432648022642e-06, x=0.8565315189278313, y=0.5019534070620206, z=6.48029060350566e-05, vx=-0.5148897674306456, vy=0.863166617746706, vz=-2.6293779909606207e-05, hash=3)
sim.add(m=3.2271560375549977e-07, x=1.2715562126894933, y=0.6373036415243706, z=-0.018032088916419903, vx=-0.3294900645141585, vy=0.7985528763425188, vz=0.024826263746752363, hash=4)
sim.add(m=0.0009547919152112404, x=2.6023564581533245, y=-4.39895228385286, z=-0.03997367088596116, vx=0.37214740802398044, vy=0.24416228247932154, vz=-0.009338160992015496, hash=5)
sim.add(m=0.0002858856727222417, x=5.175812736706754, y=-8.543298821551435, z=-0.057508280769607116, vx=0.25931529161001204, vy=0.1672141496647019, vz=-0.013233697233746537, hash=6)
sim.add(m=4.36624373583127e-05, x=15.51729607179721, y=12.265936110250092, z=-0.15547249594828197, vx=-0.14346578948245664, vy=0.1687160663166264, vz=0.0024852338660758637, hash=7)
sim.add(m=5.151383772628674e-05, x=29.41620671774543, y=-5.442979910051512, z=-0.5658378412865516, vx=0.03199226168880429, vy=0.1805233745901719, vz=-0.004454798282492593, hash=8)

sim.N_active = 9
numplanets = 9
sim.move_to_com()

#Add Asteroids
asters_theta = sim.particles[3].theta + np.random.uniform(0.0001, np.pi-.0001,size=num_asteroids_initial)              
asters_a = np.random.uniform(.99, 1.01,size=num_asteroids_initial)
asters_e = np.abs(np.random.uniform(0, .1,size=num_asteroids_initial))
asters_inc = sim.particles[3].inc + np.random.uniform(0, 2/9*np.pi,size=num_asteroids_initial)                           
asters_omega = sim.particles[3].omega + np.random.uniform(0, 2*np.pi,size=num_asteroids_initial)
asters_Omega = sim.particles[3].Omega + np.random.uniform(0, 2*np.pi,size=num_asteroids_initial)

for i in range(num_asteroids_initial):
    sim.add(m=0,a=asters_a[i],e=asters_e[i],inc=asters_inc[i],omega=asters_omega[i],
        Omega=asters_Omega[i],theta=asters_theta[i],hash=i+sim.N_active)

After running the simulation I wanted to calculate the Keplerian orbital elements, which I then seem to be slightly incorrect. I have both attempted to calculate the Keplerian orbital elements with my own functions and through REBOUND by re-initializing my sims in this manner:

#Initialize planets
sun = np.array(list(loaded_data['0'].values()))[-1]; current_hash.append(int(0))
mercury = np.array(list(loaded_data['1'].values()))[-1]; current_hash.append(int(1))
venus = np.array(list(loaded_data['2'].values()))[-1]; current_hash.append(int(2))
earth = np.array(list(loaded_data['3'].values()))[-1]; current_hash.append(int(3))
mars = np.array(list(loaded_data['4'].values()))[-1]; current_hash.append(int(4))
jupiter = np.array(list(loaded_data['5'].values()))[-1]; current_hash.append(int(5))
saturn = np.array(list(loaded_data['6'].values()))[-1]; current_hash.append(int(6))
uranus = np.array(list(loaded_data['7'].values()))[-1]; current_hash.append(int(7))
neptune = np.array(list(loaded_data['8'].values()))[-1]; current_hash.append(int(8))

sim = rebound.Simulation();
sim.integrator = "IAS15"#"whfast"#
sim.testparticle_type = 0 #test particles do not affect gravity

#Initialize planets
#values obtained from NASA ephemeris for date = "JD2458475.30035" 
sim.add(m=1.0, x=sun[0], y=sun[1], z=sun[2], vx=sun[3], vy=sun[4], vz=sun[5], hash=0)
sim.add(m=1.6601141530543488e-07, x=mercury[0], y=mercury[1], z=mercury[2], vx=mercury[3], vy=mercury[4], vz=mercury[5], hash=1)
sim.add(m=2.4478382877847715e-06, x=venus[0], y=venus[1], z=venus[2], vx=venus[3], vy=venus[4], vz=venus[5], hash=2)
sim.add(m=3.040432648022642e-06, x=earth[0], y=earth[1], z=earth[2], vx=earth[3], vy=earth[4], vz=earth[5], hash=3)
sim.add(m=3.2271560375549977e-07, x=mars[0], y=mars[1], z=mars[2], vx=mars[3], vy=mars[4], vz=mars[5], hash=4)
sim.add(m=0.0009547919152112404, x=jupiter[0], y=jupiter[1], z=jupiter[2], vx=jupiter[3], vy=jupiter[4], vz=jupiter[5], hash=5)
sim.add(m=0.0002858856727222417, x=saturn[0], y=saturn[1], z=saturn[2], vx=saturn[3], vy=saturn[4], vz=saturn[5], hash=6)
sim.add(m=4.36624373583127e-05, x=uranus[0], y=uranus[1], z=uranus[2], vx=uranus[3], vy=uranus[4], vz=uranus[5], hash=7)
sim.add(m=5.151383772628674e-05, x=neptune[0], y=neptune[1], z=neptune[2], vx=neptune[3], vy=neptune[4], vz=neptune[5], hash=8)

sim.N_active = 9
numplanets = 9
sim.move_to_com()

#Add Asteroids
for i in range(astershape[0]):
    ihash = current_hash[i+numplanets]
    sim.add(m=0, x=stable[i][0], y=stable[i][1], z=stable[i][2], \
            vx=stable[i][3], vy=stable[i][4], vz=stable[i][5], hash=ihash)

For example I find the semi-major axis of the Earth deviates from .97 au to 1.03 au, though when printing "sim.particles[hash = 'Earth').a" during an integration the semi-major axis I receive does not deviate from 1 au by more than 0.0001 au. I find this behavior to be occurring in all of my test Earth Trojan asteroids as well. Picture below shows the result:

trajectory-earth-like-1myr-a

If I change the coordinates to heliocentric it seems to make little difference and does not correct the issue.

Calculating the orbital elements manually from REBOUND's XYZ VXVYVZ results in slightly different elements but still shows the same range of oscillation in Earth's semi-major axis.

Calculating orbital elements from xyz vxvvyvz:

    mu = 1.0
    r = np.array([x,y,z])
    v = np.array([vx,vy,vz])

    rmag = np.sqrt(r.dot(r))
    vmag = np.sqrt(v.dot(v))

    h = np.cross(r,v)
    hmag = np.sqrt(h.dot(h))
    n = np.cross(np.array([0,0,1]),h)
    #Semi-major axis
    a = 1/((2/rmag)-(vmag**2)/mu)

    evector = np.cross(v,h)/(mu) - r/rmag;
   #Eccentricity
    e = np.sqrt(evector.dot(evector))
    #Inclination
    inc = np.arccos(h[2]/hmag); 

    if np.dot(r,v) > 0:
        nu = np.arccos(np.dot(evector,r)/(e*rmag))
    else:
        nu = 2*np.pi - np.arccos(np.dot(evector,r)/(e*rmag))
    if evector[2] >= 0:
        omega = np.arccos(np.dot(n,evector)/(e*np.sqrt(n.dot(n))))
    else:
        omega = 2*np.pi - np.arccos(np.dot(n,evector)/(e*np.sqrt(n.dot(n))))
    if n[1] >= 0:
        Omega = np.arccos(n[0]/np.sqrt(n.dot(n)))
    else:
        Omega = 2*np.pi - np.arccos(n[0]/np.sqrt(n.dot(n)))

    TrueLongitude = nu + omega + Omega

Minor issue: when printing these elements while the simulation is running there does seem to be slight differences between the function above and what REBOUND outputs, I am curious what is causing these discrepancies? Below are outputs of Earth's semi-major axis, eccentricity and inclination both calculated from REBOUND while the sim is integrating compared to the output from the equations above.

image

hannorein commented 3 years ago

On why your own orbit routine gives different results:

SuperdoerTrav commented 3 years ago

On why your own orbit routine gives different results:

  • mu=1 is probably not exactly right. It should be mu = G*(m_0 + m_earth), where m_0 is the mass of the central object (heliocentric) or the mass of the interior centre of mass (Jacobi)
  • xyz and vxvyvz need to be relative to the central object (heliocentric) or the interior centre of mass (Jacobi)

Appreciate your taking the time to respond to this issue. The correct mu tweaks these orbital elements at around the 5th decimal. The coordinates I am feeding into these functions are those out put by REBOUND. By default the values should be relative to the centre of mass (Jacobi) correct?

hannorein commented 3 years ago

No the cartesian coordinates REBOUND outputs are in the inertial frame you work in. For example the centre of mass frame if you make a call to move_to_com().

SuperdoerTrav commented 3 years ago

I'm sure the issue is with a failure on my part to be in the correct coordinate system. To get the coordinates for the planets in the first initialization we've used the REBOUND get NASA ephemeris tool for a specific JD date. We then let the sim run, and after a certain time we request the serialized xyz vxvyvz coordinates for all particles in the sim. These coordinates are then in the center of mass frame. To get consistent Keplerian orbital elements we would need to transform back into the Jacobi coordinates of the system?

hannorein commented 3 years ago

I'm really not sure what could be the problem in your code. Maybe you can post a complete example with two or three particles that shows the issue? It might be a misunderstanding of some sort, maybe it's not clear what the orbital elements are that REBOUND returns, but I'm pretty confident that REBOUND itself should not be the problem here. The following will work:

The orbital elements from (A) and (B) will be the same. If you are careful about the output/input (or using the REBOUND binary files), then this will be exact. Otherwise there will be a floating point rounding error at 1e-16.

SuperdoerTrav commented 3 years ago

I will look to reproduce this issue in a complete example and send it. Since having this problem I have attempted to move forward with new simulation and I'm finding the simulation archive tool to not work as I thought it should either.

I've used the archive with an time interval to save and I find over a short integration when I reopen the binary file it works as intended and I may pick up where I left off. When I use this tool in a longer simulation such as over 100 time steps when I return to open the binary file it always says there is only one index in the file and that the simulation is at time of

  1. I've tried several ways of trouble shooting this but ended up just moving to single manual snapshots that I've been indexing myself. Though for a 500 particle system these binary files seems quite large ~ 500 mb?

Thank you, Travis

On Wed, Sep 15, 2021, 8:38 AM Hanno Rein @.***> wrote:

I'm really not sure what could be the problem in your code. Maybe you can post a complete example with two or three particles that shows the issue? It might be a misunderstanding of some sort, maybe it's not clear what the orbital elements are that REBOUND returns, but I'm pretty confident that REBOUND itself should not be the problem here. The following will work:

  • Setup a simulation in REBOUND
  • Integrate it forward in time
  • Output cartesian coordinates
  • Output orbital elements (A)
  • Create a new simulation
  • Add particles to new simulation using cartesian coordintes
  • Output orbital elements (B)

The orbital elements from (A) and (B) will be the same. If you are careful about the output/input (or using the REBOUND binary files), then this will be exact. Otherwise there will be a floating point rounding error at 1e-16.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/hannorein/rebound/issues/567#issuecomment-920132254, or unsubscribe https://github.com/notifications/unsubscribe-auth/AKLHCEFYCZYE7F6HFI7EMKTUCC4YXANCNFSM5DMUZLSQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

SuperdoerTrav commented 3 years ago

Sorry, I am in error on the file sizes. The 500 Mb is the total file size of all the manual snapshots I had until the point. The single snapshots are around 500 Kb.

On Thu, Sep 16, 2021 at 2:45 PM Travis Yeager @.***> wrote:

I will look to reproduce this issue in a complete example and send it. Since having this problem I have attempted to move forward with new simulation and I'm finding the simulation archive tool to not work as I thought it should either.

I've used the archive with an time interval to save and I find over a short integration when I reopen the binary file it works as intended and I may pick up where I left off. When I use this tool in a longer simulation such as over 100 time steps when I return to open the binary file it always says there is only one index in the file and that the simulation is at time of 0. I've tried several ways of trouble shooting this but ended up just moving to single manual snapshots that I've been indexing myself. Though for a 500 particle system these binary files seems quite large ~ 500 mb?

Thank you, Travis

On Wed, Sep 15, 2021, 8:38 AM Hanno Rein @.***> wrote:

I'm really not sure what could be the problem in your code. Maybe you can post a complete example with two or three particles that shows the issue? It might be a misunderstanding of some sort, maybe it's not clear what the orbital elements are that REBOUND returns, but I'm pretty confident that REBOUND itself should not be the problem here. The following will work:

  • Setup a simulation in REBOUND
  • Integrate it forward in time
  • Output cartesian coordinates
  • Output orbital elements (A)
  • Create a new simulation
  • Add particles to new simulation using cartesian coordintes
  • Output orbital elements (B)

The orbital elements from (A) and (B) will be the same. If you are careful about the output/input (or using the REBOUND binary files), then this will be exact. Otherwise there will be a floating point rounding error at 1e-16.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/hannorein/rebound/issues/567#issuecomment-920132254, or unsubscribe https://github.com/notifications/unsubscribe-auth/AKLHCEFYCZYE7F6HFI7EMKTUCC4YXANCNFSM5DMUZLSQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

SuperdoerTrav commented 3 years ago

I have reproduced this code and fed it new data and I am now getting entirely consistent results. I am unsure what I did to cause the initial error, my previous data sets do still produce incorrect results, though I do not believe this has anything to do with REBOUND. Thank you for all your time.

hannorein commented 3 years ago

I agree that this is most likely not related to REBOUND. But if you ever find out what the problem was, post it here. Maybe it is helpful to others too.