Fixed indexing issue in the shgrav acceleration calculation and fixed…

… some issues when recording discards

examples/spherical_harmonics_cb/spherical_harmonics_cb.py

@@ -82,11 +82,10 @@
 capm_tp     = rng.uniform(0.0, 360.0, ntp)
 
 sim.add_body(name=name_tp, a=a_tp, e=e_tp, inc=inc_tp, capom=capom_tp, omega=omega_tp, capm=capm_tp)
-sim.set_parameter(tstart=0.0, tstop=5.0e2, dt=0.01, istep_out=100, dump_cadence=0, compute_conservation_values=True, mtiny=mtiny)
+sim.set_parameter(tstart=0.0, tstop=1.0, dt=0.01, istep_out=1, dump_cadence=0, compute_conservation_values=True, mtiny=mtiny)
 
 # Display the run configuration parameters.
 sim.get_parameter()
-sim.save()
 
 # Run the simulation. Arguments may be defined here or thorugh the swiftest.Simulation() method.
 sim.run()

src/shgrav/shgrav_accel.f90

@@ -158,14 +158,14 @@ module subroutine shgrav_acc(body, nbody_system)
             class is (swiftest_pl)
                 do i = 1, body%nbody
                     if (body%lmask(i)) then
-                        call shgrav_g_acc_one(cb%Gmass, cb%radius, cb%rotphase, body%rh(:,i), cb%c_lm, body%aobl, &
+                        call shgrav_g_acc_one(cb%Gmass, cb%radius, cb%rotphase, body%rh(:,i), cb%c_lm, body%aobl(:,i), &
                             GMpl=body%Gmass(i), aoblcb=cb%aobl)
                     end if
                 end do
             class is (swiftest_tp)
                 do i = 1, body%nbody
                     if (body%lmask(i)) then
-                        call shgrav_g_acc_one(cb%Gmass, cb%radius, cb%rotphase, body%rh(:,i), cb%c_lm, body%aobl)
+                        call shgrav_g_acc_one(cb%Gmass, cb%radius, cb%rotphase, body%rh(:,i), cb%c_lm, body%aobl(:,i))
                     end if
                 end do
             end select

src/swiftest/swiftest_discard.f90

@@ -49,8 +49,18 @@ module subroutine swiftest_discard_system(self, param)
          end if
 
          if (ltp_discards.or.lpl_discards) then
+            ! Advance the collision id number and save it
+            collider%maxid_collision = collider%maxid_collision + 1
+            collider%collision_id = collider%maxid_collision
+            collider%impactors%regime = COLLRESOLVE_REGIME_MERGE
+            write(idstr,*) collider%collision_id
+            call swiftest_io_log_one_message(COLLISION_LOG_OUT, "collision_id " // trim(adjustl(idstr)))
+
             if (ltp_discards) then
                allocate(ldiscard, source=tp%ldiscard(:))
+               do i = 1, ntp
+                  if (ldiscard(i)) call tp%info(i)%set_value(collision_id=collider%collision_id)
+               end do
                allocate(tpsub, mold=tp)
                call tp%spill(tpsub, ldiscard, ldestructive=.true.)
                nsub = tpsub%nbody
@@ -70,6 +80,9 @@ module subroutine swiftest_discard_system(self, param)
                                               ! simply used to trigger a snapshot.
                if (param%lenergy) call self%conservation_report(param, lterminal=.false.)
                allocate(ldiscard, source=pl%ldiscard(:))
+               do i = 1, npl
+                  if (ldiscard(i)) call pl%info(i)%set_value(collision_id=collider%collision_id)
+               end do
                allocate(plsub, mold=pl)
                call pl%spill(plsub, ldiscard, ldestructive=.false.)
                nsub = plsub%nbody
@@ -86,13 +99,7 @@ module subroutine swiftest_discard_system(self, param)
                end select
                call pl_discards%setup(0,param) 
             end if
-            ! Advance the collision id number and save it
-            collider%maxid_collision = max(collider%maxid_collision, maxval(nbody_system%pl%info(:)%collision_id))
-            collider%maxid_collision = collider%maxid_collision + 1
-            collider%collision_id = collider%maxid_collision
-            collider%impactors%regime = COLLRESOLVE_REGIME_MERGE
-            write(idstr,*) collider%collision_id
-            call swiftest_io_log_one_message(COLLISION_LOG_OUT, "collision_id " // trim(adjustl(idstr)))
+
 
             call collision_history%take_snapshot(param,nbody_system, t, "particle") 
          end if

src/swiftest/swiftest_obl.f90

@@ -246,7 +246,6 @@ module subroutine swiftest_non_spherical_cb_acc_pl(self, nbody_system)
       end associate
 
       return
-
    end subroutine swiftest_non_spherical_cb_acc_pl
 
 
@@ -289,8 +288,8 @@ module subroutine swiftest_non_spherical_cb_acc_tp(self, nbody_system)
          end do
 
       end associate
-      return
 
+      return
    end subroutine swiftest_non_spherical_cb_acc_tp