add phi display

lib/LeapFrog.py

@@ -47,6 +47,7 @@ def leapfrog(dyn_syst, bin_syst, duration, dt, recover_param=False, display=Fals
     L = np.zeros((N, 3),dtype=np.longdouble)
     sma = np.zeros(N,dtype=np.longdouble)
     ecc = np.zeros(N,dtype=np.longdouble)
+    phi = np.zeros(N,dtype=np.longdouble)
     for j in range(N):
         LP(dyn_syst,dt)
 
@@ -54,6 +55,7 @@ def leapfrog(dyn_syst, bin_syst, duration, dt, recover_param=False, display=Fals
         L[j] = dyn_syst.L
         sma[j] = bin_syst.sma
         ecc[j] = bin_syst.ecc
+        phi[j] = dyn_syst.phi
 
         if display and j % 10 == 0:
             if gif:
@@ -72,4 +74,4 @@ def leapfrog(dyn_syst, bin_syst, duration, dt, recover_param=False, display=Fals
             system("convert tmp/temp.gif -fuzz 10% -layers Optimize plots/{0:s}_dynsyst.gif".format(savename))
 
     if recover_param:
-        return E, L, sma, ecc
+        return E, L, sma, ecc, phi

lib/hermite.py

@@ -93,6 +93,7 @@ def hermite(dyn_syst, bin_syst, duration, dt, recover_param=False, display=False
     L = np.zeros((N, 3),dtype=np.longdouble)
     sma = np.zeros(N,dtype=np.longdouble)
     ecc = np.zeros(N,dtype=np.longdouble)
+    phi = np.zeros(N,dtype=np.longdouble)
 
     for j in range(N):
         HPC(dyn_syst, dt)
@@ -101,6 +102,7 @@ def hermite(dyn_syst, bin_syst, duration, dt, recover_param=False, display=False
         L[j] = dyn_syst.L
         sma[j] = bin_syst.sma
         ecc[j] = bin_syst.ecc
+        phi[j] = dyn_syst.phi
 
         if display and j % 10 == 0:
             if gif:
@@ -119,4 +121,4 @@ def hermite(dyn_syst, bin_syst, duration, dt, recover_param=False, display=False
             system("convert tmp/temp.gif -fuzz 10% -layers Optimize plots/{0:s}_dynsyst.gif".format(savename))
 
     if recover_param:
-        return E, L, sma, ecc
+        return E, L, sma, ecc, phi

lib/objects.py

@@ -5,8 +5,6 @@ Class definition for physical attribute
 """
 from os import system
 import numpy as np
-from astropy.coordinates import Angle
-from astropy import units as u
 from lib.plots import DynamicUpdate
 from lib.units import *
 
@@ -191,10 +189,6 @@ class System(Body):
         E = T + W
         return E
 
-
-
-
-
     @property
     def ecc(self): #exentricity of two body sub system
         if len(self.bodylist) == 2 :
@@ -212,15 +206,13 @@ class System(Body):
         return sma
 
     @property
-    def phi(self,body1,body2): #return angle in degree between plans formed by body1 and body2 (perurbator) trajectories
+    def phi(self): #return angle in degree between plans formed by body1 and body2 (perurbator) trajectories
         if len(self.bodylist) == 3 :
             body1 = self.bodylist[0]
             body2 = self.bodylist[2]
             n1 = np.cross(body1.q, body1.v)
             n2 = np.cross(body2.q, body2.v)
-            phi = np.arccos(np.dot(n1, n2) / (np.linalg.norm(n1) * np.linalg.norm(n2)))
-            phi = Angle(phi, u.radian)
-            phi = phi.dec
+            phi = np.arccos(np.dot(n1, n2) / (np.linalg.norm(n1) * np.linalg.norm(n2)))*180./np.pi
         else :
             phi = np.nan
         return phi

lib/plots.py

@@ -101,7 +101,7 @@ class DynamicUpdate():
         plt.close()
 
 
-def display_parameters(E,L,sma,ecc,parameters,savename=""):
+def display_parameters(E,L,sma,ecc,phi,parameters,savename=""):
     """
     """
     if savename != "":
@@ -150,5 +150,14 @@ def display_parameters(E,L,sma,ecc,parameters,savename=""):
     fig4.suptitle("Mechanical energy of the whole system "+title2)
     fig4.savefig("plots/{0:s}E.png".format(savename),bbox_inches="tight")
     
+    fig5 = plt.figure(figsize=(15,7))
+    ax5 = fig5.add_subplot(111)
+    for i in range(len(phi)):
+        ax5.plot(np.arange(phi[i].shape[0])*step[-1]/yr, phi[i], label="step of {0:.2e}s".format(step[i]))
+    ax5.set(xlabel=r"$t \, [yr]$", ylabel=r"$\phi \, [^{\circ}]$")
+    ax5.legend()
+    fig5.suptitle("Inclination angle of the perturbator's orbital plane "+title2)
+    fig5.savefig("plots/{0:s}phi.png".format(savename),bbox_inches="tight")
+    
     plt.show(block=True)
   

main.py

@@ -31,12 +31,12 @@ def main():
     step = np.sort(step)[::-1]
     integrator = "leapfrog"
     n_bodies = 3
-    display = True
+    display = False
     gif = False
     savename = "{0:d}bodies_{1:s}".format(n_bodies, integrator)
 
     #simulation start
-    E, L, ecc, sma = [], [], [], []
+    E, L, ecc, sma, phi = [], [], [], [], []
     for i,step0 in enumerate(step):
         bodylist = []
         for j in range(n_bodies):
@@ -47,16 +47,17 @@ def main():
         if i != 0:
             display = False
         if integrator.lower() in ['leapfrog', 'frogleap', 'frog']:
-            E0, L0, sma0, ecc0 = leapfrog(dyn_syst, bin_syst, duration, step0, recover_param=True, display=display, savename=savename, gif=gif)
+            E0, L0, sma0, ecc0, phi0 = leapfrog(dyn_syst, bin_syst, duration, step0, recover_param=True, display=display, savename=savename, gif=gif)
         elif integrator.lower() in ['hermite','herm']:
-            E0, L0, sma0, ecc0 = hermite(dyn_syst, bin_syst, duration, step0, recover_param=True, display=display, savename=savename, gif=gif)
+            E0, L0, sma0, ecc0, phi0 = 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)
+        phi.append(phi0)
 
     parameters = [duration, step, dyn_syst, integrator]
-    display_parameters(E, L, sma, ecc, parameters=parameters, savename=savename)
+    display_parameters(E, L, sma, ecc, phi, parameters=parameters, savename=savename)
     return 0
 
 if __name__ == '__main__':