Made sure rotation vectors were in deg/TU instead of rad/TU

examples/Fragmentation/Fragmentation_Movie.py

@@ -31,6 +31,7 @@
 import xarray as xr
 import matplotlib.pyplot as plt
 import matplotlib.animation as animation
+from scipy.spatial.transform import Rotation as R
 
 # ----------------------------------------------------------------------------------------------------------------------
 # Define the names and initial conditions of the various fragmentation simulation types
@@ -74,14 +75,14 @@
                                               np.array([ 0.05,  6.18,     0.0])] 
                }
 
-rot_vectors = {"disruption_headon"         : [np.array([0.0, 0.0, 0.0]),
-                                              np.array([0.0, 0.0, 0.0])],
-               "disruption_off_axis":        [np.array([0.0, 0.0, 0.0]),
-                                              np.array([0.0, 0.0, 0.0])],
-               "supercatastrophic_headon":   [np.array([0.0, 0.0, 0.0]),
-                                              np.array([0.0, 0.0, 0.0])],
-               "supercatastrophic_off_axis": [np.array([0.0, 0.0, 0.0]),
-                                              np.array([0.0, 0.0, 0.0])],
+rot_vectors = {"disruption_headon"         : [np.array([1.0e3, 0.0, 1.0e2]),
+                                              np.array([-1.0e2, -1.0e3, 0.0])],
+               "disruption_off_axis":        [np.array([0.0, 0.0, 2.0e3]),
+                                              np.array([0.0, 0.0, -1.0e4])],
+               "supercatastrophic_headon":   [np.array([3.0e2, 3.0e2, 1.0e3]),
+                                              np.array([1.0e3, -3.0e2, -5.0e3])],
+               "supercatastrophic_off_axis": [np.array([3.0e2, 3.0e2, 1.0e3]),
+                                              np.array([1.0e3, 3.0e3, -1.0e4])],
                "hitandrun_disrupt"         : [np.array([0.0, 0.0, 6.0e3]),
                                               np.array([0.0, 0.0, 1.0e4])],
                "hitandrun_pure"            : [np.array([0.0, 0.0, 6.0e3]),
@@ -129,7 +130,7 @@ def encounter_combiner(sim):
     """
 
     # Only keep a minimal subset of necessary data from the simulation and encounter datasets
-    keep_vars = ['name','rh','vh','Gmass','radius']
+    keep_vars = ['name','rh','vh','Gmass','radius','rot']
     data = sim.data[keep_vars]
     enc = sim.encounters[keep_vars].load()
 
@@ -158,6 +159,7 @@ def encounter_combiner(sim):
     # Add a bit of padding to the time, otherwise there are some issues with the interpolation in the last few frames.
     smooth_time = np.linspace(start=tgood.isel(time=0), stop=ds.time[-1], num=int(1.2*num_movie_frames))
     ds = ds.interp(time=smooth_time)
+    ds['rotangle'] = xr.zeros_like(ds['rot'])
 
     return ds
 
@@ -167,7 +169,6 @@ class AnimatedScatter(object):
     def __init__(self, sim, animfile, title, style, nskip=1):
 
         self.ds = encounter_combiner(sim)
-        self.sim = sim
         self.title = title
         self.body_color_list = {'Initial conditions': 'xkcd:windows blue',
                       'Disruption': 'xkcd:baby poop',
@@ -215,10 +216,6 @@ def update_plot(self, frame):
 
         # Define a function to calculate a reference frame for the animation
         # This will be based on the initial velocity of the Target body
-        def reference_frame(r_ref, v_ref, t):
-            coord_pos = r_ref + v_ref * t 
-            return coord_pos.values[0], coord_pos.values[1]
-
         def center(Gmass, x, y):
             x = x[~np.isnan(x)]
             y = y[~np.isnan(y)]
@@ -242,11 +239,14 @@ def data_stream(self, frame=0):
             Gmass = ds['Gmass'].values
             rh = ds['rh'].values
             point_rad = radius * self.ax_pt_size
-
-            # Save the initial velocity of body 1 to use as a reference
-            if frame == 0:
-                self.r_ref = ds.sel(name="Target")['rh'] 
-                self.v_ref = ds.sel(name="Target")['vh'] 
+            if frame > 0:
+                dsold = self.ds.isel(time=frame-1)
+                told = dsold['time'].values[()]
+                dt = t - told
+                ds['rotangle'] = dsold['rotangle'] + dt * dsold['rot']
+                ds['rotangle'] = np.mod(ds['rotangle'], 360.0)
+                self.ds[dict(time=frame)]['rotangle'] = ds['rotangle']
+
             yield t, Gmass, rh, point_rad
 
 if __name__ == "__main__":

python/swiftest/swiftest/init_cond.py

@@ -96,16 +96,16 @@ def solar_system_horizons(name: str,
     }
 
     planetrot = {
-       'Sun' : np.longdouble(2*np.pi / 25.05) / swiftest.JD2S, # Approximate
-       'Mercury': np.longdouble(2*np.pi / 58.646) / swiftest.JD2S,
-       'Venus': np.longdouble(2*np.pi / 243.0226 ) / swiftest.JD2S,
-       'Earth': np.longdouble(2*np.pi / 0.99726968) / swiftest.JD2S,
-       'Mars': np.longdouble(2*np.pi / 1.025957) / swiftest.JD2S,
-       'Jupiter': np.longdouble(2*np.pi / (9.9250 / 24.0) ) / swiftest.JD2S,
-       'Saturn': np.longdouble(2*np.pi / (10.656 / 24.0) ) / swiftest.JD2S,
-       'Uranus': np.longdouble(2*np.pi / 0.71833) / swiftest.JD2S,
-       'Neptune': np.longdouble(2*np.pi / 0.6713) / swiftest.JD2S,
-       'Pluto': np.longdouble(2*np.pi / 6.387230) / swiftest.JD2S
+       'Sun' : np.longdouble(360.0 / 25.05) / swiftest.JD2S, # Approximate
+       'Mercury': np.longdouble(360.0 / 58.646) / swiftest.JD2S,
+       'Venus': np.longdouble(360.0 / 243.0226 ) / swiftest.JD2S,
+       'Earth': np.longdouble(360.0 / 0.99726968) / swiftest.JD2S,
+       'Mars': np.longdouble(360.0 / 1.025957) / swiftest.JD2S,
+       'Jupiter': np.longdouble(360.0 / (9.9250 / 24.0) ) / swiftest.JD2S,
+       'Saturn': np.longdouble(360.0 / (10.656 / 24.0) ) / swiftest.JD2S,
+       'Uranus': np.longdouble(360.0 / 0.71833) / swiftest.JD2S,
+       'Neptune': np.longdouble(360.0 / 0.6713) / swiftest.JD2S,
+       'Pluto': np.longdouble(360.0 / 6.387230) / swiftest.JD2S
     }
 
     planetIpz = { # Only the polar moments of inertia are used currently. Where the quantity is unkown, we just use the value of a sphere = 0.4
@@ -262,7 +262,7 @@ def vec2xr(param: Dict, **kwargs: Any):
     rhill : float or array-like of float, optional
         Hill's radius values if these are massive bodies
     rot:  (n,3) array-like of float, optional
-        Rotation rate vectors if these are massive bodies with rotation enabled. This can be used instead of passing
+        Rotation rate vectors if these are massive bodies with rotation enabled in deg/TU
     Ip: (n,3) array-like of flaot, optional
         Principal axes moments of inertia vectors if these are massive bodies with rotation enabled. This can be used
         instead of passing Ip1, Ip2, and Ip3 separately

src/collision/collision_io.f90

@@ -410,7 +410,7 @@ module subroutine collision_io_netcdf_write_frame_snapshot(self, history, param)
                      call netcdf_io_check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[   idslot, stage, eslot]), "collision_io_netcdf_write_frame_snapshot nf90_put_var radius_varid"  )
                      if (param%lrotation) then
                         call netcdf_io_check( nf90_put_var(nc%id, nc%Ip_varid,     pl%Ip(:,i),   start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var Ip_varid"  )
-                        call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid,    pl%rot(:,i),  start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var rotx_varid"  )
+                        call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid,    pl%rot(:,i)*RAD2DEG,  start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var rotx_varid"  )
                      end if
                   end do
                end do

src/encounter/encounter_io.f90

@@ -261,7 +261,7 @@ module subroutine encounter_io_netcdf_write_frame_snapshot(self, history, param)
 
                if (param%lrotation) then
                   call netcdf_io_check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, tslot], count=[NDIM,1,1]), "encounter_io_netcdf_write_frame_snapshot nf90_put_var pl Ip_varid"  )
-                  call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i), start=[1,idslot, tslot], count=[NDIM,1,1]), "encounter_io_netcdf_write_frame_snapshot nf90_put_var pl rotx_varid"  )
+                  call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i)*RAD2DEG, start=[1,idslot, tslot], count=[NDIM,1,1]), "encounter_io_netcdf_write_frame_snapshot nf90_put_var pl rotx_varid"  )
                end if
 
                charstring = trim(adjustl(pl%info(i)%name))

src/swiftest/swiftest_io.f90

@@ -621,6 +621,7 @@ module subroutine swiftest_io_netcdf_get_t0_values_system(self, param)
 
             if (param%lrotation) then
                call netcdf_io_check( nf90_get_var(nc%id, nc%rot_varid, rot0, start=[1,1,tslot], count=[NDIM,1,1]), "netcdf_io_get_t0_values_system rot_varid" )
+               rot0(:) = rot0(:) * DEG2RAD
                call netcdf_io_check( nf90_get_var(nc%id, nc%Ip_varid, Ip0, start=[1,1,tslot], count=[NDIM,1,1]), "netcdf_io_get_t0_values_system Ip_varid" )
                cb%L0(:) = Ip0(3) * mass0 * cb%R0**2 * rot0(:)
             end if
@@ -1174,7 +1175,8 @@ module function swiftest_io_netcdf_read_frame_system(self, nc, param) result(ier
             end do
 
             call netcdf_io_check( nf90_get_var(nc%id, nc%rot_varid, vectemp, start=[1, 1, tslot], count=[NDIM,idmax,1]), "netcdf_io_read_frame_system nf90_getvar rot_varid"  )
-            cb%rot(:) = vectemp(:,1)
+            vectemp(:,:) = vectemp(:,:) * DEG2RAD
+            cb%rot(:) = vectemp(:,1) 
             do i = 1, NDIM
                if (npl > 0) pl%rot(i,:) = pack(vectemp(i,:), plmask(:))
             end do
@@ -1608,7 +1610,7 @@ module subroutine swiftest_io_netcdf_write_frame_body(self, nc, param)
                   if (param%lclose) call netcdf_io_check( nf90_put_var(nc%id, nc%radius_varid, self%radius(j), start=[idslot, tslot]), "netcdf_io_write_frame_body nf90_put_var body radius_varid"  )
                   if (param%lrotation) then
                      call netcdf_io_check( nf90_put_var(nc%id, nc%Ip_varid, self%Ip(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "netcdf_io_write_frame_body nf90_put_var body Ip_varid"  )
-                     call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid, self%rot(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "netcdf_io_write_frame_body nf90_put_var body rotx_varid"  )
+                     call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid, self%rot(:, j) * RAD2DEG, start=[1,idslot, tslot], count=[NDIM,1,1]), "netcdf_io_write_frame_body nf90_put_var body rotx_varid"  )
                   end if
                   ! if (param%ltides) then
                   !    call netcdf_io_check( nf90_put_var(nc%id, nc%k2_varid, self%k2(j), start=[idslot, tslot]), "netcdf_io_write_frame_body nf90_put_var body k2_varid"  )
@@ -1653,7 +1655,7 @@ module subroutine swiftest_io_netcdf_write_frame_cb(self, nc, param)
          call netcdf_io_check( nf90_put_var(nc%id, nc%j4rp4_varid, self%j4rp4, start=[tslot]), "netcdf_io_write_frame_cb nf90_put_var cb j4rp4_varid" )
          if (param%lrotation) then
             call netcdf_io_check( nf90_put_var(nc%id, nc%Ip_varid, self%Ip(:), start=[1, idslot, tslot], count=[NDIM,1,1]), "netcdf_io_write_frame_cb nf90_put_var cb Ip_varid"  )
-            call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid, self%rot(:), start=[1, idslot, tslot], count=[NDIM,1,1]), "netcdf_io_write_frame_cby nf90_put_var cb rot_varid"  )
+            call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid, self%rot(:) * RAD2DEG, start=[1, idslot, tslot], count=[NDIM,1,1]), "netcdf_io_write_frame_cby nf90_put_var cb rot_varid"  )
          end if
 
          call netcdf_io_check( nf90_set_fill(nc%id, old_mode, old_mode), "netcdf_io_write_frame_cb nf90_set_fill old_mode"  )
@@ -2687,6 +2689,7 @@ module subroutine swiftest_io_read_in_cb(self, param)
       if (param%lrotation) then
          read(iu, *, err = 667, iomsg = errmsg) self%Ip(1), self%Ip(2), self%Ip(3)
          read(iu, *, err = 667, iomsg = errmsg) self%rot(1), self%rot(2), self%rot(3)
+         self%rot(:) = self%rot(:) * DEG2RAD
       end if
       ierr = 0
       close(iu, err = 667, iomsg = errmsg)
@@ -2832,6 +2835,7 @@ module function swiftest_io_read_frame_body(self, iu, param) result(ierr)
                   if (param%lrotation) then
                      read(iu, *, err = 667, iomsg = errmsg) self%Ip(1, i), self%Ip(2, i), self%Ip(3, i)
                      read(iu, *, err = 667, iomsg = errmsg) self%rot(1, i), self%rot(2, i), self%rot(3, i)
+                     self%rot(:,i) = self%rot(:,i) * DEG2RAD 
                   end if
                   ! if (param%ltides) then
                   !    read(iu, *, err = 667, iomsg = errmsg) self%k2(i)