diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index d879eec53..5d92a1a3e 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -22,7 +22,7 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, real(DP), intent(inout) :: Qloss !! Energy lost during the collision logical, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? ! Internals - real(DP) :: mscale, rscale, vscale, tscale, Lscale, Escale ! Scale factors that reduce quantities to O(~1) in the collisional system + real(DP) :: mscale, dscale, vscale, tscale, Lscale, Escale ! Scale factors that reduce quantities to O(~1) in the collisional system real(DP) :: mtot real(DP), dimension(NDIM) :: xcom, vcom integer(I4B) :: ii @@ -46,6 +46,7 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, integer(I4B), parameter :: MAXTRY = 3000 integer(I4B), parameter :: TANTRY = 3 logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes, fpe_quiet_modes + class(swiftest_parameters), allocatable :: tmpparam if (nfrag < NFRAG_MIN) then write(*,*) "symba_frag_pos needs at least ",NFRAG_MIN," fragments, but only ",nfrag," were given." @@ -59,7 +60,7 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, f_spin = 0.05_DP mscale = 1.0_DP - rscale = 1.0_DP + dscale = 1.0_DP vscale = 1.0_DP tscale = 1.0_DP Lscale = 1.0_DP @@ -166,19 +167,19 @@ subroutine set_scale_factors() ! Set scale factors Escale = 0.5_DP * (mass(1) * dot_product(v(:,1), v(:,1)) + mass(2) * dot_product(v(:,2), v(:,2))) - rscale = sum(radius(:)) - mscale = sqrt(Escale * rscale / param%GU) + dscale = sum(radius(:)) + mscale = sqrt(Escale * dscale / param%GU) vscale = sqrt(Escale / mscale) - tscale = rscale / vscale - Lscale = mscale * rscale * vscale + tscale = dscale / vscale + Lscale = mscale * dscale * vscale - xcom(:) = xcom(:) / rscale + xcom(:) = xcom(:) / dscale vcom(:) = vcom(:) / vscale mtot = mtot / mscale mass = mass / mscale - radius = radius / rscale - x = x / rscale + radius = radius / dscale + x = x / dscale v = v / vscale L_spin = L_spin / Lscale do i = 1, 2 @@ -186,7 +187,7 @@ subroutine set_scale_factors() end do m_frag = m_frag / mscale - rad_frag = rad_frag / rscale + rad_frag = rad_frag / dscale Qloss = Qloss / Escale return @@ -201,13 +202,13 @@ subroutine restore_scale_factors() call ieee_set_halting_mode(IEEE_ALL,.false.) ! Restore scale factors - xcom(:) = xcom(:) * rscale + xcom(:) = xcom(:) * dscale vcom(:) = vcom(:) * vscale mtot = mtot * mscale mass = mass * mscale - radius = radius * rscale - x = x * rscale + radius = radius * dscale + x = x * dscale v = v * vscale L_spin = L_spin * Lscale do i = 1, 2 @@ -215,9 +216,9 @@ subroutine restore_scale_factors() end do m_frag = m_frag * mscale - rad_frag = rad_frag * rscale + rad_frag = rad_frag * dscale rot_frag = rot_frag / tscale - x_frag = x_frag * rscale + x_frag = x_frag * dscale v_frag = v_frag * vscale Qloss = Qloss * Escale @@ -240,7 +241,7 @@ subroutine restore_scale_factors() Lmag_after = Lmag_after * Lscale mscale = 1.0_DP - rscale = 1.0_DP + dscale = 1.0_DP vscale = 1.0_DP tscale = 1.0_DP Lscale = 1.0_DP @@ -314,7 +315,6 @@ subroutine calculate_system_energy(linclude_fragments) logical, dimension(:), allocatable :: ltmp logical :: lk_plpl class(swiftest_nbody_system), allocatable :: tmpsys - class(swiftest_parameters), allocatable :: tmpparam ! Because we're making a copy of symba_pl with the excludes/fragments appended, we need to deallocate the ! big k_plpl array and recreate it when we're done, otherwise we run the risk of blowing up the memory by @@ -344,12 +344,12 @@ subroutine calculate_system_energy(linclude_fragments) ltmp(1:npl) = .true. call tmpsys%pl%setup(npl_new, param) call tmpsys%pl%fill(pl, ltmp) - call tmpsys%rescale(tmpparam, mscale, rscale, tscale) + call tmpsys%rescale(tmpparam, mscale, dscale, tscale) if (linclude_fragments) then ! Append the fragments if they are included ! Energy calculation requires the fragments to be in the system barcyentric frame, s tmpsys%pl%mass(npl+1:npl_new) = m_frag(1:nfrag) - tmpsys%pl%Gmass(npl+1:npl_new) = m_frag(1:nfrag) * param%GU * rscale**3 * mscale / tscale**2 + tmpsys%pl%Gmass(npl+1:npl_new) = m_frag(1:nfrag) * param%GU * dscale**3 * mscale / tscale**2 tmpsys%pl%radius(npl+1:npl_new) = rad_frag(1:nfrag) tmpsys%pl%xb(:,npl+1:npl_new) = xb_frag(:,1:nfrag) tmpsys%pl%vb(:,npl+1:npl_new) = vb_frag(:,1:nfrag)