From c02ffb64269e7c03a3471b6a8761d7833668507c Mon Sep 17 00:00:00 2001 From: Andrew Benson Date: Fri, 15 Sep 2023 18:33:55 +0000 Subject: [PATCH] fix: Correct enclosed mass in dimensionfull Hernquist profiles --- .../mass_distributions.spherical.Hernquist.F90 | 4 ++-- source/tests.mass_distributions.F90 | 18 +++++++++++++++--- 2 files changed, 17 insertions(+), 5 deletions(-) diff --git a/source/mass_distributions.spherical.Hernquist.F90 b/source/mass_distributions.spherical.Hernquist.F90 index d1a09d30c1..fb07c116fa 100644 --- a/source/mass_distributions.spherical.Hernquist.F90 +++ b/source/mass_distributions.spherical.Hernquist.F90 @@ -268,9 +268,9 @@ double precision function hernquistMassEnclosedBySphere(self,radius,componentTyp fractionalRadius=radius/self%scaleLength if (fractionalRadius > fractionalRadiusLarge) then ! For very large radius approximate the mass enclosed as the total mass. - hernquistMassEnclosedBySphere=2.0d0*Pi*self%densityNormalization + hernquistMassEnclosedBySphere=self%mass else - hernquistMassEnclosedBySphere=2.0d0*Pi*self%densityNormalization*fractionalRadius**2/(1.0d0+fractionalRadius)**2 + hernquistMassEnclosedBySphere=self%mass*fractionalRadius**2/(1.0d0+fractionalRadius)**2 end if return end function hernquistMassEnclosedBySphere diff --git a/source/tests.mass_distributions.F90 b/source/tests.mass_distributions.F90 index f693dd1d71..d3b2ade72d 100644 --- a/source/tests.mass_distributions.F90 +++ b/source/tests.mass_distributions.F90 @@ -83,9 +83,21 @@ program Test_Mass_Distributions select type (massDistribution_) class is (massDistributionSpherical) position=[1.0d0,0.0d0,0.0d0] - call Assert("Half mass radius" ,massDistribution_%radiusHalfMass ( ), 1.0d0+sqrt(2.0d0),absTol=1.0d-6) - call Assert("Mass within scale radius" ,massDistribution_%massEnclosedBySphere(1.0d0 ), 0.25d0 ,absTol=1.0d-6) - call Assert("Potential at scale radius",massDistribution_%potential (position),-0.50d0 ,absTol=1.0d-6) + call Assert("Half mass radius (dimensionless)",massDistribution_%radiusHalfMass ( ), 1.0d0+sqrt(2.0d0) ,absTol=1.0d-6) + call Assert("Mass within scale radius (dimensionless)",massDistribution_%massEnclosedBySphere(1.0d0 ), 0.25d0 ,absTol=1.0d-6) + call Assert("Potential at scale radius (dimensionless)",massDistribution_%potential (position),-0.50d0 ,absTol=1.0d-6) + end select + deallocate(massDistribution_) + allocate(massDistributionHernquist :: massDistribution_) + select type (massDistribution_) + type is (massDistributionHernquist) + massDistribution_=massDistributionHernquist(dimensionless=.false.,mass=2.0d0,scaleLength=2.0d0) + end select + select type (massDistribution_) + class is (massDistributionSpherical) + position=[2.0d0,0.0d0,0.0d0] + call Assert("Half mass radius (dimensionful )",massDistribution_%radiusHalfMass ( ),2.0d0*( 1.0d0+sqrt(2.0d0)),absTol=1.0d-6) + call Assert("Mass within scale radius (dimensionful )",massDistribution_%massEnclosedBySphere(2.0d0 ),2.0d0*( 0.25d0 ),absTol=1.0d-6) end select deallocate(massDistribution_) call Unit_Tests_End_Group()