Added the missing gr_pos_kick step in the close encounter integration…

…s. Moved the vb2vh call to be consistent with the regular helio step sequence
MintonGroup · Dec 8, 2021 · 56a5b3e · 56a5b3e
1 parent 90bfc56
commit 56a5b3e
Show file tree

Hide file tree

Showing 5 changed files with 17 additions and 15 deletions.
diff --git a/examples/symba_energy_momentum/param.sun.in b/examples/symba_energy_momentum/param.sun.in
@@ -1,26 +1,24 @@
 !Parameter file for the SyMBA-RINGMOONS test
 T0             0.0 
-TSTOP          3.0e-2
+TSTOP          3.0e-1
 DT             1e-3
 CB_IN          cb.in
 PL_IN          sun.in
 TP_IN          tp.in
 IN_TYPE        ASCII
 ISTEP_OUT      1
 ISTEP_DUMP     1
-BIN_OUT        bin.sun.dat
+BIN_OUT        bin.sun.nc
 PARTICLE_OUT   particle.sun.dat
-OUT_TYPE       REAL8
-OUT_FORM       XV                ! osculating element output
+OUT_TYPE       NETCDF_DOUBLE
+OUT_FORM       XVEL              ! osculating element output
 OUT_STAT       REPLACE
 ISTEP_DUMP     1                 ! system dump cadence
 CHK_CLOSE      yes               ! check for planetary close encounters
 CHK_RMIN       0.005
 CHK_RMAX       1e2
 CHK_EJECT      -1.0              ! ignore this check
 CHK_QMIN       -1.0              ! ignore this check
-ENC_OUT        enc.sun.dat
-DISCARD_OUT    discard.sun.dat
 EXTRA_FORCE    no                ! no extra user-defined forces
 BIG_DISCARD    no                ! output all planets if anything discarded
 RHILL_PRESENT  no                ! Hill's sphere radii in input file
@@ -30,6 +28,5 @@ MU2KG          1.98908e30
 TU2S           3.1556925e7
 DU2M           1.49598e11
 ENERGY         yes
-ENERGY_OUT     energy.sun.out
 ROTATION       yes
 SEED           8 1230834 2346113 123409874 -123121105 -767545 -534058022 343309814 -12535638 
diff --git a/examples/symba_swifter_comparison/8pl_16tp_encounters/cb.in b/examples/symba_swifter_comparison/8pl_16tp_encounters/cb.in
@@ -1,5 +1,5 @@
 Sun
 0.00029591220819207774
 0.004650467260962157
-4.7535806948127355e-12
--2.2473967953572827e-18
+0.0
+0.0
diff --git a/examples/symba_swifter_comparison/8pl_16tp_encounters/cb.swiftest.in b/examples/symba_swifter_comparison/8pl_16tp_encounters/cb.swiftest.in
@@ -1,5 +1,5 @@
 Sun
 0.00029591220819207774
 0.004650467260962157
-4.7535806948127355e-12
--2.2473967953572827e-18
+0.0
+0.0
diff --git a/examples/symba_swifter_comparison/8pl_16tp_encounters/param.swiftest.in b/examples/symba_swifter_comparison/8pl_16tp_encounters/param.swiftest.in
@@ -32,6 +32,7 @@ FRAGMENTATION    NO
 ROTATION         NO
 TIDES            NO
 ENERGY           NO
-GR               NO
+GR               YES
 GMTINY           1e-12
-ENCOUNTER_CHECK  SORTSWEEP
+ENCOUNTER_CHECK  TRIANGULAR
+INTERACTION_LOOPS TRIANGULAR
diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90
@@ -71,22 +71,26 @@ module subroutine symba_step_interp_system(self, param, t, dt)
                   call pl%vh2vb(cb)
                   call pl%lindrift(cb, dth, lbeg=.true.)
                   call pl%kick(system, param, t, dth, lbeg=.true.)
+                  if (param%lgr) call pl%gr_pos_kick(param, dth)
                   call pl%drift(system, param, dt)
 
                   call tp%vh2vb(vbcb = -cb%ptbeg)
                   call tp%lindrift(cb, dth, lbeg=.true.)
                   call tp%kick(system, param, t, dth, lbeg=.true.)
+                  if (param%lgr) call tp%gr_pos_kick(param, dth)
                   call tp%drift(system, param, dt)
 
                   call system%recursive_step(param, t, 0)
 
                   call pl%kick(system, param, t, dth, lbeg=.false.)
-                  call pl%vb2vh(cb)
+                  if (param%lgr) call pl%gr_pos_kick(param, dth)
                   call pl%lindrift(cb, dth, lbeg=.false.)
+                  call pl%vb2vh(cb)
 
                   call tp%kick(system, param, t, dth, lbeg=.false.)
-                  call tp%vb2vh(vbcb = -cb%ptend)
+                  if (param%lgr) call tp%gr_pos_kick(param, dth)
                   call tp%lindrift(cb, dth, lbeg=.false.)
+                  call tp%vb2vh(vbcb = -cb%ptend)
                end select
             end select
          end select