diff --git a/src/io/io.f90 b/src/io/io.f90 index 37a69a6ec..670d98db8 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -14,11 +14,10 @@ module subroutine io_conservation_report(self, param, lterminal) ! Internals real(DP), dimension(NDIM) :: Ltot_now, Lorbit_now, Lspin_now real(DP), dimension(NDIM), save :: Ltot_last, Lorbit_last, Lspin_last - real(DP), save :: ke_orbit_last, ke_spin_last, pe_last, Eorbit_last real(DP) :: ke_orbit_now, ke_spin_now, pe_now, Eorbit_now real(DP) :: Eorbit_error, Etotal_error, Ecoll_error - real(DP) :: GMtot_now, Merror - real(DP) :: Lmag_now, Lerror + real(DP) :: GMtot_now + real(DP) :: Lerror, Merror character(len=STRMAX) :: errmsg 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")' @@ -28,12 +27,9 @@ module subroutine io_conservation_report(self, param, lterminal) "; D(Eorbit+Ecollisions)/|E0| = ", ES12.5, & "; DM/M0 = ", ES12.5)' - associate(system => self, pl => self%pl, cb => self%cb, npl => self%pl%nbody, Ecollisions => self%Ecollisions, Lescape => self%Lescape, GMescape => self%GMescape, & - Euntracked => self%Euntracked, Eorbit_orig => param%Eorbit_orig, GMtot_orig => param%GMtot_orig, & - Ltot_orig => param%Ltot_orig(:), Lmag_orig => param%Lmag_orig, Lorbit_orig => param%Lorbit_orig(:), Lspin_orig => param%Lspin_orig(:), & - lfirst => param%lfirstenergy) + associate(system => self, pl => self%pl, cb => self%cb, npl => self%pl%nbody) if (param%energy_out /= "") then - if (lfirst .and. (param%out_stat /= "OLD")) then + if (param%lfirstenergy .and. (param%out_stat /= "OLD")) then open(unit = EGYIU, file = param%energy_out, form = "formatted", status = "replace", action = "write", err = 667, iomsg = errmsg) write(EGYIU,EGYHEADER, err = 667, iomsg = errmsg) else @@ -46,48 +42,34 @@ module subroutine io_conservation_report(self, param, lterminal) ke_orbit_now = system%ke_orbit ke_spin_now = system%ke_spin pe_now = system%pe - Lorbit_now = system%Lorbit - Lspin_now = system%Lspin + Lorbit_now(:) = system%Lorbit(:) + Lspin_now(:) = system%Lspin(:) Eorbit_now = ke_orbit_now + ke_spin_now + pe_now - Ltot_now(:) = Lorbit_now(:) + Lspin_now(:) + Lescape(:) - GMtot_now = cb%Gmass + sum(pl%Gmass(1:npl)) + system%GMescape - if (lfirst) then - Eorbit_orig = Eorbit_now - GMtot_orig = GMtot_now - Lorbit_orig(:) = Lorbit_now(:) - Lspin_orig(:) = Lspin_now(:) - Ltot_orig(:) = Ltot_now(:) - Lmag_orig = norm2(Ltot_orig(:)) - lfirst = .false. + Ltot_now(:) = system%Ltot(:) + param%Lescape(:) + GMtot_now = system%GMtot + param%GMescape + + if (param%lfirstenergy) then + param%Eorbit_orig = Eorbit_now + param%GMtot_orig = GMtot_now + param%Lorbit_orig(:) = Lorbit_now(:) + param%Lspin_orig(:) = Lspin_now(:) + param%Ltot_orig(:) = Ltot_now(:) + param%lfirstenergy = .false. end if if (param%energy_out /= "") then - write(EGYIU,EGYFMT, err = 667, iomsg = errmsg) param%t, Eorbit_now, Ecollisions, Ltot_now, GMtot_now + write(EGYIU,EGYFMT, err = 667, iomsg = errmsg) param%t, Eorbit_now, param%Ecollisions, Ltot_now, GMtot_now close(EGYIU, err = 667, iomsg = errmsg) end if - if (.not.lfirst .and. lterminal) then - Lmag_now = norm2(Ltot_now) - Lerror = norm2(Ltot_now - Ltot_orig) / Lmag_orig - Eorbit_error = (Eorbit_now - Eorbit_orig) / abs(Eorbit_orig) - Ecoll_error = Ecollisions / abs(Eorbit_orig) - Etotal_error = (Eorbit_now - Ecollisions - Eorbit_orig - Euntracked) / abs(Eorbit_orig) - Merror = (GMtot_now - GMtot_orig) / GMtot_orig + + if (.not.param%lfirstenergy .and. lterminal) then + Lerror = norm2(Ltot_now(:) - param%Ltot_orig(:)) / norm2(param%Ltot_orig(:)) + Eorbit_error = (Eorbit_now - param%Eorbit_orig) / abs(param%Eorbit_orig) + Ecoll_error = param%Ecollisions / abs(param%Eorbit_orig) + Etotal_error = (Eorbit_now - param%Ecollisions - param%Eorbit_orig - param%Euntracked) / abs(param%Eorbit_orig) + Merror = (GMtot_now - param%GMtot_orig) / param%GMtot_orig write(*, egytermfmt) Lerror, Ecoll_error, Etotal_error, Merror - ! if (Ecoll_error > 0.0_DP) then - ! write(*,*) 'Something has gone wrong! Collisional energy should not be positive!' - ! write(*,*) 'dke_orbit: ',(ke_orbit_now - ke_orbit_last) / abs(Eorbit_orig) - ! write(*,*) 'dke_spin : ',(ke_spin_now - ke_spin_last) / abs(Eorbit_orig) - ! write(*,*) 'dpe : ',(pe_now - pe_last) / abs(Eorbit_orig) - ! write(*,*) - ! end if end if - ke_orbit_last = ke_orbit_now - ke_spin_last = ke_spin_now - pe_last = pe_now - Eorbit_last = Eorbit_now - Lorbit_last(:) = Lorbit_now(:) - Lspin_last(:) = Lspin_now(:) - Ltot_last(:) = Ltot_now(:) end associate return @@ -212,6 +194,7 @@ module subroutine io_dump_system(self, param) dump_param%out_stat = 'APPEND' dump_param%in_type = REAL8_TYPE dump_param%T0 = param%t + call dump_param%dump(param_file_name) call self%cb%dump(dump_param) @@ -536,7 +519,6 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) param_value = io_get_token(line, ifirst, ilast, iostat) read(param_value, *, err = 667, iomsg = iomsg) param%Ltot_orig(i) end do - param%Lmag_orig = norm2(param%Ltot_orig(:)) case("LORBIT_ORIG") read(param_value, *, err = 667, iomsg = iomsg) param%Lorbit_orig(1) do i = 2, NDIM diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 465170151..f08e8eb05 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -61,11 +61,9 @@ module swiftest_classes ! Initial values to pass to the energy report subroutine (usually only used in the case of a restart, otherwise these will be updated with initial conditions values) real(DP) :: Eorbit_orig = 0.0_DP !! Initial orbital energy real(DP) :: GMtot_orig = 0.0_DP !! Initial system mass - real(DP) :: Lmag_orig = 0.0_DP !! Initial total angular momentum magnitude real(DP), dimension(NDIM) :: Ltot_orig = 0.0_DP !! Initial total angular momentum vector real(DP), dimension(NDIM) :: Lorbit_orig = 0.0_DP !! Initial orbital angular momentum real(DP), dimension(NDIM) :: Lspin_orig = 0.0_DP !! Initial spin angular momentum vector - real(DP), dimension(NDIM) :: Ltot = 0.0_DP !! System angular momentum vector real(DP), dimension(NDIM) :: Lescape = 0.0_DP !! Angular momentum of bodies that escaped the system (used for bookeeping) real(DP) :: GMescape = 0.0_DP !! Mass of bodies that escaped the system (used for bookeeping) real(DP) :: Ecollisions = 0.0_DP !! Energy lost from system due to collisions @@ -283,17 +281,14 @@ module swiftest_classes class(swiftest_tp), allocatable :: tp !! Test particle data structure class(swiftest_tp), allocatable :: tp_discards !! Discarded test particle data structure class(swiftest_pl), allocatable :: pl_discards !! Discarded massive body particle data structure - real(DP) :: Gmtot = 0.0_DP !! Total system mass - used for barycentric coordinate conversion + real(DP) :: GMtot = 0.0_DP !! Total system mass - used for barycentric coordinate conversion real(DP) :: ke_orbit = 0.0_DP !! System orbital kinetic energy real(DP) :: ke_spin = 0.0_DP !! System spin kinetic energy real(DP) :: pe = 0.0_DP !! System potential energy real(DP) :: te = 0.0_DP !! System total energy real(DP), dimension(NDIM) :: Lorbit = 0.0_DP !! System orbital angular momentum vector real(DP), dimension(NDIM) :: Lspin = 0.0_DP !! System spin angular momentum vector - real(DP), dimension(NDIM) :: Lescape = 0.0_DP !! Angular momentum of bodies that escaped the system (used for bookeeping) - real(DP) :: GMescape = 0.0_DP !! Mass of bodies that escaped the system (used for bookeeping) - real(DP) :: Ecollisions = 0.0_DP !! Energy lost from system due to collisions - real(DP) :: Euntracked = 0.0_DP !! Energy gained from system due to escaped bodies + real(DP), dimension(NDIM) :: Ltot = 0.0_DP !! System angular momentum vector logical :: lbeg !! True if this is the beginning of a step. This is used so that test particle steps can be calculated !! separately from massive bodies. Massive body variables are saved at half steps, and passed to !! the test particles diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 9384a11a6..f2b9134d8 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -215,14 +215,14 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param) implicit none class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions end subroutine symba_collision_resolve_fragmentations module subroutine symba_collision_resolve_mergers(self, system, param) implicit none class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions end subroutine symba_collision_resolve_mergers module subroutine symba_collision_resolve_plplenc(self, system, param, t, dt, irec) @@ -309,7 +309,7 @@ end function symba_encounter_check_tp module function symba_collision_casedisruption(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) result(status) implicit none class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision real(DP), dimension(:), intent(inout) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision @@ -321,7 +321,7 @@ end function symba_collision_casedisruption module function symba_collision_casehitandrun(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) result(status) implicit none class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision real(DP), dimension(:), intent(inout) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision @@ -333,7 +333,7 @@ end function symba_collision_casehitandrun module function symba_collision_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip) result(status) implicit none class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision real(DP), dimension(:,:), intent(in) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision real(DP), dimension(:), intent(in) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collisio @@ -343,7 +343,7 @@ end function symba_collision_casemerge module function symba_collision_casesupercatastrophic(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) result(status) implicit none class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision real(DP), dimension(:), intent(inout) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index b00a4a8e6..5ba08558c 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -13,7 +13,7 @@ module function symba_collision_casedisruption(system, param, family, x, v, mass implicit none ! Arguments class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision real(DP), dimension(:), intent(inout) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision @@ -106,7 +106,7 @@ module function symba_collision_casehitandrun(system, param, family, x, v, mass, implicit none ! Arguments class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision real(DP), dimension(:), intent(inout) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision @@ -213,7 +213,7 @@ module function symba_collision_casemerge(system, param, family, x, v, mass, rad implicit none ! Arguments class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision real(DP), dimension(:,:), intent(in) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision real(DP), dimension(:), intent(in) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision @@ -269,7 +269,7 @@ module function symba_collision_casemerge(system, param, family, x, v, mass, rad ! Assume prinicpal axis rotation on 3rd Ip axis rot_frag(:,1) = L_spin_new(:) / (Ip_frag(3,1) * m_frag(1) * rad_frag(1)**2) else ! If spin is not enabled, we will consider the lost pre-collision angular momentum as "escaped" and add it to our bookkeeping variable - system%Lescape(:) = system%Lescape(:) + L_orb_old(:) + param%Lescape(:) = param%Lescape(:) + L_orb_old(:) end if ! Keep track of the component of potential energy due to the pre-impact family for book-keeping @@ -280,8 +280,8 @@ module function symba_collision_casemerge(system, param, family, x, v, mass, rad pe = pe - pl%Gmass(i) * pl%mass(j) / norm2(pl%xb(:, i) - pl%xb(:, j)) end do end do - system%Ecollisions = system%Ecollisions + pe - system%Euntracked = system%Euntracked - pe + param%Ecollisions = param%Ecollisions + pe + param%Euntracked = param%Euntracked - pe status = MERGED call symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, status) @@ -300,7 +300,7 @@ module function symba_collision_casesupercatastrophic(system, param, family, x, implicit none ! Arguments class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision real(DP), dimension(:), intent(inout) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision @@ -930,7 +930,7 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param) ! Arguments class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions ! Internals ! Internals integer(I4B), dimension(:), allocatable :: family !! List of indices of all bodies inovlved in the collision @@ -1023,7 +1023,7 @@ module subroutine symba_collision_resolve_mergers(self, system, param) ! Arguments class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions ! Internals integer(I4B), dimension(:), allocatable :: family !! List of indices of all bodies inovlved in the collision integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision @@ -1119,7 +1119,7 @@ module subroutine symba_collision_resolve_plplenc(self, system, param, t, dt, ir if (param%lenergy) then call system%get_energy_and_momentum(param) Eorbit_after = system%te - system%Ecollisions = system%Ecollisions + (Eorbit_after - Eorbit_before) + param%Ecollisions = param%Ecollisions + (Eorbit_after - Eorbit_before) end if end select diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index 32879a9a5..193461bc7 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -96,7 +96,7 @@ subroutine symba_discard_conserve_mtm(pl, system, param, ipl, lescape_body) ! Add the pre-collision ke of the central body to the records ! Add planet mass to central body accumulator if (lescape_body) then - system%GMescape = system%GMescape + pl%Gmass(ipl) + param%GMescape = param%GMescape + pl%Gmass(ipl) do i = 1, pl%nbody if (i == ipl) cycle pe = pe - pl%Gmass(i) * pl%mass(ipl) / norm2(pl%xb(:, ipl) - pl%xb(:, i)) @@ -119,8 +119,8 @@ subroutine symba_discard_conserve_mtm(pl, system, param, ipl, lescape_body) Ltot(:) = Ltot(:) - Lpl(:) end do Ltot(:) = Ltot(:) - cb%mass * (cb%xb(:) .cross. cb%vb(:)) - system%Lescape(:) = system%Lescape(:) + Ltot(:) - if (param%lrotation) system%Lescape(:) = system%Lescape + pl%mass(ipl) * pl%radius(ipl)**2 * pl%Ip(3, ipl) * pl%rot(:, ipl) + param%Lescape(:) = param%Lescape(:) + Ltot(:) + if (param%lrotation) param%Lescape(:) = param%Lescape + pl%mass(ipl) * pl%radius(ipl)**2 * pl%Ip(3, ipl) * pl%rot(:, ipl) else xcom(:) = (pl%mass(ipl) * pl%xb(:, ipl) + cb%mass * cb%xb(:)) / (cb%mass + pl%mass(ipl)) @@ -154,11 +154,11 @@ subroutine symba_discard_conserve_mtm(pl, system, param, ipl, lescape_body) ! 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 - system%Ecollisions = system%Ecollisions + ke_orbit + ke_spin + pe - system%Euntracked = system%Euntracked - (ke_orbit + ke_spin + pe) + param%Ecollisions = param%Ecollisions + ke_orbit + ke_spin + pe + param%Euntracked = param%Euntracked - (ke_orbit + ke_spin + pe) else - system%Ecollisions = system%Ecollisions + pe - system%Euntracked = system%Euntracked - pe + param%Ecollisions = param%Ecollisions + pe + param%Euntracked = param%Euntracked - pe end if end select @@ -313,7 +313,7 @@ module subroutine symba_discard_pl(self, system, param) if (param%lenergy) then call system%get_energy_and_momentum(param) Eorbit_after = system%te - system%Ecollisions = system%Ecollisions + (Eorbit_after - Eorbit_before) + param%Ecollisions = param%Ecollisions + (Eorbit_after - Eorbit_before) end if end associate diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index dffdc29ae..4243070ac 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -28,6 +28,7 @@ module subroutine util_get_energy_momentum_system(self, param) associate(system => self, pl => self%pl, npl => self%pl%nbody, cb => self%cb) system%Lorbit(:) = 0.0_DP system%Lspin(:) = 0.0_DP + system%Ltot(:) = 0.0_DP system%ke_orbit = 0.0_DP system%ke_spin = 0.0_DP @@ -38,8 +39,10 @@ module subroutine util_get_energy_momentum_system(self, param) Lplspinx(:) = 0.0_DP Lplspiny(:) = 0.0_DP Lplspinz(:) = 0.0_DP + lstatus(1:npl) = pl%status(1:npl) /= INACTIVE + system%GMtot = cb%Gmass + sum(pl%Gmass(1:npl), lstatus(1:npl)) kecb = cb%mass * dot_product(cb%vb(:), cb%vb(:)) Lcborbit(:) = cb%mass * (cb%xb(:) .cross. cb%vb(:)) @@ -121,6 +124,7 @@ module subroutine util_get_energy_momentum_system(self, param) end if system%te = system%ke_orbit + system%ke_spin + system%pe + system%Ltot(:) = system%Lorbit(:) + system%Lspin(:) end associate return