Fixed bug that was preventing discards with the planet happening corr…

…ectly. Collisions due to small pericenter had to be passed back and forth between the heliocentric and planetocentric test particle structures.
MintonGroup · Jul 7, 2021 · bb05501 · bb05501
1 parent 62ba8c7
commit bb05501
Showing 1 changed file with 5 additions and 1 deletion.
diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90
@@ -509,13 +509,15 @@ subroutine rmvs_make_planetocentric(system, param)
                               tpenci%cb_heliocentric = cb
                               tpenci%ipleP = i
                               tpenci%status(:) = ACTIVE
-                              tpenci%name(:) = pack(tp%name(:), encmask(:)) 
                               ! Grab all the encountering test particles and convert them to a planetocentric frame
+                              tpenci%name(:) = pack(tp%name(:), encmask(:)) 
                               do j = 1, NDIM 
                                  tpenci%xheliocentric(j, :) = pack(tp%xh(j,:), encmask(:)) 
                                  tpenci%xh(j, :) = tpenci%xheliocentric(j, :) - pl%inner(0)%x(j, i)
                                  tpenci%vh(j, :) = pack(tp%vh(j,:), encmask(:)) - pl%inner(0)%v(j, i)
                               end do
+                              tpenci%lperi(:) = pack(tp%lperi(:), encmask(:)) 
+                              tpenci%plperP(:) = pack(tp%plperP(:), encmask(:)) 
                               ! Make sure that the test particles get the planetocentric value of mu 
                               allocate(cbenci%inner(0:NTPHENC))
                               do inner_index = 0, NTPHENC 
@@ -590,6 +592,8 @@ subroutine rmvs_end_planetocentric(system)
                                  tp%xh(j, tpind(1:pl%nenc(i))) = tpenci%xh(j,1:pl%nenc(i)) + pl%inner(NTPHENC)%x(j, i)
                                  tp%vh(j, tpind(1:pl%nenc(i))) = tpenci%vh(j,1:pl%nenc(i)) + pl%inner(NTPHENC)%v(j, i)
                               end do
+                              tp%lperi(tpind(1:pl%nenc(i))) = tpenci%lperi(tpind(1:pl%nenc(i)))
+                              tp%plperP(tpind(1:pl%nenc(i))) = tpenci%plperP(tpind(1:pl%nenc(i)))
                               deallocate(pl%planetocentric(i)%tp)
                               deallocate(cbenci%inner)
                               do inner_index = 0, NTPHENC