diff --git a/src/io/io.f90 b/src/io/io.f90 index 28ef8d143..f225c6523 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -46,29 +46,29 @@ module subroutine io_conservation_report(self, param, lterminal) Lorbit_now(:) = system%Lorbit(:) Lspin_now(:) = system%Lspin(:) Eorbit_now = ke_orbit_now + ke_spin_now + pe_now - Ltot_now(:) = system%Ltot(:) + param%Lescape(:) - GMtot_now = system%GMtot + param%GMescape + Ltot_now(:) = system%Ltot(:) + system%Lescape(:) + GMtot_now = system%GMtot + system%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(:) + system%Eorbit_orig = Eorbit_now + system%GMtot_orig = GMtot_now + system%Lorbit_orig(:) = Lorbit_now(:) + system%Lspin_orig(:) = Lspin_now(:) + system%Ltot_orig(:) = Ltot_now(:) param%lfirstenergy = .false. end if if (((param%out_type == REAL4_TYPE) .or. (param%out_type == REAL8_TYPE)) .and. (param%energy_out /= "")) then - write(EGYIU,EGYFMT, err = 667, iomsg = errmsg) param%t, Eorbit_now, param%Ecollisions, Ltot_now, GMtot_now + write(EGYIU,EGYFMT, err = 667, iomsg = errmsg) param%t, Eorbit_now, system%Ecollisions, Ltot_now, GMtot_now close(EGYIU, err = 667, iomsg = errmsg) end if if (.not.param%lfirstenergy) 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 + Lerror = norm2(Ltot_now(:) - system%Ltot_orig(:)) / norm2(system%Ltot_orig(:)) + Eorbit_error = (Eorbit_now - system%Eorbit_orig) / abs(system%Eorbit_orig) + Ecoll_error = system%Ecollisions / abs(system%Eorbit_orig) + Etotal_error = (Eorbit_now - system%Ecollisions - system%Eorbit_orig - system%Euntracked) / abs(system%Eorbit_orig) + Merror = (GMtot_now - system%GMtot_orig) / system%GMtot_orig if (lterminal) write(*, EGYTERMFMT) Lerror, Ecoll_error, Etotal_error, Merror if (abs(Merror) > 100 * epsilon(Merror)) then write(*,*) "Severe error! Mass not conserved! Halting!" @@ -281,6 +281,17 @@ module subroutine io_dump_system(self, param) dump_param%incbfile = trim(adjustl(DUMP_CB_FILE(idx))) dump_param%inplfile = trim(adjustl(DUMP_PL_FILE(idx))) dump_param%intpfile = trim(adjustl(DUMP_TP_FILE(idx))) + + dump_param%Eorbit_orig = self%Eorbit_orig + dump_param%GMtot_orig = self%GMtot_orig + dump_param%Ltot_orig(:) = self%Ltot_orig(:) + dump_param%Lorbit_orig(:) = self%Lorbit_orig(:) + dump_param%Lspin_orig(:) = self%Lspin_orig(:) + dump_param%GMescape = self%GMescape + dump_param%Ecollisions = self%Ecollisions + dump_param%Euntracked = self%Euntracked + dump_param%Lescape(:) = self%Lescape + else if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then dump_param%in_type = NETCDF_DOUBLE_TYPE dump_param%in_netcdf = trim(adjustl(DUMP_NC_FILE(idx))) @@ -995,15 +1006,17 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg) if (param%lenergy) then call io_param_writer_one("FIRSTENERGY", param%lfirstenergy, unit) - call io_param_writer_one("EORBIT_ORIG", param%Eorbit_orig, unit) - call io_param_writer_one("GMTOT_ORIG", param%GMtot_orig, unit) - call io_param_writer_one("LTOT_ORIG", param%Ltot_orig(:), unit) - call io_param_writer_one("LORBIT_ORIG", param%Lorbit_orig(:), unit) - call io_param_writer_one("LSPIN_ORIG", param%Lspin_orig(:), unit) - call io_param_writer_one("LESCAPE", param%Lescape(:), unit) - call io_param_writer_one("GMESCAPE",param%GMescape, unit) - call io_param_writer_one("ECOLLISIONS",param%Ecollisions, unit) - call io_param_writer_one("EUNTRACKED",param%Euntracked, unit) + if ((param%out_type == REAL8_TYPE) .or. (param%out_type == REAL4_TYPE)) then + call io_param_writer_one("EORBIT_ORIG", param%Eorbit_orig, unit) + call io_param_writer_one("GMTOT_ORIG", param%GMtot_orig, unit) + call io_param_writer_one("LTOT_ORIG", param%Ltot_orig(:), unit) + call io_param_writer_one("LORBIT_ORIG", param%Lorbit_orig(:), unit) + call io_param_writer_one("LSPIN_ORIG", param%Lspin_orig(:), unit) + call io_param_writer_one("LESCAPE", param%Lescape(:), unit) + call io_param_writer_one("GMESCAPE",param%GMescape, unit) + call io_param_writer_one("ECOLLISIONS",param%Ecollisions, unit) + call io_param_writer_one("EUNTRACKED",param%Euntracked, unit) + end if end if call io_param_writer_one("FIRSTKICK",param%lfirstkick, unit) call io_param_writer_one("MAXID",param%maxid, unit) @@ -1372,13 +1385,22 @@ module subroutine io_read_in_system(self, param) allocate(tmp_param, source=param) tmp_param%outfile = param%in_netcdf tmp_param%out_form = param%in_form - ierr = self%read_frame(tmp_param%nciu, tmp_param) + ierr = self%read_frame(tmp_param%nciu, tmp_param) deallocate(tmp_param) if (ierr /=0) call util_exit(FAILURE) else call self%cb%read_in(param) call self%pl%read_in(param) call self%tp%read_in(param) + ! Copy over param file variable inputs + self%Eorbit_orig = param%Eorbit_orig + self%GMtot_orig = param%GMtot_orig + self%Ltot_orig(:) = param%Ltot_orig(:) + self%Lorbit_orig(:) = param%Lorbit_orig(:) + self%Lspin_orig(:) = param%Lspin_orig(:) + self%Lescape(:) = param%Lescape(:) + self%Ecollisions = param%Ecollisions + self%Euntracked = param%Euntracked end if return diff --git a/src/main/swiftest_driver.f90 b/src/main/swiftest_driver.f90 index eb2a375a1..5fdb43c8b 100644 --- a/src/main/swiftest_driver.f90 +++ b/src/main/swiftest_driver.f90 @@ -65,8 +65,8 @@ program swiftest_driver old_t_final = nbody_system%get_old_t_final_netcdf(param) else old_t_final = t0 - if (istep_out > 0) call nbody_system%write_frame(param) if (param%lenergy) call nbody_system%conservation_report(param, lterminal=.false.) ! This will save the initial values of energy and momentum + if (istep_out > 0) call nbody_system%write_frame(param) end if !> Define the maximum number of threads diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 65e05ae87..643d045f9 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -152,12 +152,12 @@ 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) :: GMtot_orig = 0.0_DP !! Initial system mass 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) :: 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) :: 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 logical :: lfirstenergy = .true. !! This is the first time computing energe @@ -424,6 +424,15 @@ module swiftest_classes 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) :: Ltot = 0.0_DP !! System angular momentum vector + real(DP) :: Eorbit_orig = 0.0_DP !! Initial orbital energy + real(DP) :: GMtot_orig = 0.0_DP !! Initial system mass + 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) :: 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 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 @@ -442,7 +451,6 @@ module swiftest_classes procedure :: read_frame_netcdf => netcdf_read_frame_system !! Read in a frame of input data from file procedure :: write_frame_netcdf => netcdf_write_frame_system !! Write a frame of input data from file procedure :: write_hdr_bin => io_write_hdr_system !! Write a header for an output frame in Fortran binary format - !procedure :: read_hdr_bin => io_read_hdr !! Read a header for an output frame in Fortran binary format procedure :: read_hdr_netcdf => netcdf_read_hdr_system !! Read a header for an output frame in NetCDF format procedure :: write_hdr_netcdf => netcdf_write_hdr_system !! Write a header for an output frame in NetCDF format procedure :: read_in => io_read_in_system !! Reads the initial conditions for an nbody system @@ -961,7 +969,7 @@ end subroutine netcdf_flush module function netcdf_get_old_t_final_system(self, param) result(old_t_final) implicit none - class(swiftest_nbody_system), intent(in) :: self !! Swiftest nbody system object + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP) :: old_t_final !! Final time from last run end function netcdf_get_old_t_final_system diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index 384a7f437..f8363d976 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -58,23 +58,58 @@ module function netcdf_get_old_t_final_system(self, param) result(old_t_final) !! implicit none ! Arguments - class(swiftest_nbody_system), intent(in) :: self + class(swiftest_nbody_system), intent(inout) :: self class(swiftest_parameters), intent(inout) :: param ! Result real(DP) :: old_t_final ! Internals - integer(I4B) :: itmax - real(DP), dimension(:), allocatable :: tvals - + integer(I4B) :: itmax, idmax + real(DP), dimension(:), allocatable :: vals + real(DP), dimension(1) :: val + real(DP) :: KE_orb_orig, KE_spin_orig, PE_orig, Ltmp call param%nciu%open(param) call check( nf90_inquire_dimension(param%nciu%ncid, param%nciu%time_dimid, len=itmax) ) - allocate(tvals(itmax)) - call check( nf90_get_var(param%nciu%ncid, param%nciu%time_varid, tvals(:)) ) + call check( nf90_inquire_dimension(param%nciu%ncid, param%nciu%id_dimid, len=idmax) ) + allocate(vals(idmax)) + call check( nf90_get_var(param%nciu%ncid, param%nciu%time_varid, val, start=[itmax], count=[1]) ) + + old_t_final = val(1) + + if (param%lenergy) then + call check( nf90_get_var(param%nciu%ncid, param%nciu%KE_orb_varid, val, start=[1], count=[1]) ) + KE_orb_orig = val(1) + + call check( nf90_get_var(param%nciu%ncid, param%nciu%KE_spin_varid, val, start=[1], count=[1]) ) + KE_spin_orig = val(1) + + call check( nf90_get_var(param%nciu%ncid, param%nciu%PE_varid, val, start=[1], count=[1]) ) + PE_orig = val(1) - old_t_final = tvals(itmax) + self%Eorbit_orig = KE_orb_orig + KE_spin_orig + PE_orig + + call check( nf90_get_var(param%nciu%ncid, param%nciu%L_orbx_varid, val, start=[1], count=[1]) ) + self%Lorbit_orig(1) = val(1) + call check( nf90_get_var(param%nciu%ncid, param%nciu%L_orby_varid, val, start=[1], count=[1]) ) + self%Lorbit_orig(2) = val(1) + call check( nf90_get_var(param%nciu%ncid, param%nciu%L_orbz_varid, val, start=[1], count=[1]) ) + self%Lorbit_orig(3) = val(1) + + call check( nf90_get_var(param%nciu%ncid, param%nciu%L_spinx_varid, val, start=[1], count=[1]) ) + self%Lspin_orig(1) = val(1) + call check( nf90_get_var(param%nciu%ncid, param%nciu%L_spiny_varid, val, start=[1], count=[1]) ) + self%Lspin_orig(2) = val(1) + call check( nf90_get_var(param%nciu%ncid, param%nciu%L_spinz_varid, val, start=[1], count=[1]) ) + self%Lspin_orig(3) = val(1) + + self%Ltot_orig(:) = self%Lorbit_orig(:) + self%Lspin_orig(:) + + call check( nf90_get_var(param%nciu%ncid, param%nciu%Gmass_varid, vals, start=[1,1], count=[idmax,1]) ) + self%GMtot_orig = vals(1) + sum(vals(2:idmax), vals(2:idmax) == vals(2:idmax)) + + end if - deallocate(tvals) + deallocate(vals) return end function netcdf_get_old_t_final_system @@ -607,21 +642,21 @@ module subroutine netcdf_read_hdr_system(self, iu, param) call check( nf90_get_var(iu%ncid, iu%ntp_varid, self%tp%nbody, start=[tslot]) ) if (param%lenergy) then - call check( nf90_get_var(iu%ncid, iu%KE_orb_varid, self%ke_orbit, start=[tslot]) ) - call check( nf90_get_var(iu%ncid, iu%KE_spin_varid, self%ke_spin, start=[tslot]) ) - call check( nf90_get_var(iu%ncid, iu%PE_varid, self%pe, start=[tslot]) ) - call check( nf90_get_var(iu%ncid, iu%L_orbx_varid, self%Lorbit(1), start=[tslot]) ) - call check( nf90_get_var(iu%ncid, iu%L_orby_varid, self%Lorbit(2), start=[tslot]) ) - call check( nf90_get_var(iu%ncid, iu%L_orbz_varid, self%Lorbit(3), start=[tslot]) ) - call check( nf90_get_var(iu%ncid, iu%L_spinx_varid, self%Lspin(1), start=[tslot]) ) - call check( nf90_get_var(iu%ncid, iu%L_spiny_varid, self%Lspin(2), start=[tslot]) ) - call check( nf90_get_var(iu%ncid, iu%L_spinz_varid, self%Lspin(3), start=[tslot]) ) - call check( nf90_get_var(iu%ncid, iu%L_escapex_varid, param%Lescape(1), start=[tslot]) ) - call check( nf90_get_var(iu%ncid, iu%L_escapey_varid, param%Lescape(2), start=[tslot]) ) - call check( nf90_get_var(iu%ncid, iu%L_escapez_varid, param%Lescape(3), start=[tslot]) ) - call check( nf90_get_var(iu%ncid, iu%Ecollisions_varid, param%Ecollisions, start=[tslot]) ) - call check( nf90_get_var(iu%ncid, iu%Euntracked_varid, param%Euntracked, start=[tslot]) ) - call check( nf90_get_var(iu%ncid, iu%GMescape_varid, param%GMescape, start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%KE_orb_varid, self%ke_orbit, start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%KE_spin_varid, self%ke_spin, start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%PE_varid, self%pe, start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%L_orbx_varid, self%Lorbit(1), start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%L_orby_varid, self%Lorbit(2), start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%L_orbz_varid, self%Lorbit(3), start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%L_spinx_varid, self%Lspin(1), start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%L_spiny_varid, self%Lspin(2), start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%L_spinz_varid, self%Lspin(3), start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%L_escapex_varid, self%Lescape(1), start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%L_escapey_varid, self%Lescape(2), start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%L_escapez_varid, self%Lescape(3), start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%Ecollisions_varid, self%Ecollisions, start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%Euntracked_varid, self%Euntracked, start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%GMescape_varid, self%GMescape, start=[tslot]) ) end if return @@ -1077,12 +1112,12 @@ module subroutine netcdf_write_hdr_system(self, iu, param) call check( nf90_put_var(iu%ncid, iu%L_spinx_varid, self%Lspin(1), start=[tslot]) ) call check( nf90_put_var(iu%ncid, iu%L_spiny_varid, self%Lspin(2), start=[tslot]) ) call check( nf90_put_var(iu%ncid, iu%L_spinz_varid, self%Lspin(3), start=[tslot]) ) - call check( nf90_put_var(iu%ncid, iu%L_escapex_varid, param%Lescape(1), start=[tslot]) ) - call check( nf90_put_var(iu%ncid, iu%L_escapey_varid, param%Lescape(2), start=[tslot]) ) - call check( nf90_put_var(iu%ncid, iu%L_escapez_varid, param%Lescape(3), start=[tslot]) ) - call check( nf90_put_var(iu%ncid, iu%Ecollisions_varid, param%Ecollisions, start=[tslot]) ) - call check( nf90_put_var(iu%ncid, iu%Euntracked_varid, param%Euntracked, start=[tslot]) ) - call check( nf90_put_var(iu%ncid, iu%GMescape_varid, param%GMescape, start=[tslot]) ) + call check( nf90_put_var(iu%ncid, iu%L_escapex_varid, self%Lescape(1), start=[tslot]) ) + call check( nf90_put_var(iu%ncid, iu%L_escapey_varid, self%Lescape(2), start=[tslot]) ) + call check( nf90_put_var(iu%ncid, iu%L_escapez_varid, self%Lescape(3), start=[tslot]) ) + call check( nf90_put_var(iu%ncid, iu%Ecollisions_varid, self%Ecollisions, start=[tslot]) ) + call check( nf90_put_var(iu%ncid, iu%Euntracked_varid, self%Euntracked, start=[tslot]) ) + call check( nf90_put_var(iu%ncid, iu%GMescape_varid, self%GMescape, start=[tslot]) ) end if return diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 74a12413b..3a319d820 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -192,8 +192,8 @@ module function symba_collision_casemerge(system, param, colliders, frag) resul pe = pe - pl%Gmass(i) * pl%mass(j) / norm2(pl%xb(:, i) - pl%xb(:, j)) end do end do - param%Ecollisions = param%Ecollisions + pe - param%Euntracked = param%Euntracked - pe + system%Ecollisions = system%Ecollisions + pe + system%Euntracked = system%Euntracked - pe ! Update any encounter lists that have the removed bodies in them so that they instead point to the new do k = 1, system%plplenc_list%nenc @@ -1025,7 +1025,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 - param%Ecollisions = param%Ecollisions + (Eorbit_after - Eorbit_before) + system%Ecollisions = system%Ecollisions + (Eorbit_after - Eorbit_before) end if end select diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index 592b146b8..0ae89587f 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -114,7 +114,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 - param%GMescape = param%GMescape + pl%Gmass(ipl) + system%GMescape = system%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)) @@ -137,8 +137,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(:)) - 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) + 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) else xcom(:) = (pl%mass(ipl) * pl%xb(:, ipl) + cb%mass * cb%xb(:)) / (cb%mass + pl%mass(ipl)) @@ -172,11 +172,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 - param%Ecollisions = param%Ecollisions + ke_orbit + ke_spin + pe - param%Euntracked = param%Euntracked - (ke_orbit + ke_spin + pe) + system%Ecollisions = system%Ecollisions + ke_orbit + ke_spin + pe + system%Euntracked = system%Euntracked - (ke_orbit + ke_spin + pe) else - param%Ecollisions = param%Ecollisions + pe - param%Euntracked = param%Euntracked - pe + system%Ecollisions = system%Ecollisions + pe + system%Euntracked = system%Euntracked - pe end if end select @@ -358,7 +358,7 @@ module subroutine symba_discard_pl(self, system, param) if (param%lenergy) then call system%get_energy_and_momentum(param) Eorbit_after = system%te - param%Ecollisions = param%Ecollisions + (Eorbit_after - Eorbit_before) + system%Ecollisions = system%Ecollisions + (Eorbit_after - Eorbit_before) end if end associate