Custom operator on eccentricity and semimajor axis

[huhell]

Hello,
I am trying to add tidal dissipation as a custom operator, which would update the eccentricity, and semimajor axis as e += dedt * dt and a += dadt * dt.
I use Equations (1) and (2) from Jackson et al. (2008), ApJ, 678, 1396. Here is my code

sim = rebound.Simulation()
sim.units = ('yr', 'AU', 'Msun')
sim.integrator = 'whfast'
sim.add(m=1.0)
sim.add(m=constants.Mjup/constants.Msun, a=a0 / constants.AU, e=e0)
sim.move_to_com()
ps = sim.particles

def dissipation(reb_sim, rebx_operator, dt):
        sim = reb_sim.contents
        dissip = rebx_operator.contents
        ps = sim.particles

        Rp = 1.0
        Qp = 1e4
        Rs = 1.0
        Qs = 1e6

        e1 = (63.0/4.0) * (constants.GG * (ps[0].m*constants.Msun)**3)**0.5 * ((Rp*constants.Rjup)**5 / (Qp * ps[1].m*constants.Mjup))
        e2 = (171.0/16.0) * (constants.GG / (ps[0].m*constants.Msun))**0.5 * ((Rs*constants.Rsun)**5 * ps[1].m*constants.Mjup / Qs)
        de_dt = -(e1 + e2) * (ps[1].a*constants.AU)**(-13.0/2.0) * e

        a1 = (63.0/2.0) * (constants.GG * (ps[0].m*constants.Msun)**3)**0.5 * ((Rp*constants.Rjup)**5 / (Qp * ps[1].m*constants.Mjup))
        a2 = (9.0/2.0) * (constants.GG / (ps[0].m*constants.Msun))**0.5 * ((Rs*constants.Rsun)**5 * ps[1].m*constants.Mjup / Qs)
        da_dt = -(a1*ps[1].e**2 + a2) * (ps[1].a*constants.AU)**(-11.0/2.0)

        #ps[1].e += de_dt * dt
        #ps[1].a += (da_dt * dt) / constants.AU

rebx = reboundx.Extras(sim)
myop = rebx.create_operator("dissipation")
myop.operator_type = "updater"
myop.step_function = dissipation
rebx.add_operator(myop)

But when I integrate using sim.integrate(), I get the following error:

Traceback (most recent call last):
  File "_ctypes/callbacks.c", line 315, in 'calling callback function'
  File "test.py", line 143, in dissipation
    ps[1].a += (da_dt * dt) / constants.AU
AttributeError: can't set attribute

Is it impossible to update the eccentricity and semimajor axis ?

Thanks a lot !

[dtamayo]

Hi,

Unfortunately it's not well defined to just set the eccentricity or semimajor axis (you have to choose what to do with the remaining elements), so REBOUND doesn't let you do that.

If you just want exponential damping of the eccentricity and semimajor axis, you can just use modify_orbits_forces that's already implemented, and just set the appropriate tau_a and tau_e using your tidal parameters (tau_e = 1/de/dt).

If that's not good enough, you're doing the right thing by implementing an operator, but you have to set xyzzy and vxvyvz for an operator. One option would be to follow appendix A of Lee and Peale 2002.

[huhell]

Hi,
Thank you very much for the clarification !