Got close encounter saves in place. Now just troubleshooting

examples/Fragmentation/Fragmentation_Movie.py

@@ -41,7 +41,7 @@
 movie_titles = dict(zip(available_movie_styles, movie_title_list))
 
 # These initial conditions were generated by trial and error
-pos_vectors = {"disruption_headon"         : [np.array([1.0, -5.0e-05, 0.0]),
+pos_vectors = {"disruption_headon"         : [np.array([1.0005, -5.0e-05, 0.0]),
                                               np.array([1.0,  5.0e-05 ,0.0])],
                "supercatastrophic_off_axis": [np.array([1.0, -4.2e-05,      0.0]),
                                               np.array([1.0,  4.2e-05,      0.0])],
@@ -202,7 +202,7 @@ def data_stream(self, frame=0):
         # Set fragmentation parameters
         minimum_fragment_gmass = 0.2 * body_Gmass[style][1] # Make the minimum fragment mass a fraction of the smallest body
         gmtiny = 0.99 * body_Gmass[style][1] # Make GMTINY just smaller than the smallest original body. This will prevent runaway collisional cascades
-        sim.set_parameter(fragmentation=True, encounter_save="TRAJECTORY", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False)
+        sim.set_parameter(fragmentation=True, encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False)
         sim.run(dt=1e-4, tstop=1.0e-3, istep_out=1, dump_cadence=1)
 
         print("Generating animation")

python/swiftest/swiftest/simulation_class.py

@@ -214,10 +214,12 @@ def __init__(self,read_param: bool = False, read_old_output: bool = False, simdi
             Check for close encounters between bodies. If set to True, then the radii of massive bodies must be included
             in initial conditions.
             Parameter input file equivalent: `CHK_CLOSE`
-        encounter_save : {"NONE","TRAJECTORY","CLOSEST"}, default "NONE"
-            Indicate if and how encounter data should be saved. If set to "TRAJECTORY" the full close encounter
-            trajectories are saved to file. If set to "CLOSEST" only the trajectories at the time of closest approach
-            are saved. If set to "NONE" no trajectory information is saved.
+        encounter_save : {"NONE","TRAJECTORY","CLOSEST", "BOTH"}, default "NONE"
+            Indicate if and how encounter data should be saved. If set to "TRAJECTORY", the position and velocity vectors
+            of all bodies undergoing close encounters are saved at each intermediate step to the encounter files.
+            If set to "CLOSEST", the position  and velocities at the point of closest approach between pairs of bodies are 
+            computed and stored to the encounter files. If set to "BOTH", then this stores the values that would be computed
+            in "TRAJECTORY" and "CLOSEST". If set to "NONE" no trajectory information is saved.
             *WARNING*: Enabling this feature could lead to very large files.
         general_relativity : bool, default True
             Include the post-Newtonian correction in acceleration calculations.
@@ -1023,7 +1025,7 @@ def set_feature(self,
                     tides: bool | None = None,
                     interaction_loops: Literal["TRIANGULAR", "FLAT", "ADAPTIVE"] | None = None,
                     encounter_check_loops: Literal["TRIANGULAR", "SORTSWEEP", "ADAPTIVE"] | None = None,
-                    encounter_save: Literal["NONE", "TRAJECTORY", "CLOSEST"] | None = None,
+                    encounter_save: Literal["NONE", "TRAJECTORY", "CLOSEST", "BOTH"] | None = None,
                     verbose: bool | None = None,
                     **kwargs: Any
                     ):
@@ -1035,10 +1037,12 @@ def set_feature(self,
         close_encounter_check : bool, optional
             Check for close encounters between bodies. If set to True, then the radii of massive bodies must be included
             in initial conditions.
-        encounter_save : {"NONE","TRAJECTORY","CLOSEST"}, default "NONE"
-            Indicate if and how encounter data should be saved. If set to "TRAJECTORY" the full close encounter
-            trajectories are saved to file. If set to "CLOSEST" only the trajectories at the time of closest approach
-            are saved. If set to "NONE" no trajectory information is saved.
+        encounter_save : {"NONE","TRAJECTORY","CLOSEST","BOTH"}, default "NONE"
+            Indicate if and how encounter data should be saved. If set to "TRAJECTORY", the position and velocity vectors
+            of all bodies undergoing close encounters are saved at each intermediate step to the encounter files.
+            If set to "CLOSEST", the position  and velocities at the point of closest approach between pairs of bodies are 
+            computed and stored to the encounter files. If set to "BOTH", then this stores the values that would be computed
+            in "TRAJECTORY" and "CLOSEST". If set to "NONE" no trajectory information is saved.
             *WARNING*: Enabling this feature could lead to very large files.
         general_relativity : bool, optional
             Include the post-Newtonian correction in acceleration calculations.
@@ -1199,7 +1203,7 @@ def set_feature(self,
                 update_list.append("encounter_check_loops")
 
         if encounter_save is not None:
-            valid_vals = ["NONE", "TRAJECTORY", "CLOSEST"]
+            valid_vals = ["NONE", "TRAJECTORY", "CLOSEST", "BOTH"]
             encounter_save = encounter_save.upper()
             if encounter_save not in valid_vals:
                 msg = f"{encounter_save} is not a valid option for encounter_save."

src/encounter/encounter_util.f90

@@ -570,9 +570,9 @@ module subroutine encounter_util_snapshot_encounter(self, param, system, t, arg)
       character(*),                 intent(in), optional :: arg    !! Optional argument (needed for extended storage type used in collision snapshots)
       ! Arguments
       class(encounter_snapshot), allocatable :: snapshot
-      integer(I4B) :: i, j, k, npl_snap, ntp_snap
+      integer(I4B) :: i, j, k, npl_snap, ntp_snap, iflag
       real(DP), dimension(NDIM) :: rrel, vrel, rcom, vcom
-      real(DP) :: mu, a, q, capm, tperi
+      real(DP) :: Gmtot, a, q, capm, tperi
       real(DP), dimension(NDIM,2) :: rb,vb
 
       if (.not.present(t)) then
@@ -621,6 +621,8 @@ module subroutine encounter_util_snapshot_encounter(self, param, system, t, arg)
                                  ntp_snap = count(tp%lmask(1:ntp))
                               end if
 
+                              if (npl_snap + ntp_snap == 0) return ! Nothing to snapshot
+
                               pl_snap%nbody = npl_snap
 
                               ! Take snapshot of the currently encountering massive bodies
@@ -674,7 +676,7 @@ module subroutine encounter_util_snapshot_encounter(self, param, system, t, arg)
                                           pl_snap%id(:) = pl%id([i,j])
                                           pl_snap%info(:) = pl%info([i,j])
                                           pl_snap%Gmass(:) = pl%Gmass([i,j])
-                                          mu = sum(pl_snap%Gmass(:))
+                                          Gmtot = sum(pl_snap%Gmass(:))
                                           if (param%lclose) pl_snap%radius(:) = pl%radius([i,j])
                                           if (param%lrotation) then
                                              do i = 1, NDIM
@@ -686,15 +688,36 @@ module subroutine encounter_util_snapshot_encounter(self, param, system, t, arg)
                                           ! Compute pericenter passage time to get the closest approach parameters
                                           rrel(:) = plplenc_list%r2(:,k) - plplenc_list%r1(:,k)
                                           vrel(:) = plplenc_list%v2(:,k) - plplenc_list%v1(:,k)
-                                          call orbel_xv2aqt(mu, rrel(1), rrel(2), rrel(3), vrel(1), vrel(2), vrel(3), a, q, capm, tperi)
+                                          call orbel_xv2aqt(Gmtot, rrel(1), rrel(2), rrel(3), vrel(1), vrel(2), vrel(3), a, q, capm, tperi)
                                           snapshot%t = t + tperi
 
                                           ! Computer the center mass of the pair
-
-                                          ! do i = 1, NDIM
-                                          !    pl_snap%rh(i,:) = pl%rh(i,[i,j])
-                                          !    pl_snap%vh(i,:) = pl%vb(i,1:npl), pl%lmask(1:npl))
-                                          ! end do
+                                          rcom(:) = (plplenc_list%r1(:,k) * pl_snap%Gmass(1) + plplenc_list%r2(:,k) * pl_snap%Gmass(2)) / Gmtot
+                                          vcom(:) = (plplenc_list%v1(:,k) * pl_snap%Gmass(1) + plplenc_list%v2(:,k) * pl_snap%Gmass(2)) / Gmtot
+                                          rb(:,1) = plplenc_list%r1(:,k) - rcom(:)
+                                          rb(:,2) = plplenc_list%r2(:,k) - rcom(:)
+                                          vb(:,1) = plplenc_list%v1(:,k) - vcom(:)
+                                          vb(:,2) = plplenc_list%v2(:,k) - vcom(:)
+
+                                          ! Drift the relative orbit to get the new relative position and velocity
+                                          call drift_one(Gmtot, rrel(1), rrel(2), rrel(3), vrel(1), vrel(2), vrel(3), tperi, iflag)
+                                          if (iflag /= 0) write(*,*) "Danby error in encounter_util_snapshot_encounter. Closest approach positions and vectors may not be accurate."
+
+                                          ! Get the new position and velocity vectors
+                                          rb(:,1) = -(pl_snap%Gmass(2) / Gmtot) * rrel(:)
+                                          rb(:,2) =  (pl_snap%Gmass(1)) / Gmtot * rrel(:)
+
+                                          vb(:,1) = -(pl_snap%Gmass(2) / Gmtot) * vrel(:)
+                                          vb(:,2) =  (pl_snap%Gmass(1)) / Gmtot * vrel(:)
+
+                                          ! Move the CoM assuming constant velocity over the time it takes to reach periapsis
+                                          !rcom(:) = rcom(:) + vcom(:) * tperi
+
+                                          ! Compute the heliocentric position and velocity vector at periapsis
+                                          pl_snap%rh(:,1) = rb(:,1) + rcom(:)
+                                          pl_snap%rh(:,2) = rb(:,2) + rcom(:)
+                                          pl_snap%vh(:,1) = vb(:,1) + vcom(:)
+                                          pl_snap%vh(:,2) = vb(:,2) + vcom(:)
 
                                           call pl_snap%sort("id", ascending=.true.)
                                           call encounter_util_save_encounter(param%encounter_history,snapshot,snapshot%t)

src/symba/symba_io.f90

@@ -118,15 +118,18 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms
          ! All reporting of collision information in SyMBA (including mergers) is now recorded in the Fraggle logfile
          call io_log_start(param, FRAGGLE_LOG_OUT, "Fraggle logfile")
 
-         if ((param%encounter_save /= "NONE") .and. (param%encounter_save /= "TRAJECTORY") .and. (param%encounter_save /= "CLOSEST")) then
+         if ((param%encounter_save /= "NONE")       .and. &
+             (param%encounter_save /= "TRAJECTORY") .and. &
+             (param%encounter_save /= "CLOSEST")    .and. &
+             (param%encounter_save /= "BOTH")) then
             write(iomsg,*) 'Invalid encounter_save parameter: ',trim(adjustl(param%out_type))
-            write(iomsg,*) 'Valid options are NONE, TRAJECTORY, or CLOSEST'
+            write(iomsg,*) 'Valid options are NONE, TRAJECTORY, CLOSEST, or BOTH'
             iostat = -1
             return
          end if
 
-         param%lenc_save_trajectory = (param%encounter_save == "TRAJECTORY") 
-         param%lenc_save_closest = (param%encounter_save == "CLOSEST") .or. param%lenc_save_trajectory ! Closest approaches are always saved when trajectories are saved
+         param%lenc_save_trajectory = (param%encounter_save == "TRAJECTORY") .or. (param%encounter_save == "BOTH")
+         param%lenc_save_closest = (param%encounter_save == "CLOSEST") .or. (param%encounter_save == "BOTH")
 
          ! Call the base method (which also prints the contents to screen)
          call io_param_reader(param, unit, iotype, v_list, iostat, iomsg)