diff --git a/lib/LeapFrog.py b/lib/LeapFrog.py index fb6570f..7a6b423 100644 --- a/lib/LeapFrog.py +++ b/lib/LeapFrog.py @@ -22,7 +22,7 @@ def Kick(dyn_syst, dt): for otherbody in dyn_syst.bodylist: if body != otherbody: rij = np.linalg.norm(body.q - otherbody.q) - body.a = body.a - (body.q - otherbody.q) * Ga * otherbody.m / (rij ** 3) + body.a = body.a - (body.q - otherbody.q) * G * otherbody.m / (rij ** 3) body.v = body.v + dt * body.a diff --git a/lib/hermite.py b/lib/hermite.py index 73ce822..994898e 100644 --- a/lib/hermite.py +++ b/lib/hermite.py @@ -17,7 +17,7 @@ def Update_a(dyn_syst): # update acceleration of bodies in system for otherbody in dyn_syst.bodylist: if body != otherbody: rij = np.linalg.norm(body.q - otherbody.q) - body.a = body.a - (body.q - otherbody.q) * Ga * otherbody.m / (rij ** 3) + body.a = body.a - (body.q - otherbody.q) * G * otherbody.m / (rij ** 3) def Update_j(dyn_syst): # update jerk of bodies in system @@ -29,7 +29,7 @@ def Update_j(dyn_syst): # update jerk of bodies in system deltav = (body.v - otherbody.v) deltar = (body.q - otherbody.q) vr = deltav + 3. * deltar * np.inner(deltav, deltar) / (rij ** 2) - body.j = body.j - Ga * otherbody.m / (rij ** 3) * vr + body.j = body.j - G * otherbody.m / (rij ** 3) * vr def Predict(dyn_syst, dt): # update predicted position and velocities of bodies in system @@ -44,7 +44,7 @@ def Update_ap(dyn_syst): # update acceleration of bodies in system for otherbody in dyn_syst.bodylist: if body != otherbody: rij = np.linalg.norm(body.qp - otherbody.qp) - body.ap = body.ap - (body.qp - otherbody.qp) * Ga * otherbody.m / (rij ** 3) + body.ap = body.ap - (body.qp - otherbody.qp) * G * otherbody.m / (rij ** 3) def Update_jp(dyn_syst): # update jerk of bodies in system @@ -56,7 +56,7 @@ def Update_jp(dyn_syst): # update jerk of bodies in system deltav = (body.vp - otherbody.vp) deltar = (body.qp - otherbody.qp) vr = deltav + 3. * deltar * np.inner(deltav, deltar) / (rij ** 2) - body.jp = body.jp - Ga * otherbody.m / (rij ** 3) * vr + body.jp = body.jp - G * otherbody.m / (rij ** 3) * vr def Correct(dyn_syst, dt): # correct position and velocities of bodies in system diff --git a/lib/objects.py b/lib/objects.py index 6f90f68..9f8517e 100755 --- a/lib/objects.py +++ b/lib/objects.py @@ -137,7 +137,7 @@ class System(Body): for otherbody in self.bodylist: if body != otherbody: rij = np.linalg.norm(body.qb-otherbody.qb) - W = W - Ga*body.m*otherbody.m/rij + W = W - G*body.m*otherbody.m/rij E = T + W return E @@ -150,7 +150,7 @@ class System(Body): for otherbody in self.bodylist: if body != otherbody: rij = np.linalg.norm(body.q-otherbody.q) - W = W - Ga*otherbody.m*body.m/(2.*rij) + W = W - G*otherbody.m*body.m/(2.*rij) E = T + W return E @@ -171,7 +171,7 @@ class System(Body): for otherbody in self.bodylist: if body != otherbody: rij = np.linalg.norm(body.q-otherbody.q) - W = W - Ga*otherbody.m*body.m/(2.*rij) + W = W - G*otherbody.m*body.m/(2.*rij) E = T + W return E @@ -185,7 +185,7 @@ class System(Body): @property def eccCOM(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.ECOM*(np.linalg.norm(self.LCOM)**2))/(G**2*self.M**2*self.mu**3) + 1. else : ecc = np.nan @@ -194,7 +194,7 @@ class System(Body): @property def smaCOM(self): #semi major axis of two body sub system if len(self.bodylist) == 2 : - sma = -Ga*self.mu*self.bodylist[0].m/(2.*self.ECOM) + sma = -G*self.mu*self.bodylist[0].m/(2.*self.ECOM) else : sma = np.nan return sma @@ -202,7 +202,7 @@ class System(Body): @property def ecc(self): #exentricity of two body sub system if len(self.bodylist) == 2 : - ecc = (2.*self.EBIN*(np.linalg.norm(self.LBIN)**2))/(Ga**2*self.M**2*self.mu**3) + 1. + ecc = (2.*self.EBIN*(np.linalg.norm(self.LBIN)**2))/(G**2*self.M**2*self.mu**3) + 1. else : ecc = np.nan return ecc @@ -210,7 +210,7 @@ class System(Body): @property def sma(self): #semi major axis of two body sub system if len(self.bodylist) == 2 : - sma = -Ga*self.mu*self.bodylist[0].m/(2.*self.EBIN) + sma = -G*self.mu*self.bodylist[0].m/(2.*self.EBIN) else : sma = np.nan return sma diff --git a/lib/units.py b/lib/units.py index 2256156..e23b866 100644 --- a/lib/units.py +++ b/lib/units.py @@ -7,5 +7,4 @@ Units used in the project. 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 \ No newline at end of file +globals()['yr'] = 3.15576e7 #year in seconds \ No newline at end of file diff --git a/main.py b/main.py index 99805dd..4f8c251 100755 --- a/main.py +++ b/main.py @@ -22,9 +22,9 @@ def main(): x3 = np.array([np.cos(psi[2]), 0., np.sin(psi[2])],dtype=np.longdouble)*a[2]*(1.+e[2]) q = np.array([x1, x2, x3],dtype=np.longdouble) - v1 = np.array([np.sqrt(Ga*m[0]*m[1]/((m[0]+m[1])*np.sqrt(np.sum((q[0]-q[1])**2)))), 0., 0.],dtype=np.longdouble) - v2 = np.array([-np.sqrt(Ga*m[0]*m[1]/((m[0]+m[1])*np.sqrt(np.sum((q[0]-q[1])**2)))), 0., 0.],dtype=np.longdouble) - v3 = np.array([0., np.sqrt(Ga*(m[0]+m[1])*(2./np.sqrt(np.sum(q[2]**2))-1./a[2])), 0.],dtype=np.longdouble) + v1 = np.array([np.sqrt(G*m[0]*m[1]/((m[0]+m[1])*np.sqrt(np.sum((q[0]-q[1])**2)))), 0., 0.],dtype=np.longdouble) + v2 = np.array([-np.sqrt(G*m[0]*m[1]/((m[0]+m[1])*np.sqrt(np.sum((q[0]-q[1])**2)))), 0., 0.],dtype=np.longdouble) + v3 = np.array([0., np.sqrt(G*(m[0]+m[1])*(2./np.sqrt(np.sum(q[2]**2))-1./a[2])), 0.],dtype=np.longdouble) v = np.array([v1, v2, v3],dtype=np.longdouble) #integration parameters