diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 0f6b8da0d..f0c25d196 100755 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/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__": diff --git a/python/swiftest/swiftest/init_cond.py b/python/swiftest/swiftest/init_cond.py index caeae710f..30fb46c6d 100644 --- a/python/swiftest/swiftest/init_cond.py +++ b/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 diff --git a/src/collision/collision_io.f90 b/src/collision/collision_io.f90 index bd556fd7a..d0e8f30ae 100644 --- a/src/collision/collision_io.f90 +++ b/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 diff --git a/src/encounter/encounter_io.f90 b/src/encounter/encounter_io.f90 index 32cabe8b5..ceb3194fe 100644 --- a/src/encounter/encounter_io.f90 +++ b/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)) diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index fc43abe54..a2b604357 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/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)