diff --git a/lib/objects.py b/lib/objects.py index 33962b8..6d7db93 100755 --- a/lib/objects.py +++ b/lib/objects.py @@ -101,7 +101,7 @@ class System(Body): def COMV(self): #return center of mass velocity in cartesian np_array coord = np.zeros(3) for body in self.bodylist: - coord = coord + body.v + coord = coord + body.m*body.v coord = coord/self.M return coord @@ -113,25 +113,15 @@ class System(Body): @property def LCOM(self): #return angular momentum of the center of mass LCOM = np.zeros(3) - dr = np.zeros(3) - dv = np.zeros(3) - for body in self.bodylist: - for otherbody in self.bodylist: - if body != otherbody: - dr = body.q-otherbody.q - dv = body.v-otherbody.v + 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) return LCOM @property def ECOM(self): #return mechanical energy of the center of mass - dr = np.zeros(3) - dv = np.zeros(3) - for body in self.bodylist: - for otherbody in self.bodylist: - if body != otherbody: - dr = body.q-otherbody.q - dv = body.v-otherbody.v + 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) return ECOM diff --git a/main.py b/main.py index 14504c5..fbf0b0b 100755 --- a/main.py +++ b/main.py @@ -16,33 +16,33 @@ def main(): e = np.array([0., 0., 1./4.]) # Eccentricity psi = np.array([0., 0., 0.])*np.pi/180. # Inclination of the orbital plane in degrees - x1 = np.array([0., -1., 0.])*a[0] - x2 = np.array([0., 1., 0.])*a[1] - x3 = np.array([np.cos(psi[2]), 0., np.sin(psi[2])])*a[2] + x1 = np.array([0., -1., 0.])*a[0]*(1.+e[0]) + x2 = np.array([0., 1., 0.])*a[1]*(1.+e[1]) + x3 = np.array([np.cos(psi[2]), 0., np.sin(psi[2])])*a[2]*(1.+e[2]) q = np.array([x1, x2, x3]) - v1 = np.array([np.sqrt(Ga*m[1]**2/((m[0]+m[1])*np.sqrt(np.sum((q[0]-q[1])**2)))), 0., 0.]) - v2 = np.array([-np.sqrt(Ga*m[0]**2/((m[0]+m[1])*np.sqrt(np.sum((q[0]-q[1])**2)))), 0., 0.]) + 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.]) + 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.]) v3 = np.array([0., np.sqrt(Ga*(m[0]+m[1])*(2./np.sqrt(np.sum(q[2]**2))-1./a[2])), 0.]) v = np.array([v1, v2, v3]) #integration parameters duration, step = 100*yr/yr, np.array([1./(365.25*2.), 1./(365.25*1.), 5./(365.25*1.)])*yr/yr #integration time and step in years - step = np.sort(step)[::-1] + step = np.sort(step) integrator = "leapfrog" - n_bodies = 2 - display = True - savename = "{0:d}bodies_{1:s}".format(n_bodies, integrator) + n_bodies = 3 + display = False + savename = "{0:d}bodies_mass_{1:s}".format(n_bodies, integrator) #simulation start - bodylist = [] - for i in range(n_bodies): - bodylist.append(Body(m[i], q[i], v[i])) - bin_syst = System(bodylist[0:2]) - dyn_syst = System(bodylist, main=True) - E, L = [], [] for i,step0 in enumerate(step): + bodylist = [] + for j in range(n_bodies): + bodylist.append(Body(m[j], q[j], v[j])) + bin_syst = System(bodylist[0:2]) + dyn_syst = System(bodylist, main=True) + if i != 0: display = False if integrator.lower() in ['leapfrog', 'frogleap', 'frog']: diff --git a/plots/2bodies_leapfrog_a_e.png b/plots/2bodies_leapfrog_a_e.png index bd08fa1..76d3b93 100644 Binary files a/plots/2bodies_leapfrog_a_e.png and b/plots/2bodies_leapfrog_a_e.png differ diff --git a/plots/2bodies_leapfrog_dEm.png b/plots/2bodies_leapfrog_dEm.png index d13cda3..a2e9517 100644 Binary files a/plots/2bodies_leapfrog_dEm.png and b/plots/2bodies_leapfrog_dEm.png differ diff --git a/plots/2bodies_leapfrog_dL2.png b/plots/2bodies_leapfrog_dL2.png index 096e0d1..cab3beb 100644 Binary files a/plots/2bodies_leapfrog_dL2.png and b/plots/2bodies_leapfrog_dL2.png differ diff --git a/plots/2bodies_leapfrog_dynsyst.gif b/plots/2bodies_leapfrog_dynsyst.gif index a781e09..0807821 100644 Binary files a/plots/2bodies_leapfrog_dynsyst.gif and b/plots/2bodies_leapfrog_dynsyst.gif differ diff --git a/plots/3bodies_leapfrog_a_e.png b/plots/3bodies_leapfrog_a_e.png new file mode 100644 index 0000000..add8bb5 Binary files /dev/null and b/plots/3bodies_leapfrog_a_e.png differ diff --git a/plots/3bodies_leapfrog_dEm.png b/plots/3bodies_leapfrog_dEm.png index fc05309..02652f8 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 ba84c86..d44a3d1 100644 Binary files a/plots/3bodies_leapfrog_dL2.png and b/plots/3bodies_leapfrog_dL2.png differ diff --git a/plots/3bodies_leapfrog_dynsyst.gif b/plots/3bodies_leapfrog_dynsyst.gif index fde4db1..8a579dd 100644 Binary files a/plots/3bodies_leapfrog_dynsyst.gif and b/plots/3bodies_leapfrog_dynsyst.gif differ diff --git a/plots/3bodies_mass_leapfrog_a_e.png b/plots/3bodies_mass_leapfrog_a_e.png new file mode 100644 index 0000000..89e4df4 Binary files /dev/null and b/plots/3bodies_mass_leapfrog_a_e.png differ diff --git a/plots/3bodies_mass_leapfrog_dEm.png b/plots/3bodies_mass_leapfrog_dEm.png new file mode 100644 index 0000000..b18b912 Binary files /dev/null and b/plots/3bodies_mass_leapfrog_dEm.png differ diff --git a/plots/3bodies_mass_leapfrog_dL2.png b/plots/3bodies_mass_leapfrog_dL2.png new file mode 100644 index 0000000..7a0dd99 Binary files /dev/null and b/plots/3bodies_mass_leapfrog_dL2.png differ diff --git a/plots/3bodies_psi80deg_leapfrog_a_e.png b/plots/3bodies_psi80deg_leapfrog_a_e.png new file mode 100644 index 0000000..34c8b94 Binary files /dev/null and b/plots/3bodies_psi80deg_leapfrog_a_e.png differ diff --git a/plots/3bodies_psi80deg_leapfrog_dEm.png b/plots/3bodies_psi80deg_leapfrog_dEm.png new file mode 100644 index 0000000..1d04692 Binary files /dev/null and b/plots/3bodies_psi80deg_leapfrog_dEm.png differ diff --git a/plots/3bodies_psi80deg_leapfrog_dL2.png b/plots/3bodies_psi80deg_leapfrog_dL2.png new file mode 100644 index 0000000..2c19732 Binary files /dev/null and b/plots/3bodies_psi80deg_leapfrog_dL2.png differ