Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
r2 = dx*dx + dy*dy + dz*dz
E_pot -= particles[i].m*particles[j].m/np.sqrt(r2)
return E_kin+E_pot
times = np.logspace(np.log10(100.*dt),np.log10(tmax),Ngrid)
if integrator=="wh" or integrator=="mercury":
move_to_heliocentric()
else:
rebound.move_to_center_of_momentum()
ei = energy()
es = []
runtime = 0.
for t in times:
rebound.integrate(t,exactFinishTime=0)
ef = energy()
e = np.fabs((ei-ef)/ei)+1.1e-16
es.append(e)
runtime += rebound.get_timing()
integrator, run, trial = par
print integrator.ljust(13) + " %9.5fs"%(runtime) + "\t Error: %e" %( e)
es = np.array(es)
return [times, es]
# Setup particles (data taken from NASA Horizons)
# This could also be easily read in from a file.
rebound.add_particle( m=1.00000597682, x=-4.06428567034226e-3, y=-6.08813756435987e-3, z=-1.66162304225834e-6, vx=+6.69048890636161e-6, vy=-6.33922479583593e-6, vz=-3.13202145590767e-9 ) # Sun
rebound.add_particle( m=1./1047.355, x=+3.40546614227466e+0, y=+3.62978190075864e+0, z=+3.42386261766577e-2, vx=-5.59797969310664e-3, vy=+5.51815399480116e-3, vz=-2.66711392865591e-6 ) # Jupiter
rebound.add_particle( m=1./3501.6, x=+6.60801554403466e+0, y=+6.38084674585064e+0, z=-1.36145963724542e-1, vx=-4.17354020307064e-3, vy=+3.99723751748116e-3, vz=+1.67206320571441e-5 ) # Saturn
rebound.add_particle( m=1./22869., x=+1.11636331405597e+1, y=+1.60373479057256e+1, z=+3.61783279369958e-1, vx=-3.25884806151064e-3, vy=+2.06438412905916e-3, vz=-2.17699042180559e-5 ) # Uranus
rebound.add_particle( m=1./19314., x=-3.01777243405203e+1, y=+1.91155314998064e+0, z=-1.53887595621042e-1, vx=-2.17471785045538e-4, vy=-3.11361111025884e-3, vz=+3.58344705491441e-5 ) # Neptune
# Set the center of momentum to be at the origin
rebound.move_to_center_of_momentum()
rebound.set_integrator("mercury")
#rebound.set_integrator("mikkola")
rebound.set_dt(40.)
rebound.integrate(1e7)
print(rebound.get_timing())
print(rebound.status())
rebound.integrate(2e7)
print(rebound.get_timing())
print(rebound.status())
E_pot -= G*particles[i].m*particles[j].m/np.sqrt(r2)
return E_kin+E_pot
times = np.logspace(np.log10(100.*dt),np.log10(tmax),1000)
if integrator=="wh" or integrator=="mercury":
move_to_heliocentric()
else:
rebound.move_to_center_of_momentum()
ei = energy()
es = []
timing = 0.
for t in times:
rebound.integrate(t,exactFinishTime=0)
ef = energy()
e = np.fabs((ei-ef)/ei)+1.1e-16
es.append(e)
timing += rebound.get_timing()
es = np.array(es)
print integrator + " done. %.5fs"%(timing)
return [times, es]
# Move particles so that the center of mass is (and stays) at the origin
rebound.move_to_center_of_momentum()
# Integrate until t=100 (roughly 16 orbits)
rebound.integrate(100.)
# Modify particles
# As an example, we are reverting the velocities
particles = rebound.particles_get()
for i in range(rebound.get_N()):
particles[i].vx *= -1.
particles[i].vy *= -1.
particles[i].vz *= -1.
# Integrate another 100 time units, until t=200
rebound.integrate(200.)
# Get particles back and print positions
# Since we're integrating forward and then backward we end up with the
# particles exactly where they started out from (note that we moved to the
# center of momentum frame)
for i in range(rebound.get_N()):
print(particles[i].x, particles[i].y, particles[i].z)
# Set variables (defaults are G=1, t=0, dt=0.01)
#rebound.set_G(1.)
#rebound.set_t(0.)
#rebound.set_dt(0.01)
# Add particles
# All parameters omitted are set to 0 by default.
rebound.particle_add( m=1. ) # Star
rebound.particle_add( m=1e-3, x=1., vy=1. ) # Planet
# Move particles so that the center of mass is (and stays) at the origin
rebound.move_to_center_of_momentum()
# Integrate until t=100 (roughly 16 orbits)
rebound.integrate(100.)
# Modify particles
# As an example, we are reverting the velocities
particles = rebound.particles_get()
for i in range(rebound.get_N()):
particles[i].vx *= -1.
particles[i].vy *= -1.
particles[i].vz *= -1.
# Integrate another 100 time units, until t=200
rebound.integrate(200.)
# Get particles back and print positions
# Since we're integrating forward and then backward we end up with the
# particles exactly where they started out from (note that we moved to the
# center of momentum frame)
e = 1.-pow(10.,e)
dt = pow(10.,dt)*torb
rebound.reset()
rebound.integrator = integrator
rebound.force_is_velocitydependent = 0
rebound.dt = dt
rebound.add(m=1.)
rebound.add(m=0., x=(1.-e), vy=np.sqrt((1.+e)/(1.-e)))
particles = rebound.particles
Ei = -1./np.sqrt(particles[1].x*particles[1].x+particles[1].y*particles[1].y+particles[1].z*particles[1].z) + 0.5 * (particles[1].vx*particles[1].vx+particles[1].vy*particles[1].vy+particles[1].vz*particles[1].vz)
rebound.integrate(tmax,exactFinishTime=0,keepSynchronized=1)
Ef = -1./np.sqrt(particles[1].x*particles[1].x+particles[1].y*particles[1].y+particles[1].z*particles[1].z) + 0.5 * (particles[1].vx*particles[1].vx+particles[1].vy*particles[1].vy+particles[1].vz*particles[1].vz)
return [float(rebound.iter)/rebound.t*dt, np.fabs((Ef-Ei)/Ei)+1e-16, rebound.timing/rebound.t*dt*1e6/2., (Ef-Ei)/Ei]
rebound.dt = dt
rebound.add(m=1.)
rebound.add(m=0.,a=1.,e=e0)
#rebound.move_to_com()
#rebound.init_megno(1.e-16)
particles = rebound.particles
def starkforce(): # need to put inside simulation(par) to have access to S and particles
particles[1].ax += -S
rebound.additional_forces = starkforce
rebound.integrate(50000.*np.pi)
return [rebound.megno, rebound.t]
dx = particles[i].x-particles[j].x
dy = particles[i].y-particles[j].y
dz = particles[i].z-particles[j].z
r2 = dx*dx + dy*dy + dz*dz
E_pot -= G*particles[i].m*particles[j].m/np.sqrt(r2)
return E_kin+E_pot
if integrator=="wh" or integrator=="mercury" or integrator[0:7]=="swifter":
move_to_heliocentric()
else:
rebound.move_to_com()
ei = energy()
runtime = 0.
rebound.integrate(tmax,exact_finish_time=0)
ef = energy()
e = np.fabs((ei-ef)/ei)+1.1e-16
runtime += rebound.timing
integrator, dt, run = par
print integrator.ljust(13) + " %9.5fs"%(runtime) + "\t Error: %e" %( e)
return [runtime, e]