fix: Correct enclosed mass in dimensionfull Hernquist profiles

galacticusorg · Sep 15, 2023 · c02ffb6 · c02ffb6
1 parent 90736f9
commit c02ffb6
Show file tree

Hide file tree

Showing 2 changed files with 17 additions and 5 deletions.
diff --git 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
@@ -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()