Question: Additional Forces

Hello,

I'm trying to run some simulations by adding the entire spherical harmonics simulations, but I'm encountering some problems with the program using REBOUND. I have tested my code with a simple RK45, and it works fine. However, now I'm trying to implement this using REBOUND, and I'm sure that I'm doing something wrong so far.


sim = rebound.Simulation()
date = "2023-07-24 00:00"
sim.units = ('s', 'm', 'kg')

x0 = 6478.136000000000*1e3
y0 =-2.194237229530863e-13*1e3
z0 = 1.340073885192083e-15*1e3
vx0 =-1.833551089721986e-17*1e3
vy0 = 7.844114000000000*1e3
vz0 = 6.845279345644194e-12*1e3

sim.add("399",date=date)
sim.status()
sim.move_to_com()
sim.add(m=0.0,x=x0,y=y0,z=z0,vx=vx0,vy=vy0,vz=vz0) # initial condition on CM
ps = sim.particles
sim.integrator = 'IAS15'
x = np.zeros((1,Noutputs))
y = np.zeros((1,Noutputs))
z = np.zeros((1,Noutputs))
t0 = []
w_earth = 7.292115e-5

def harm(reb_sim):
    sim = reb_sim.contents

    vec_r = [sim.particles[1].x,sim.particles[1].y,sim.particles[1].z] # meters
    gst = w_earth * sim.t
    rotat_vect = mxv(M1, vec_r)

    ACG = [0, 0, 0]
    ACG = egm2008_kuga(radius[0], gm[0], c, s, nm, gst, rotat_vect)
    #print(ACG,sim.t)
    sim.particles[1].ax += ACG[0]
    sim.particles[1].ay += ACG[1]
    sim.particles[1].az += ACG[2]

sim.additional_forces = harm
sim.force_is_position_dependent = 1

a_s = np.zeros(len(times))
for i,time in enumerate(times):
    sim.integrate(time)
    t0.append(sim.t)
    x[0][i] = ps[1].x
    y[0][i] = ps[1].y
    z[0][i] = ps[1].z
    a_s[i] = sim.particles[1].ax

I've integrated this orbit around Earth for 1 day, but it escapes and at the end 1e6 meters distance from the Earth. I suspect that the acceleration is updating in every single step of the integrator, but I only want it to update in the successful steps.

I appreciate any help.

hannorein / rebound

Question: Additional Forces #712