Fixed J2 and J4 bug

MintonGroup · May 15, 2023 · 172f437 · 172f437
1 parent 6f724c1
commit 172f437
Show file tree

Hide file tree

Showing 3 changed files with 16 additions and 7 deletions.
diff --git a/python/swiftest/swiftest/init_cond.py b/python/swiftest/swiftest/init_cond.py
@@ -277,8 +277,8 @@ def vec2xr(param: Dict, **kwargs: Any):
     space_coords = np.array(["x","y","z"])
 
     vector_vars = ["rh","vh","Ip","rot"]
-    scalar_vars = ["name","a","e","inc","capom","omega","capm","Gmass","radius","rhill","J2","J4"]
-    time_vars =  ["rh","vh","Ip","rot","a","e","inc","capom","omega","capm","Gmass","radius","rhill","J2","J4"]
+    scalar_vars = ["name","a","e","inc","capom","omega","capm","Gmass","radius","rhill","j2rp2","j4rp4"]
+    time_vars =  ["rh","vh","Ip","rot","a","e","inc","capom","omega","capm","Gmass","radius","rhill","j2rp2","j4rp4"]
 
     # Check for valid keyword arguments
     kwargs = {k:kwargs[k] for k,v in kwargs.items() if v is not None}

diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py
@@ -2223,9 +2223,9 @@ def add_solar_system_body(self,
             values = list(np.hsplit(np.array(body_list[0],dtype=np.dtype(object)),17))
         else:
             values = list(np.squeeze(np.hsplit(np.array(body_list,np.dtype(object)),17)))
-        keys = ["id","name","a","e","inc","capom","omega","capm","rh","vh","Gmass","radius","rhill","Ip","rot","J2","J4"]
+        keys = ["id","name","a","e","inc","capom","omega","capm","rh","vh","Gmass","radius","rhill","Ip","rot","j2rp2","j4rp4"]
         kwargs = dict(zip(keys,values))
-        scalar_floats = ["a","e","inc","capom","omega","capm","Gmass","radius","rhill","J2","J4"]
+        scalar_floats = ["a","e","inc","capom","omega","capm","Gmass","radius","rhill","j2rp2","j4rp4"]
         vector_floats = ["rh","vh","Ip","rot"]
         scalar_ints = ["id"]
 
@@ -2434,7 +2434,7 @@ def add_body(self,
         Hill's radius values if these are massive bodies
         rot: (3) or (n,3) array-like of float, optional
         Rotation rate vectors if these are massive bodies with rotation enabled.
-        Ip: (3) or (n,3) array-like of flaot, optional
+        Ip: (3) or (n,3) array-like of float, optional
         Principal axes moments of inertia vectors if these are massive bodies with rotation enabled.
 
         Returns
@@ -2451,7 +2451,7 @@ def input_to_array(val,t,n=None):
             elif t == "i":
                 t = np.int64
             elif t == "s":
-                t = np.str
+                t = str
 
             if val is None:
                 return None, n
@@ -2551,7 +2551,7 @@ def input_to_array_3d(val,n=None):
                 Gmass = self.GU * mass
 
         dsnew = init_cond.vec2xr(self.param, name=name, a=a, e=e, inc=inc, capom=capom, omega=omega, capm=capm, id=id,
-                                 Gmass=Gmass, radius=radius, rhill=rhill, Ip=Ip, rh=rh, vh=vh,rot=rot, J2=J2, J4=J4, time=time)
+                                 Gmass=Gmass, radius=radius, rhill=rhill, Ip=Ip, rh=rh, vh=vh,rot=rot, j2rp2=J2, j4rp4=J4, time=time)
 
         dsnew = self._combine_and_fix_dsnew(dsnew)
         self.save(verbose=False)

diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90
@@ -2930,6 +2930,15 @@ module subroutine swiftest_io_read_in_system(self, nc, param)
       if (.not.param%loblatecb) then
          if (allocated(self%pl%aobl)) deallocate(self%pl%aobl)
          if (allocated(self%tp%aobl)) deallocate(self%tp%aobl)
+      else
+         if (self%pl%nbody > 0) then
+            if (.not. allocated(self%pl%aobl)) allocate(self%pl%aobl(NDIM,self%pl%nbody))
+            self%pl%aobl(:,:) = 0.0_DP
+         end if
+         if (self%tp%nbody > 0) then
+            if (.not. allocated(self%tp%aobl)) allocate(self%tp%aobl(NDIM,self%tp%nbody))
+            self%tp%aobl(:,:) = 0.0_DP
+         end if
       end if
 
       return