diff --git a/lib/LeapFrog.py b/lib/LeapFrog.py index 46a1046..982f25b 100644 --- a/lib/LeapFrog.py +++ b/lib/LeapFrog.py @@ -55,7 +55,7 @@ def leapfrog(dyn_syst, bin_syst, duration, dt, recover_param=False, display=Fals sma[j] = bin_syst.sma ecc[j] = bin_syst.ecc - if display and j % 50 == 0: + if display and j % 10 == 0: # display progression if len(dyn_syst.bodylist) == 1: d.on_running(dyn_syst, step=j, label="step {0:d}/{1:d}".format(j, N)) diff --git a/lib/hermite.py b/lib/hermite.py index fbcd85e..62022b4 100644 --- a/lib/hermite.py +++ b/lib/hermite.py @@ -102,7 +102,7 @@ def hermite(dyn_syst, duration, dt, recover_param=False, display=False, savename sma[j] = bin_syst.sma ecc[j] = bin_syst.ecc - if display and j % 100 == 0: + if display and j % 10 == 0: # display progression if len(dyn_syst.bodylist) == 1: d.on_running(dyn_syst, step=j, label="step {0:d}/{1:d}".format(j, N)) diff --git a/lib/objects.py b/lib/objects.py index aa516b7..33962b8 100755 --- a/lib/objects.py +++ b/lib/objects.py @@ -81,6 +81,14 @@ class System(Body): mass = mass + body.m return mass + @property + def mu(self): + prod = 1 + for body in self.bodylist: + prod = prod * body.m + mu = prod/self.M + return mu + @property def COM(self): #return center of mass in cartesian np_array coord = np.zeros(3) @@ -102,6 +110,31 @@ class System(Body): body.q = body.q - self.COM body.v = body.v - self.COMV + @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 + 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 + ECOM = self.mu/2.*np.linalg.norm(dv)**2 - Ga*self.M*self.mu/np.linalg.norm(dr) + return ECOM + @property def L(self): #return angular momentum of bodies in system L = np.zeros(3) @@ -122,18 +155,10 @@ class System(Body): E = T + W return E - @property - def mu(self): - prod = 1 - for body in self.bodylist: - prod = prod * body.m - mu = prod/self.M - return mu - @property def ecc(self): #exentricity of two body sub system if len(self.bodylist) == 2 : - ecc = (2.*self.E*(np.linalg.norm(self.L)**2))/((Ga**2)*(self.M**2)*(self.mu**3)) + 1. + ecc = (2.*self.ECOM*(np.linalg.norm(self.LCOM)**2))/((Ga**2)*(self.M**2)*(self.mu**3)) + 1. else : ecc = np.nan return ecc @@ -141,7 +166,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.E) + sma = -Ga*self.M*self.mu/(2.*self.ECOM) else : sma = np.nan return sma diff --git a/lib/plots.py b/lib/plots.py index b573d75..c9a727f 100755 --- a/lib/plots.py +++ b/lib/plots.py @@ -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))/au + 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() #Update data (with the new _and_ the old points) for i,body in enumerate(dyn_syst.bodylist): - x, y, z = body.q/au + x, y, z = body.q self.lines[i].set_data_3d([x], [y], [z]) if not label is None: if self.blackstyle: @@ -94,7 +94,7 @@ class DynamicUpdate(): #We need to draw *and* flush self.fig.canvas.draw() self.fig.canvas.flush_events() - if not step is None and step%100==0: + if not step is None and step%50==0: self.fig.savefig("tmp/{0:06d}.png".format(step),bbox_inches="tight") def close(self): diff --git a/main.py b/main.py index 87ca220..14504c5 100755 --- a/main.py +++ b/main.py @@ -31,7 +31,7 @@ def main(): step = np.sort(step)[::-1] integrator = "leapfrog" n_bodies = 2 - display = False + display = True savename = "{0:d}bodies_{1:s}".format(n_bodies, integrator) #simulation start @@ -42,7 +42,9 @@ def main(): dyn_syst = System(bodylist, main=True) E, L = [], [] - for step0 in step: + for i,step0 in enumerate(step): + 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) elif integrator.lower() in ['hermite','herm']: @@ -52,7 +54,6 @@ def main(): parameters = [duration, step, dyn_syst, integrator] display_parameters(E, L, sma, ecc, parameters=parameters, savename=savename) - print(sma,ecc) return 0 if __name__ == '__main__': diff --git a/plots/2bodies_leapfrog_a_e.png b/plots/2bodies_leapfrog_a_e.png index 0306913..bd08fa1 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_dynsyst.gif b/plots/2bodies_leapfrog_dynsyst.gif index a0fdd11..a781e09 100644 Binary files a/plots/2bodies_leapfrog_dynsyst.gif and b/plots/2bodies_leapfrog_dynsyst.gif differ