Merge branch 'debug'

python/swiftest/swiftest/io.py

@@ -820,6 +820,7 @@ def process_netcdf_input(ds, param):
     elif param['OUT_TYPE'] == "NETCDF_FLOAT":
         ds = fix_types(ds,ftype=np.float32)
 
+    ds = ds.where(ds['id']>=0, drop=True)
     return ds
 
 def swiftest2xr(param, verbose=True):

src/collision/collision_regime.f90

@@ -65,7 +65,7 @@ subroutine collision_regime_LS12(impactors, nbody_system, param)
       integer(I4B) :: jtarg, jproj
       real(DP), dimension(2) :: radius_si, mass_si, density_si
       real(DP) :: min_mfrag_si, Mcb_si
-      real(DP), dimension(NDIM)  :: x1_si, v1_si, x2_si, v2_si, runit
+      real(DP), dimension(NDIM)  :: x1_si, v1_si, x2_si, v2_si
       real(DP) :: mlr, mslr, mtot, dentot
       integer(I4B), parameter :: NMASS_DIST = 3   !! Number of mass bins returned by the regime calculation (largest fragment, second largest, and remainder)  
 
@@ -149,7 +149,7 @@ subroutine collision_regime_LS12_SI(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2,
       real(DP), parameter   :: CRUFU = 2.0_DP - 3 * MU_BAR ! central potential variable from Rufu and Aharonson (2019)
       real(DP), parameter   :: SUPERCAT_QRATIO = 1.8_DP ! See Section 4.1 of LS12
       ! Internals
-      real(DP)           :: a1, alpha, aint, b, bcrit, c_star, egy, zeta, l, lint, mu, phi, theta
+      real(DP)           :: a1, alpha, aint, b, bcrit, c_star, egy, zeta, l, lint, mu, phi, theta, ke, pe, Qmerge
       real(DP)           :: Qr, Qrd_pstar, Qr_erosion, Qr_supercat
       real(DP)           :: Vhr, Verosion, Vescp, Vhill, Vimp, Vsupercat
       real(DP)           :: Mint, Mtot, Mtmp
@@ -213,32 +213,35 @@ subroutine collision_regime_LS12_SI(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2,
       Vhr = Vescp * (C1 * zeta**2 * theta**(2.5_DP) + C2 * zeta**2 + C3 * theta**(2.5_DP) + C4) ! Kokubo & Genda (2010) eq. (3)
       bcrit = rad1 / (rad1 + rad2)
       Qloss = 0.0_DP
-      U_binding = (3.0_DP * Mtot) / (5.0_DP * Rp) ! LS12 eq. 27
+      U_binding = (3 * GC * Mtot**2) / (5 * Rp) ! LS12 eq. 27
+      ke = 0.5_DP * Mtot * Vimp**2
+      pe = - GC * m1 * m2 / norm2(rh2 - rh1)
+      Qmerge = ke + pe + U_binding 
 
       if ((m1 < min_mfrag).or.(m2 < min_mfrag)) then 
          regime = COLLRESOLVE_REGIME_MERGE !perfect merging regime
          Mlr = Mtot
          Mslr = 0.0_DP
-         Qloss = 0.0_DP
+         Qloss = Qmerge
          call swiftest_io_log_one_message(COLLISION_LOG_OUT, &
                                  "Fragments would have mass below the minimum. Converting this collision into a merger.")
       else 
          if( Vimp < Vescp) then
             regime = COLLRESOLVE_REGIME_MERGE !perfect merging regime
             Mlr = Mtot
             Mslr = 0.0_DP
-            Qloss = 0.0_DP
+            Qloss = Qmerge
          else if (Vimp < Verosion) then 
             if (b < bcrit) then
                regime = COLLRESOLVE_REGIME_MERGE !partial accretion regime"
                Mlr = Mtot
                Mslr = 0.0_DP
-               Qloss = 0.0_DP
+               Qloss = Qmerge
             else if ((b > bcrit) .and. (Vimp < Vhr)) then
                regime = COLLRESOLVE_REGIME_MERGE ! graze and merge
                Mlr = Mtot
                Mslr = 0.0_DP
-               Qloss = 0.0_DP
+               Qloss = Qmerge
             else
                Mlr = m1
                Mslr = max(calc_Qrd_rev(m2, m1, Mint, den1, den2, Vimp, c_star), min_mfrag)
@@ -250,7 +253,7 @@ subroutine collision_regime_LS12_SI(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2,
                regime = COLLRESOLVE_REGIME_MERGE !cratering regime"
                Mlr = Mtot
                Mslr = 0.0_DP
-               Qloss = 0.0_DP
+               Qloss = Qmerge
             else 
                Mslr = max(Mtot * (3.0_DP - BETA) * (1.0_DP - N1 * Mlr / Mtot) / (N2 * BETA), min_mfrag)  ! LS12 eq (37)
                regime = COLLRESOLVE_REGIME_DISRUPTION !disruption

src/collision/collision_resolve.f90

@@ -325,7 +325,6 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status)
       integer(I4B) :: merge_text_length 
       character(len=NAMELEN) :: merge_text
       character(len=NAMELEN) :: newname, origin_type
-      character(len=STRMAX)  :: message
       real(DP) :: volume
 
       select type(nbody_system)

src/encounter/encounter_io.f90

@@ -232,7 +232,7 @@ module subroutine encounter_io_netcdf_write_frame_snapshot(self, history, param)
 
       ! Internals
       integer(I4B)           :: i, idslot, old_mode, npl, ntp
-      character(len=:), allocatable :: charstring
+      character(len=STRMAX) :: charstring
 
       select type(param)
       class is (swiftest_parameters)

src/fraggle/fraggle_generate.f90

@@ -420,7 +420,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
       real(DP), parameter     :: vmax_initial_factor = 5.0_DP ! For the initial "guess" of fragment velocities, this is the maximum velocity relative to escape velocity that the fragments will have
       integer(I4B), parameter :: MAXLOOP = 200
       integer(I4B), parameter :: MAXTRY  = 20
-      real(DP), parameter :: SUCCESS_METRIC = 1.0_DP
+      real(DP), parameter :: SUCCESS_METRIC = 1.0e-3_DP
       class(collision_fraggle), allocatable :: collider_local
       character(len=STRMAX) :: message
 

src/swiftest/swiftest_gr.f90

@@ -200,6 +200,10 @@ pure module subroutine swiftest_gr_vel2pseudovel(param, mu, rh, vh, pv)
          pv(1:NDIM) = vh(1:NDIM) ! Initial guess
          rterm = 3 * mu / norm2(rh(:))
          v2 = dot_product(vh(:), vh(:))
+         if (v2 < TINY(1.0_DP)) then
+            pv(:) = 0.0_DP
+            return
+         end if
          do n = 1, MAXITER
             pv2 = dot_product(pv(:), pv(:))
             G = 1.0_DP - inv_c2 * (0.5_DP * pv2 + rterm)

src/swiftest/swiftest_orbel.f90

@@ -939,6 +939,11 @@ pure module subroutine swiftest_orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, in
       capom = 0.0_DP
       omega = 0.0_DP
       capm = 0.0_DP
+      varpi = 0.0_DP
+      lam = 0.0_DP
+      f = 0.0_DP
+      cape = 0.0_DP
+      capf = 0.0_DP
       x = [px, py, pz]
       v = [vx, vy, vz]
       r = .mag. x(:)