diff --git a/examples/symba_energy_momentum/escape.in b/examples/symba_energy_momentum/escape.in index 1cbcc07db..a050488ee 100644 --- a/examples/symba_energy_momentum/escape.in +++ b/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 diff --git a/examples/symba_energy_momentum/param.escape.in b/examples/symba_energy_momentum/param.escape.in index ee8628cca..6a6ba4367 100644 --- a/examples/symba_energy_momentum/param.escape.in +++ b/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,8 +14,8 @@ 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 @@ -23,7 +23,7 @@ CHK_QMIN -1.0 ! ignore this check 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 diff --git a/examples/symba_energy_momentum/param.sun.in b/examples/symba_energy_momentum/param.sun.in index 0e26d328f..dfc26b7e0 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 1.0 +TSTOP 3.0e-2 DT 1e-3 PL_IN sun.in TP_IN tp.in diff --git a/examples/symba_energy_momentum/sun.in b/examples/symba_energy_momentum/sun.in index a84e3b125..d281c7811 100644 --- a/examples/symba_energy_momentum/sun.in +++ b/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 diff --git a/src/coord/coord_h2b.f90 b/src/coord/coord_h2b.f90 index ba1569498..2ede3b080 100644 --- a/src/coord/coord_h2b.f90 +++ b/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(:) diff --git a/src/io/io_conservation_report.f90 b/src/io/io_conservation_report.f90 index 6a666e0ff..d240eb90a 100644 --- a/src/io/io_conservation_report.f90 +++ b/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 diff --git a/src/modules/module_interfaces.f90 b/src/modules/module_interfaces.f90 index 7b39d79b7..12ec7d02b 100644 --- a/src/modules/module_interfaces.f90 +++ b/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 diff --git a/src/symba/symba_discard_conserve_mtm.f90 b/src/symba/symba_discard_conserve_mtm.f90 index f0d54fb44..df1e2474e 100644 --- a/src/symba/symba_discard_conserve_mtm.f90 +++ b/src/symba/symba_discard_conserve_mtm.f90 @@ -12,9 +12,9 @@ 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, & @@ -22,7 +22,8 @@ subroutine symba_discard_conserve_mtm(param, swiftest_plA, ipl, lescape_body) 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