From 56a5b3e08101b5f04d9c03d3f88f4979f5d0b144 Mon Sep 17 00:00:00 2001
From: David A Minton <daminton@purdue.edu>
Date: Wed, 8 Dec 2021 17:08:42 -0500
Subject: [PATCH] Added the missing gr_pos_kick step in the close encounter
 integrations. Moved the vb2vh call to be consistent with the regular helio
 step sequence

---
 examples/symba_energy_momentum/param.sun.in           | 11 ++++-------
 .../8pl_16tp_encounters/cb.in                         |  4 ++--
 .../8pl_16tp_encounters/cb.swiftest.in                |  4 ++--
 .../8pl_16tp_encounters/param.swiftest.in             |  5 +++--
 src/symba/symba_step.f90                              |  8 ++++++--
 5 files changed, 17 insertions(+), 15 deletions(-)

diff --git a/examples/symba_energy_momentum/param.sun.in b/examples/symba_energy_momentum/param.sun.in
index 351d8adb2..eaf0839f7 100644
--- a/examples/symba_energy_momentum/param.sun.in
+++ b/examples/symba_energy_momentum/param.sun.in
@@ -1,6 +1,6 @@
 !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
@@ -8,10 +8,10 @@ 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
@@ -19,8 +19,6 @@ 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
index 2df47f957..e3318aec0 100644
--- 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
index 2df47f957..e3318aec0 100644
--- 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
index bbf9e7089..7f4d21c9b 100644
--- 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
index eac112dcd..330f47a56 100644
--- 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