Fixed angular momentum jumps due to escaped bodies. Now accounting fo…

…r the shift in the system barycenter due to the loss of bodies to the outer limit

examples/symba_energy_momentum/escape.in

@@ -3,15 +3,15 @@
 0.0 0.0 0.0       ! x y z
 0.0 0.0 0.0       ! vx vy vz
 0.0 0.0 0.07      ! ip
-11.2093063  -38.75937204  82.25088158     ! rot (radian / year)
+0.0 0.0 0.0 !11.2093063  -38.75937204  82.25088158     ! rot (radian / year)
 2 1e-07 0.0009
-7e-06
+7e-05
 99.9 0.0 0.0
 100.00 10.00 0.0
 0.4 0.4 0.4 !Ip 
-0.0 0.0 0.0 !rot
+0.0 0.0 1000.0 !rot
 3 1e-08 0.0004
-3.25e-06
+3.25e-05
 1.0 4.20E-05 0.0
 0.00 -6.28 0.0
 0.4 0.4 0.4 !Ip 

examples/symba_energy_momentum/param.escape.in

@@ -1,6 +1,6 @@
 T0             0.0e0
-TSTOP          0.00100           ! simulation length in seconds = 100 years
-DT             0.00010              ! stepsize in seconds
+TSTOP          1e2           ! simulation length in seconds = 100 years
+DT             1.00              ! stepsize in seconds
 PL_IN          escape.in
 TP_IN          tp.in
 IN_TYPE        ASCII
@@ -14,16 +14,16 @@ ISTEP_DUMP     1                 ! system dump cadence
 J2             0.0               ! no J2 term
 J4             0.0               ! no J4 term
 CHK_CLOSE      yes               ! check for planetary close encounters
-CHK_RMIN       0.005
-CHK_RMAX       1e2
+CHK_RMIN       0.00465047       ! check for close solar encounters in AU
+CHK_RMAX       10000.0          ! discard outside of
 CHK_EJECT      -1.0              ! ignore this check
 CHK_QMIN       -1.0              ! ignore this check
 !CHK_QMIN_COORD HELIO            ! commented out here
 !CHK_QMIN_RANGE 1.0 1000.0       ! commented out here
 ENC_OUT        enc.escape.dat
 EXTRA_FORCE    no                ! no extra user-defined forces
 BIG_DISCARD    no                ! output all planets if anything discarded
-RHILL_PRESENT  yes               ! Hill's sphere radii in input file
+RHILL_PRESENT  no                ! Hill's sphere radii in input file
 MTINY          1.0e-16
 FRAGMENTATION  yes
 MU2KG          1.98908e30

examples/symba_energy_momentum/param.sun.in

@@ -1,6 +1,6 @@
 !Parameter file for the SyMBA-RINGMOONS test
 T0             0.0 
-TSTOP          1.0
+TSTOP          3.0e-2
 DT             1e-3
 PL_IN          sun.in
 TP_IN          tp.in

examples/symba_energy_momentum/sun.in

@@ -5,11 +5,11 @@
 0.0 0.0 0.07      ! ip
 11.2093063  -38.75937204  82.25088158     ! rot (radian / year)
 2 2e-08 
-3e-06
+3e-04
 5e-2 0.0 0.0
 0.00 10.00 0.0
 0.4 0.4 0.4 !Ip 
-0.0 0.0 2300.0 !rot
+100.0 100000.0 -2300.0 !rot
 3 2e-08 
 3e-06
 1.0 0.00E-05 0.0

src/coord/coord_h2b.f90

@@ -11,12 +11,13 @@ subroutine coord_h2b(npl, swiftest_plA, msys)
 
 ! arguments
    integer(I4B), intent(in)  :: npl
-   real(DP), intent(out)   :: msys
+   real(DP), intent(out), optional   :: msys
    type(swiftest_pl),intent(inout) :: swiftest_plA
 
 ! internals
    integer(I4B)          :: i
    logical, dimension(npl) :: lstatus
+   real(DP)              :: mtot
 
 ! executable code
    associate(vbcb => swiftest_plA%vb(:,1), xbcb => swiftest_plA%xb(:,1), &
@@ -35,10 +36,11 @@ subroutine coord_h2b(npl, swiftest_plA, msys)
          vbcb(:) = vbcb(:) + mass(i) * vh(:,i)
       end do
 
-      msys = dMcb + sum(mass(2:npl), lstatus(2:npl)) + Mcb_initial
+      mtot = dMcb + sum(mass(2:npl), lstatus(2:npl)) + Mcb_initial
+      if (present(msys)) msys = mtot
 
-      xbcb(:) = -xbcb(:) / msys                      
-      vbcb(:) = -vbcb(:) / msys    
+      xbcb(:) = -xbcb(:) / mtot
+      vbcb(:) = -vbcb(:) / mtot
 
       do i = 2,npl 
          xb(:,i) = xh(:,i) + xbcb(:)

src/io/io_conservation_report.f90

@@ -22,10 +22,10 @@ module subroutine io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param,
       real(DP)                        :: Eorbit_error, Etotal_error, Ecoll_error
       real(DP)                        :: Mtot_now, Merror
       real(DP)                        :: Lmag_now, Lerror
-      character(len=*), parameter     :: egyfmt = '(ES23.16,10(",",ES23.16,:))' ! Format code for all simulation output
-      character(len=*), parameter     :: egyheader = '("t,Eorbit,Ecollisions,Lx,Ly,Lz,Mtot")'
-      integer(I4B), parameter         :: egyiu = 72
-      character(len=*), parameter     :: egytermfmt = '("  DL/L0 = ", ES12.5 &
+      character(len=*), parameter     :: EGYFMT = '(ES23.16,10(",",ES23.16,:))' ! Format code for all simulation output
+      character(len=*), parameter     :: EGYHEADER = '("t,Eorbit,Ecollisions,Lx,Ly,Lz,Mtot")'
+      integer(I4B), parameter         :: EGYIU = 72
+      character(len=*), parameter     :: EGYTERMFMT = '("  DL/L0 = ", ES12.5 &
                                                          "; DEcollisions/|E0| = ", ES12.5, &
                                                          "; D(Eorbit+Ecollisions)/|E0| = ", ES12.5, &
                                                          "; DM/M0 = ", ES12.5)'
@@ -37,10 +37,10 @@ module subroutine io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param,
                  lfirst => param%lfirstenergy)
          if (lfirst) then
             if (param%out_stat == "OLD") then
-               open(unit = egyiu, file = energy_file, form = "formatted", status = "old", action = "write", position = "append")
+               open(unit = EGYIU, file = ENERGY_FILE, form = "formatted", status = "old", action = "write", position = "append")
             else 
-               open(unit = egyiu, file = energy_file, form = "formatted", status = "replace", action = "write")
-               write(egyiu,egyheader)
+               open(unit = EGYIU, file = ENERGY_FILE, form = "formatted", status = "replace", action = "write")
+               write(EGYIU,EGYHEADER)
             end if
          end if
          call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_now, ke_spin_now, pe_now, Lorbit_now, Lspin_now)
@@ -57,8 +57,8 @@ module subroutine io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param,
             lfirst = .false.
          end if
 
-         write(egyiu,egyfmt) t, Eorbit_now, Ecollisions, Ltot_now, Mtot_now
-         flush(egyiu)
+         write(EGYIU,EGYFMT) t, Eorbit_now, Ecollisions, Ltot_now, Mtot_now
+         flush(EGYIU)
          if (.not.lfirst .and. lterminal) then 
             Lmag_now = norm2(Ltot_now)
             Lerror = norm2(Ltot_now - Ltot_orig) / Lmag_orig

src/modules/module_interfaces.f90

@@ -58,7 +58,7 @@ SUBROUTINE coord_h2b(npl, swiftest_plA, msys)
                USE swiftest_data_structures
                IMPLICIT NONE
                INTEGER(I4B), INTENT(IN)         :: npl
-               REAL(DP), INTENT(OUT)            :: msys
+               REAL(DP), INTENT(OUT), optional            :: msys
                TYPE(swiftest_pl), INTENT(INOUT) :: swiftest_plA
           END SUBROUTINE coord_h2b
      END INTERFACE

src/symba/symba_discard_conserve_mtm.f90

@@ -12,17 +12,18 @@ subroutine symba_discard_conserve_mtm(param, swiftest_plA, ipl, lescape_body)
    integer(I4B), intent(in)    :: ipl
    logical, intent(in)         :: lescape_body
    ! Internals
-   real(DP), dimension(NDIM) :: Lpl, Lcb, xcom, vcom
+   real(DP), dimension(NDIM) :: Lpl, Ltot, Lcb, xcom, vcom
    real(DP)                  :: pe, ke_orbit, ke_spin
-   integer(I4B)              :: i
+   integer(I4B)              :: i, oldstat
 
 
    associate(npl => swiftest_plA%nbody, xb => swiftest_plA%xb, vb => swiftest_plA%vb, &
       rot => swiftest_plA%rot, Ip => swiftest_plA%Ip, radius => swiftest_plA%radius, mass => swiftest_plA%mass, &
       Lcb_initial => swiftest_plA%Lcb_initial, dLcb => swiftest_plA%dLcb, Ecollisions => swiftest_plA%Ecollisions, &
       Rcb_initial => swiftest_plA%Rcb_initial, dRcb => swiftest_plA%dRcb, &
       Mcb_initial => swiftest_plA%Mcb_initial, dMcb => swiftest_plA%dMcb, &
-      Mescape => swiftest_plA%Mescape, Lescape => swiftest_plA%Lescape, Euntracked => swiftest_plA%Euntracked)
+      Mescape => swiftest_plA%Mescape, Lescape => swiftest_plA%Lescape, Euntracked => swiftest_plA%Euntracked, &
+      status => swiftest_plA%status)
 
       ! Add the potential and kinetic energy of the lost body to the records
       pe = -mass(1) * mass(ipl) / norm2(xb(:, ipl) - xb(:, 1))
@@ -33,23 +34,38 @@ subroutine symba_discard_conserve_mtm(param, swiftest_plA, ipl, lescape_body)
       ! Add planet mass to central body accumulator
       if (lescape_body) then
          Mescape = Mescape + mass(ipl)
-         call util_crossproduct(xb(:,ipl), vb(:,ipl), Lpl)
-         Lpl(:) = mass(ipl) * (Lpl(:) + radius(ipl)**2 * Ip(3,ipl) * rot(:, ipl))
-         Lescape(:) = Lescape(:) + Lpl(:) 
          do i = 2, npl
             if (i == ipl) cycle
             pe = pe - mass(i) * mass(ipl) / norm2(xb(:, ipl) - xb(:, i))
          end do
+
+         Ltot(:) = 0.0_DP
+         do i = 1, npl
+            call util_crossproduct(mass(i) * xb(:,i), vb(:,i), Lpl)
+            Ltot(:) = Ltot(:) + Lpl(:)
+         end do
+         call coord_b2h(npl, swiftest_plA)
+         oldstat = status(ipl)
+         status(ipl) = INACTIVE
+         call coord_h2b(npl, swiftest_plA)
+         status(ipl) = oldstat
+         do i = 1, npl
+            if (i == ipl) cycle
+            call util_crossproduct(mass(i) * xb(:,i), vb(:,i), Lpl)
+            Ltot(:) = Ltot(:) - Lpl(:) 
+         end do 
+         Lescape(:) = Lescape(:) + Ltot(:) + mass(ipl) * radius(ipl)**2 * Ip(3, ipl) * rot(:, ipl)
+
       else
          xcom(:) = (mass(ipl) * xb(:, ipl) + mass(1) * xb(:,1)) / (mass(1) + mass(ipl))
          vcom(:) = (mass(ipl) * vb(:, ipl) + mass(1) * vb(:,1)) / (mass(1) + mass(ipl))
-
          call util_crossproduct(xb(:,ipl) - xcom(:), vb(:,ipl) - vcom(:), Lpl)
          Lpl(:) = mass(ipl) * (Lpl(:) + radius(ipl)**2 * Ip(3,ipl) * rot(:, ipl))
 
          call util_crossproduct(xb(:, 1) - xcom(:), vb(:, 1) - vcom(:), Lcb)
          Lcb(:) = mass(1) * Lcb(:) 
-         ke_orbit = 0.5_DP * mass(1) * dot_product(vb(:, 1), vb(:, 1)) + ke_orbit
+
+         ke_orbit = ke_orbit + 0.5_DP * mass(1) * dot_product(vb(:, 1), vb(:, 1)) 
          ke_spin = ke_spin + 0.5_DP * mass(1) * radius(1)**2 * Ip(3, 1) * dot_product(rot(:, 1), rot(:, 1))
          ! Update mass of central body to be consistent with its total mass
          dMcb = dMcb + mass(ipl)
@@ -66,7 +82,7 @@ subroutine symba_discard_conserve_mtm(param, swiftest_plA, ipl, lescape_body)
          vb(:, 1) = vcom(:)
          ke_orbit = ke_orbit - 0.5_DP * mass(1) * dot_product(vb(:, 1), vb(:, 1)) 
       end if
-      ! Update position and velocity of central body
+      call coord_b2h(npl, swiftest_plA)
 
       ! We must do this for proper book-keeping, since we can no longer track this body's contribution to energy directly
       if (lescape_body) then