debugg hermite and force use of longdouble type

lib/LeapFrog.py

@@ -18,7 +18,7 @@ def Drift(dyn_syst, dt):
 
 def Kick(dyn_syst, dt):
     for body in dyn_syst.bodylist:
-        body.a = np.zeros(3)
+        body.a = np.zeros(3,dtype=np.longdouble)
         for otherbody in dyn_syst.bodylist:
             if body != otherbody:
                 rij = np.linalg.norm(body.q - otherbody.q)
@@ -43,10 +43,10 @@ def leapfrog(dyn_syst, bin_syst, duration, dt, recover_param=False, display=Fals
         d.launch(dyn_syst.blackstyle)
 
     N = np.ceil(duration / dt).astype(int)
-    E = np.zeros(N)
-    L = np.zeros((N, 3))
-    sma = np.zeros(N)
-    ecc = np.zeros(N)
+    E = np.zeros(N,dtype=np.longdouble)
+    L = np.zeros((N, 3),dtype=np.longdouble)
+    sma = np.zeros(N,dtype=np.longdouble)
+    ecc = np.zeros(N,dtype=np.longdouble)
     for j in range(N):
         LP(dyn_syst,dt)
 
@@ -58,9 +58,9 @@ def leapfrog(dyn_syst, bin_syst, duration, dt, recover_param=False, display=Fals
         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))
+                d.on_running(dyn_syst, step=j, label="{0:.2f} years".format(j*dt))
             else:
-                d.on_running(dyn_syst, step=j, label="step {0:d}/{1:d}".format(j, N))
+                d.on_running(dyn_syst, step=j, label="{0:.2f} years".format(j*dt))
     if display:
         d.close()
         if not savename is None:

lib/hermite.py

@@ -13,23 +13,23 @@ from lib.units import *
 
 def Update_a(dyn_syst):  # update acceleration of bodies in system
     for body in dyn_syst.bodylist:
-        body.a = np.zeros(3)
+        body.a = np.zeros(3,dtype=np.longdouble)
         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) * G * otherbody.m / (rij ** 3)
+                body.a = body.a - (body.q - otherbody.q) * Ga * otherbody.m / (rij ** 3)
 
 
 def Update_j(dyn_syst):  # update jerk of bodies in system
     for body in dyn_syst.bodylist:
-        body.j = np.zeros(3)
+        body.j = np.zeros(3,dtype=np.longdouble)
         for otherbody in dyn_syst.bodylist:
             if body != otherbody:
                 rij = np.linalg.norm(body.q - otherbody.q)
                 deltav = (body.v - otherbody.v)
                 deltar = (body.q - otherbody.q)
                 vr = deltav + 3. * deltar * np.inner(deltav, deltar) / (rij ** 2)
-                body.j = body.j - G * otherbody.m / (rij ** 3) * vr
+                body.j = body.j - Ga * otherbody.m / (rij ** 3) * vr
 
 
 def Predict(dyn_syst, dt):  # update predicted position and velocities of bodies in system
@@ -40,23 +40,23 @@ def Predict(dyn_syst, dt):  # update predicted position and velocities of bodies
 
 def Update_ap(dyn_syst):  # update acceleration of bodies in system
     for body in dyn_syst.bodylist:
-        body.ap = np.zeros(3)
+        body.ap = np.zeros(3,dtype=np.longdouble)
         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) * G * otherbody.m / (rij ** 3)
+                body.ap = body.ap - (body.qp - otherbody.qp) * Ga * otherbody.m / (rij ** 3)
 
 
 def Update_jp(dyn_syst):  # update jerk of bodies in system
     for body in dyn_syst.bodylist:
-        body.jp = np.zeros(3)
+        body.jp = np.zeros(3,dtype=np.longdouble)
         for otherbody in dyn_syst.bodylist:
             if body != otherbody:
                 rij = np.linalg.norm(body.qp - otherbody.qp)
                 deltav = (body.vp - otherbody.vp)
                 deltar = (body.qp - otherbody.qp)
                 vr = deltav + 3. * deltar * np.inner(deltav, deltar) / (rij ** 2)
-                body.jp = body.jp - G * otherbody.m / (rij ** 3) * vr
+                body.jp = body.jp - Ga * otherbody.m / (rij ** 3) * vr
 
 
 def Correct(dyn_syst, dt):  # correct position and velocities of bodies in system
@@ -69,7 +69,7 @@ def Correct(dyn_syst, dt):  # correct position and velocities of bodies in syste
 
 
 def HPC(dyn_syst, dt):  # update position and velocities of bodies in system with hermite predictor corrector
-    COMShift(dyn_syst)
+    dyn_syst.COMShift()
     Update_a(dyn_syst)
     Update_j(dyn_syst)
     Predict(dyn_syst, dt)
@@ -79,7 +79,7 @@ def HPC(dyn_syst, dt):  # update position and velocities of bodies in system wit
     dyn_syst.time = dyn_syst.time + dt
 
 
-def hermite(dyn_syst, duration, dt, recover_param=False, display=False, savename=None):
+def hermite(dyn_syst, bin_syst, duration, dt, recover_param=False, display=False, savename=None):
     if display:
         try:
             system("mkdir tmp")
@@ -89,10 +89,10 @@ def hermite(dyn_syst, duration, dt, recover_param=False, display=False, savename
         d.launch(dyn_syst.blackstyle)
 
     N = np.ceil(duration / dt).astype(int)
-    E = np.zeros(N)
-    L = np.zeros((N, 3))
-    sma = np.zeros(N)
-    ecc = np.zeros(N)
+    E = np.zeros(N,dtype=np.longdouble)
+    L = np.zeros((N, 3),dtype=np.longdouble)
+    sma = np.zeros(N,dtype=np.longdouble)
+    ecc = np.zeros(N,dtype=np.longdouble)
 
     for j in range(N):
         HPC(dyn_syst, dt)
@@ -105,9 +105,9 @@ def hermite(dyn_syst, duration, dt, recover_param=False, display=False, savename
         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))
+                d.on_running(dyn_syst, step=j, label="{0:.2f} yearq".format(j*dt))
             else:
-                d.on_running(dyn_syst, step=j, label="step {0:d}/{1:d}".format(j, N))
+                d.on_running(dyn_syst, step=j, label="{0:.2f} years".format(j*dt))
     if display:
         d.close()
         if not savename is None:

lib/objects.py

@@ -14,12 +14,12 @@ class Body:
         self.m = mass
         self.q = position
         self.v = velocity
-        self.a = np.zeros(3)
-        self.ap = np.zeros(3)
-        self.j = np.zeros(3)
-        self.jp = np.zeros(3)
-        self.qp = np.zeros(3)
-        self.vp = np.zeros(3)
+        self.a = np.zeros(3,dtype=np.longdouble)
+        self.ap = np.zeros(3,dtype=np.longdouble)
+        self.j = np.zeros(3,dtype=np.longdouble)
+        self.jp = np.zeros(3,dtype=np.longdouble)
+        self.qp = np.zeros(3,dtype=np.longdouble)
+        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)
@@ -51,27 +51,27 @@ class System(Body):
 
     @property
     def get_masses(self): #return the masses of each object
-        return np.array([body.m for body in self.bodylist])
+        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])
-        ydata = np.array([body.q[1] for body in self.bodylist])
-        zdata = np.array([body.q[2] for body in self.bodylist])
+        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])
-        vydata = np.array([body.v[1] for body in self.bodylist])
-        vzdata = np.array([body.v[2] for body in self.bodylist])
+        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])
-        pydata = np.array([body.p[1] for body in self.bodylist])
-        pzdata = np.array([body.p[2] for body in self.bodylist])
+        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)
+        pzdata = np.array([body.p[2] for body in self.bodylist],dtype=np.longdouble)
         return pxdata, pydata, pzdata
 
     @property
@@ -91,7 +91,7 @@ class System(Body):
 
     @property
     def COM(self): #return center of mass in cartesian np_array
-        coord = np.zeros(3)
+        coord = np.zeros(3,dtype=np.longdouble)
         for body in self.bodylist:
             coord = coord + body.m*body.q
         coord = coord/self.M
@@ -99,7 +99,7 @@ class System(Body):
 
     @property
     def COMV(self): #return center of mass velocity in cartesian np_array
-        coord = np.zeros(3)
+        coord = np.zeros(3,dtype=np.longdouble)
         for body in self.bodylist:
             coord = coord + body.m*body.v
         coord = coord/self.M
@@ -112,7 +112,7 @@ class System(Body):
 
     @property
     def LCOM(self): #return angular momentum of the center of mass
-        LCOM = np.zeros(3)
+        LCOM = np.zeros(3,dtype=np.longdouble)
         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)
@@ -127,7 +127,7 @@ class System(Body):
 
     @property
     def L(self): #return angular momentum of bodies in system
-        L = np.zeros(3)
+        L = np.zeros(3,dtype=np.longdouble)
         for body in self.bodylist:
             L = L + np.cross(body.q,body.p)
         return L

lib/plots.py

@@ -115,7 +115,7 @@ def display_parameters(E,L,sma,ecc,parameters,savename=""):
     fig1 = plt.figure(figsize=(15,7))
     ax1 = fig1.add_subplot(111)
     for i in range(len(E)):
-        ax1.plot(np.arange(E[i].shape[0])*step[i], np.abs((E[i]-E[i][0])/E[i][0]), label="step of {0:.2e}yr".format(step[i]))
+        ax1.plot(np.arange(E[i].shape[0]-1)*step[i], np.abs((E[i][1:]-E[i][0])/E[i][0]), label="step of {0:.2e}yr".format(step[i]))
     ax1.set(xlabel=r"$t \, [yr]$", ylabel=r"$\left|\frac{\delta E_m}{E_m(t=0)}\right|$", yscale='log')
     ax1.legend()
     fig1.suptitle(title1.format("mechanical energy")+title2)
@@ -126,7 +126,7 @@ def display_parameters(E,L,sma,ecc,parameters,savename=""):
     for i in range(len(L)):
         dL = ((L[i]-L[i][0])/L[i][0])
         dL[np.isnan(dL)] = 0.
-        ax2.plot(np.arange(L[i].shape[0])*step[i], np.abs(np.sum(dL,axis=1)), label="step of {0:.2e}yr".format(step[i]))
+        ax2.plot(np.arange(L[i].shape[0]-1)*step[i], np.abs(np.sum(dL[1:],axis=1)), label="step of {0:.2e}yr".format(step[i]))
     ax2.set(xlabel=r"$t \, [yr]$", ylabel=r"$\left|\frac{\delta \vec{L}}{\vec{L}(t=0)}\right|$",yscale='log')
     ax2.legend()
     fig2.suptitle(title1.format("kinetic moment")+title2)

lib/units.py

@@ -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*yr**2/au**3 #Gravitational constant adimensionned
+globals()['Ga'] = G*Ms*yr**2/au**3 #Gravitational constant dimensionless

main.py

@@ -11,28 +11,28 @@ from lib.units import *
 
 def main():
     #initialisation
-    m = np.array([1., 1., 1e-5])*Ms/Ms  # Masses in Solar mass
-    a = np.array([1., 1., 5.])*au/au   # Semi-major axis in astronomical units
-    e = np.array([0., 0., 1./4.])   # Eccentricity
-    psi = np.array([0., 0., 0.])*np.pi/180.    # Inclination of the orbital plane in degrees
+    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
 
-    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])
+    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])
+    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.])
-    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])
+    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)
+    v = np.array([v1, v2, v3],dtype=np.longdouble)
 
     #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)
-    integrator = "leapfrog"
+    duration, step = 100*yr/yr, np.array([1./(365.25*2.), 1./(365.25*1.), 5./(365.25*1.)],dtype=np.longdouble)*yr/yr #integration time and step in years
+    step = np.sort(step)[::-1]
+    integrator = "hermite"
     n_bodies = 3
     display = False
-    savename = "{0:d}bodies_mass_{1:s}".format(n_bodies, integrator)
+    savename = "{0:d}bodies_{1:s}".format(n_bodies, integrator)
 
     #simulation start
     E, L = [], []