diff --git a/lib/objects.py b/lib/objects.py index d0e2893..2883b12 100755 --- a/lib/objects.py +++ b/lib/objects.py @@ -14,6 +14,8 @@ class Body: self.m = mass self.q = position self.v = velocity + self.qb = position + self.vb = velocity self.a = np.zeros(3,dtype=np.longdouble) self.ap = np.zeros(3,dtype=np.longdouble) self.j = np.zeros(3,dtype=np.longdouble) @@ -31,6 +33,10 @@ class Body: def p(self): return self.v*self.m + @property + def pb(self): + return self.vb*self.m + class System(Body): def __init__(self, bodylist, main = False, blackstyle=True): @@ -49,25 +55,25 @@ class System(Body): def __str__(self): # Called upon "str(system)" return str([str(body) for body in self.bodylist]) - @property + def get_masses(self): #return the masses of each object return np.array([body.m for body in self.bodylist],dtype=np.longdouble) - - @property + + def get_positions(self): #return the positions of the bodies xdata = np.array([body.q[0] for body in self.bodylist],dtype=np.longdouble) ydata = np.array([body.q[1] for body in self.bodylist],dtype=np.longdouble) zdata = np.array([body.q[2] for body in self.bodylist],dtype=np.longdouble) return xdata, ydata, zdata - - @property + + def get_velocities(self): #return the positions of the bodies vxdata = np.array([body.v[0] for body in self.bodylist],dtype=np.longdouble) vydata = np.array([body.v[1] for body in self.bodylist],dtype=np.longdouble) vzdata = np.array([body.v[2] for body in self.bodylist],dtype=np.longdouble) return vxdata, vydata, vzdata - - @property + + def get_momenta(self): #return the momenta of the bodies pxdata = np.array([body.p[0] for body in self.bodylist],dtype=np.longdouble) pydata = np.array([body.p[1] for body in self.bodylist],dtype=np.longdouble) @@ -106,9 +112,40 @@ class System(Body): return coord def COMShift(self): #Shift coordinates of bodies in system to COM frame and set COM at rest + COM = self.COM + COMV = self.COMV for body in self.bodylist: - body.q = body.q - self.COM - body.v = body.v - self.COMV + body.q = body.q - COM + body.v = body.v - COMV + + def COMShiftBin(self): #Shift coordinates of inner binary system to COM frame and set COM at rest + COM = self.COM + COMV = self.COMV + for body in self.bodylist: + body.qb = body.qb - COM + body.vb = body.vb - COMV + + @property + def LBIN(self): #return angular momentum of inner binary + self.COMShiftBin() + L = np.zeros(3,dtype=np.longdouble) + for body in self.bodylist: + L = L + np.cross(body.qb,body.pb) + return L + + @property + def EBIN(self): #return total energy of inner binary + self.COMShiftBin() + T = 0 + W = 0 + for body in self.bodylist: + T = T + 1./2.*body.m*np.linalg.norm(body.vb)**2 + for otherbody in self.bodylist: + if body != otherbody: + rij = np.linalg.norm(body.qb-otherbody.qb) + W = W - Ga*body.m*otherbody.m/rij + E = T + W + return E @property def LCOM(self): #return angular momentum of the center of mass @@ -116,6 +153,9 @@ class System(Body): dr = self.bodylist[0].m/self.mu*self.bodylist[0].q dv = self.bodylist[0].m/self.mu*self.bodylist[0].v LCOM = self.mu*np.cross(dr,dv) + + LCOM = self.L + return LCOM @property @@ -123,6 +163,9 @@ class System(Body): dr = self.bodylist[0].m/self.mu*self.bodylist[0].q dv = self.bodylist[0].m/self.mu*self.bodylist[0].v ECOM = self.mu/2.*np.linalg.norm(dv)**2 - Ga*self.M*self.mu/np.linalg.norm(dr) + + ECOM = self.E + return ECOM @property @@ -148,7 +191,7 @@ class System(Body): @property def ecc(self): #exentricity of two body sub system if len(self.bodylist) == 2 : - ecc = (2.*self.ECOM*(np.linalg.norm(self.LCOM)**2))/((Ga**2)*(self.M**2)*(self.mu**3)) + 1. + ecc = (2.*self.EBIN*(np.linalg.norm(self.LBIN)**2))/((Ga**2)*(self.M**2)*(self.mu**3)) + 1. else : ecc = np.nan return ecc @@ -156,7 +199,7 @@ class System(Body): @property def sma(self): #semi major axis of two body sub system if len(self.bodylist) == 2 : - sma = -Ga*self.M*self.mu/(2.*self.ECOM) + sma = -Ga*self.M*self.mu/(2.*self.EBIN) else : sma = np.nan return sma diff --git a/lib/plots.py b/lib/plots.py index 29f793f..5a04a8b 100755 --- a/lib/plots.py +++ b/lib/plots.py @@ -74,7 +74,7 @@ class DynamicUpdate(): self.ax.set_zlabel('AU') def on_running(self, dyn_syst, step=None, label=None): - xdata, ydata, zdata = dyn_syst.get_positions + xdata, ydata, zdata = dyn_syst.get_positions() values = np.sqrt(np.sum((np.array((xdata,ydata,zdata))**2).T,axis=1)) 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() diff --git a/main.py b/main.py index a1245c2..9d9ce57 100755 --- a/main.py +++ b/main.py @@ -14,7 +14,7 @@ def main(): 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 - psi = np.array([0., 0., 0.],dtype=np.longdouble)*np.pi/180. # Inclination of the orbital plane in degrees + psi = np.array([0., 0., 45.],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]) @@ -27,16 +27,16 @@ def main(): v = np.array([v1, v2, v3],dtype=np.longdouble) #integration parameters - duration, step = 100*yr, np.array([60.],dtype=np.longdouble) #integration time and step in seconds + duration, step = 10*yr, np.array([60000.],dtype=np.longdouble) #integration time and step in seconds step = np.sort(step)[::-1] integrator = "leapfrog" n_bodies = 3 - display = False + display = True gif = False savename = "{0:d}bodies_{1:s}".format(n_bodies, integrator) #simulation start - E, L = [], [] + E, L, ecc, sma = [], [], [], [] for i,step0 in enumerate(step): bodylist = [] for j in range(n_bodies): @@ -47,11 +47,13 @@ def main(): if i != 0: display = False if integrator.lower() in ['leapfrog', 'frogleap', 'frog']: - E0, L0, sma, ecc = leapfrog(dyn_syst, bin_syst, duration, step0, recover_param=True, display=display, savename=savename, gif=gif) + E0, L0, sma0, ecc0 = leapfrog(dyn_syst, bin_syst, duration, step0, recover_param=True, display=display, savename=savename, gif=gif) elif integrator.lower() in ['hermite','herm']: - E0, L0, sma, ecc = hermite(dyn_syst, bin_syst, duration, step0, recover_param=True, display=display, savename=savename, gif=gif) + E0, L0, sma0, ecc0 = hermite(dyn_syst, bin_syst, duration, step0, recover_param=True, display=display, savename=savename, gif=gif) E.append(E0) L.append(L0) + ecc.append(ecc0) + sma.append(sma0) parameters = [duration, step, dyn_syst, integrator] display_parameters(E, L, sma, ecc, parameters=parameters, savename=savename) diff --git a/plots/3bodies_leapfrog_a_e.png b/plots/3bodies_leapfrog_a_e.png index 6e9f87f..4bff667 100644 Binary files a/plots/3bodies_leapfrog_a_e.png and b/plots/3bodies_leapfrog_a_e.png differ diff --git a/plots/3bodies_leapfrog_dEm.png b/plots/3bodies_leapfrog_dEm.png index 00d9102..d0db2d3 100644 Binary files a/plots/3bodies_leapfrog_dEm.png and b/plots/3bodies_leapfrog_dEm.png differ diff --git a/plots/3bodies_leapfrog_dL2.png b/plots/3bodies_leapfrog_dL2.png index f5f5ff2..4bfd99a 100644 Binary files a/plots/3bodies_leapfrog_dL2.png and b/plots/3bodies_leapfrog_dL2.png differ