hannorein / rebound

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

Help to understand some changes in orbital elements #522

Closed gacarita closed 3 years ago

gacarita commented 3 years ago

I am currently running some initial condition grids for stability maps in resonant cases, in these cases I need to ensure that my "initial resonant angle" is librating around the angle I am studying. I'm setting the values of Omega omega and M and varying it, and to make sure that happens, but in a few cases the angles that I defined initially are not the ones that appear in the initial conditions I'm printing, the results at the end of the simulations seem to be correct because I am checking from the Cartesian coordinates, but I am not really understanding why the angles that I define initially are being redefined.

Would anyone know the reason for this? could this be happening due to the changing parameters of the orbital elements as I define the retrograde orbit (described here : https://rebound.readthedocs.io/en/latest/ipython/OrbitalElements.html), which could be confusing when looking at the initial conditions?

    sim.add(m=(m1),r=0.0046422) # Star
    sim.add(m=m2, a=1, e=0,r=0.00046607333333333)
    sim.move_to_com()
    sim.add(m=m3, a=a, e=e, inc=np.pi,Omega=0,omega=0,M=0)
    sim.integrator ="ias15"
    sim.dt = 0.01*sim.particles[1].P
    ps=sim.particles
    print("ptest : Omega,omega,M,inc ",ps[2].Omega*180/np.pi,ps[2].omega*180/np.pi,ps[2].M*180/np.pi,ps[2].inc*180/np.pi)
    print("planet : Omega,omega,M,inc ",ps[1].Omega*180/np.pi,ps[1].omega*180/np.pi,ps[1].M*180/np.pi)
    def phi_21(l,lp,pomega): return -q2*(l) - p2*(lp) + (p2+q2)*(pomega)
    print("ressonant angle:",phi_12(ps[2].l,ps[1].l,ps[2].pomega)*180/np.pi)

Print's


a,e : (1.5624388893445624, 0.0)
ptest : Omega,omega,M,inc  0.0 0.0 0.0 180.0
planet : Omega,omega,M,inc  0.0 180.0 180.0
ressonant angle: -1.5580271224017512e-30

a,e : (1.56354133377078, 0.0)
ptest : Omega,omega,M,inc  0.0 0.0 50.273536827068774 180.0
planet : Omega,omega,M,inc  0.0 180.0 180.0
ressonant angle: -1.5580271224017512e-30

a,e : (1.5646437781969975, 0.0)
ptest : Omega,omega,M,inc  0.0 -180.0 180.0 180.0
planet : Omega,omega,M,inc  0.0 180.0 180.0
ressonant angle: 540.0
hannorein commented 3 years ago

Some angles are degenerate for e=0.

hannorein commented 3 years ago

The coordinate transformations in REBOUND are written in a way that everything is still self consistent even in those cases. But you should not expect to get back the omega you use to initialize an orbit if the eccentricity is 0 (because it's not well defined).

gacarita commented 3 years ago

I see.

Thank you @hannorein.