Created
April 29, 2021 20:54
-
-
Save alekseygenerozov/3cbfe6ad0a65c5e98eee216ae65b2e63 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import rebound | |
import reboundx | |
%pylab inline | |
import astropy.constants as const | |
##Constants | |
G=const.G.cgs.value | |
M_sun=const.M_sun.cgs.value | |
c=const.c.cgs.value | |
pc=const.pc.cgs.value | |
##Converting speed of light to code units. Units are pc, solar masses, and G=1 | |
cc=c/(G*M_sun/pc)**0.5 | |
##Initializing simulation. | |
end_r=array([ 2.43429405e-05, -6.34018455e-05, -2.87889574e-04]) | |
end_v=array([-43247.28772886, 42257.14066544, 152805.55464892]) | |
sim=rebound.Simulation() | |
sim.add(m=4e6) | |
sim.add(m=100.05, x=end_r[0], y=end_r[1], z=end_r[2], vx=end_v[0], vy=end_v[1], vz=end_v[2]) | |
sim.move_to_com() | |
rebx = reboundx.Extras(sim) | |
gr=rebx.add("gr") | |
gr.params["c"]=cc | |
##Various quantities of interest. Rg is the gravitational radius of the central black hole. Rp is the pericenter | |
rg=sim.particles[0].m/cc**2. | |
print("Rg: {0}".format(rg)) | |
print("Initial radius/Rg: {0}".format(np.linalg.norm(end_r/rg))) | |
orbs=sim.calculate_orbits(primary=sim.particles[0]) | |
aa=orbs[0].a | |
rp=(1.-orbs[0].e)*aa | |
print("Initial sma:",aa) | |
print("Rp/Rg:",rp/rg) | |
##Integrate simulation | |
p_in=2.0*np.pi*(aa**3.0/sim.particles[0].m)**0.5 | |
sim.integrate(p_in) | |
##Final semimajor axis: negative! | |
orbs=sim.calculate_orbits(primary=sim.particles[0]) | |
print("Final sma:",orbs[0].a) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment