get back to full dimensions
@@ -26,10 +26,10 @@ class Body:
|
||||
self.vp = np.zeros(3,dtype=np.longdouble)
|
||||
|
||||
def __repr__(self): # Called upon "print(body)"
|
||||
return r"Body of mass: {0:.1e} $M_\odot$, position: {1}, velocity: {2}".format(self.m, self.q, self.v)
|
||||
return r"Body of mass: {0:.1e} $M_\odot$, position: {1}, velocity: {2}".format(self.m/Ms, self.q, self.v)
|
||||
|
||||
def __str__(self): # Called upon "str(body)"
|
||||
return r"Body of mass: {0:.1e} $M_\odot$".format(self.m)
|
||||
return r"Body of mass: {0:.1e} $M_\odot$".format(self.m/Ms)
|
||||
|
||||
@property
|
||||
def p(self):
|
||||
|
||||
@@ -75,12 +75,12 @@ class DynamicUpdate():
|
||||
|
||||
def on_running(self, dyn_syst, step=None, label=None):
|
||||
xdata, ydata, zdata = dyn_syst.get_positions()
|
||||
values = np.sqrt(np.sum((np.array((xdata,ydata,zdata))**2).T,axis=1))
|
||||
values = np.sqrt(np.sum((np.array((xdata,ydata,zdata))**2).T,axis=1))/au
|
||||
self.min_x, self.max_x = -np.max([np.abs(values).max(),self.max_x]), np.max([np.abs(values).max(),self.max_x])
|
||||
self.set_lims()
|
||||
#Update data (with the new _and_ the old points)
|
||||
for i,body in enumerate(dyn_syst.bodylist):
|
||||
x, y, z = body.q
|
||||
x, y, z = body.q/au
|
||||
self.lines[i].set_data_3d([x], [y], [z])
|
||||
if not label is None:
|
||||
if self.blackstyle:
|
||||
|
||||
@@ -8,4 +8,4 @@ globals()['G'] = 6.67e-11 #Gravitational constant in SI units
|
||||
globals()['Ms'] = 2e30 #Solar mass in kg
|
||||
globals()['au'] = 1.5e11 #Astronomical unit in m
|
||||
globals()['yr'] = 3.15576e7 #year in seconds
|
||||
globals()['Ga'] = G*Ms/au**3 #Gravitational constant dimensionless
|
||||
globals()['Ga'] = G#*Ms/au**3 #Gravitational constant dimensionless
|
||||
14
main.py
@@ -11,14 +11,10 @@ from lib.units import *
|
||||
|
||||
def main():
|
||||
#initialisation
|
||||
m = np.array([1., 1., 1e-1],dtype=np.longdouble)*Ms/Ms # Masses in Solar mass
|
||||
a = np.array([1., 1., 5.],dtype=np.longdouble)*au/au # Semi-major axis in astronomical units
|
||||
e = np.array([0., 0., 0.],dtype=np.longdouble) # Eccentricity
|
||||
<<<<<<< HEAD
|
||||
psi = np.array([0., 0., 80.],dtype=np.longdouble)*np.pi/180. # Inclination of the orbital plane in degrees
|
||||
=======
|
||||
psi = np.array([0., 0., 0.],dtype=np.longdouble)*np.pi/180. # Inclination of the orbital plane in degrees
|
||||
>>>>>>> 22fa187 (add Energy display)
|
||||
m = np.array([1., 1., 1e-1],dtype=np.longdouble)*Ms#/Ms # Masses in Solar mass
|
||||
a = np.array([1., 1., 10.],dtype=np.longdouble)*au#/au # Semi-major axis in astronomical units
|
||||
e = np.array([0., 0., 0.25],dtype=np.longdouble) # Eccentricity
|
||||
psi = np.array([0., 0., 60.],dtype=np.longdouble)*np.pi/180. # Inclination of the orbital plane in degrees
|
||||
|
||||
x1 = np.array([0., -1., 0.],dtype=np.longdouble)*a[0]*(1.+e[0])
|
||||
x2 = np.array([0., 1., 0.],dtype=np.longdouble)*a[1]*(1.+e[1])
|
||||
@@ -31,7 +27,7 @@ def main():
|
||||
v = np.array([v1, v2, v3],dtype=np.longdouble)
|
||||
|
||||
#integration parameters
|
||||
duration, step = 100*yr, np.array([600000.],dtype=np.longdouble) #integration time and step in seconds
|
||||
duration, step = 1000*yr, np.array([10.*86400.],dtype=np.longdouble) #integration time and step in seconds
|
||||
step = np.sort(step)[::-1]
|
||||
integrator = "leapfrog"
|
||||
n_bodies = 3
|
||||
|
||||
BIN
plots/3bodies_leapfrog_E.png
Normal file
|
After Width: | Height: | Size: 177 KiB |
|
Before Width: | Height: | Size: 42 KiB After Width: | Height: | Size: 62 KiB |
|
Before Width: | Height: | Size: 69 KiB After Width: | Height: | Size: 113 KiB |
|
Before Width: | Height: | Size: 68 KiB After Width: | Height: | Size: 72 KiB |