phi derivation

lib/objects.py

@@ -5,6 +5,8 @@ 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 *
 
@@ -203,3 +205,11 @@ class System(Body):
         else :
             sma = np.nan
         return sma
+
+    def phi(self,body1,body2): #return angle in degree between plans formed by body1 and body2 trajectories
+        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
+        return phi

main.py

@@ -14,7 +14,7 @@ def main():
     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., 45.],dtype=np.longdouble)*np.pi/180.    # Inclination of the orbital plane in degrees
+    psi = np.array([0., 0., 80.],dtype=np.longdouble)*np.pi/180.    # Inclination of the orbital plane in degrees
 
     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])
@@ -27,7 +27,7 @@ def main():
     v = np.array([v1, v2, v3],dtype=np.longdouble)
 
     #integration parameters
-    duration, step = 10*yr, np.array([60000.],dtype=np.longdouble) #integration time and step in seconds
+    duration, step = 100*yr, np.array([600000.],dtype=np.longdouble) #integration time and step in seconds
     step = np.sort(step)[::-1]
     integrator = "leapfrog"
     n_bodies = 3