Skip to content

Commit

Permalink
Merge pull request #652 from galacticusorg/fixBetaProfileMoments
Browse files Browse the repository at this point in the history 
Implement missing calculation of beta-profile radial moment for radii of 0 and ∞
  • Loading branch information
@abensonca
abensonca authored Jul 26, 2024
2 parents d08a189 + 285c094 commit deaa21a
Show file tree
Hide file tree
Showing 2 changed files with 63 additions and 8 deletions.
34 changes: 29 additions & 5 deletions source/mass_distributions.spherical.beta_profile.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -463,9 +463,9 @@ double precision function betaProfilePotential(self,coordinates,componentType,ma
& )
end if
end if
if (.not.self%isDimensionless()) &
& betaProfilePotential= &
& betaProfilePotential &
if (.not.self%isDimensionless()) &
& betaProfilePotential= &
& betaProfilePotential &
& *gravitationalConstantGalacticus
return
end function betaProfilePotential
Expand All @@ -476,6 +476,8 @@ double precision function betaProfileDensityRadialMoment(self,moment,radiusMinim
!!}
use :: Hypergeometric_Functions, only : Hypergeometric_2F1
use :: Numerical_Comparison , only : Values_Agree
use :: Error , only : Error_Report
use :: Gamma_Functions , only : Gamma_Function
implicit none
class (massDistributionBetaProfile ), intent(inout) :: self
double precision , intent(in ) :: moment
Expand Down Expand Up @@ -512,7 +514,7 @@ double precision function betaProfileDensityRadialMoment(self,moment,radiusMinim
& +radialMomentTwoThirds(specialCaseMoment,fractionalRadiusMaximum)
else
! General solution.
betaProfileDensityRadialMoment= &
betaProfileDensityRadialMoment= &
& +fractionalRadiusMaximum**(moment+1.0d0) &
& *Hypergeometric_2F1 ( &
& [(moment+1.0d0)/2.0d0,1.5d0*self%beta], &
Expand All @@ -522,12 +524,34 @@ double precision function betaProfileDensityRadialMoment(self,moment,radiusMinim
& / (moment+1.0d0)
end if
else
betaProfileDensityRadialMoment=0.0d0
if (moment >= 3.0d0*self%beta-1.0d0) then
betaProfileDensityRadialMoment=0.0d0
if (present(isInfinite)) then
isInfinite=.true.
return
else
call Error_Report('radial moment is infinite at r=0 for m ≥ 3β-1'//{introspection:location})
end if
else
betaProfileDensityRadialMoment=+Gamma_Function((+3.0d0+moment )/2.0d0) &
& *Gamma_Function((-1.0d0-moment+3.0d0*self%beta)/2.0d0) &
& /Gamma_Function(( +3.0d0*self%beta)/2.0d0) &
& / (+1.0d0+moment )
end if
end if
if (present(radiusMinimum)) then
fractionalRadiusMinimum=radiusMinimum/self%coreRadius
else
fractionalRadiusMinimum=0.0d0
if (moment <= -1.0d0) then
betaProfileDensityRadialMoment=0.0d0
if (present(isInfinite)) then
isInfinite=.true.
return
else
call Error_Report('radial moment is infinite at r=0 for m ≤ -1'//{introspection:location})
end if
end if
end if
if (specialCaseMoment /= -huge(0)) then
! Special case for 0ᵗʰ, 1ˢᵗ, 2ⁿᵈ, and 3ʳᵈ moments of a β=2/3 distribution.
Expand Down
37 changes: 34 additions & 3 deletions source/tests.mass_distributions.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -137,6 +137,7 @@ program Test_Mass_Distributions

! Beta profile.
call Unit_Tests_Begin_Group("Beta profile")
call Unit_Tests_Begin_Group("β = 2/3")
allocate(massDistributionBetaProfile :: massDistribution_)
select type (massDistribution_)
type is (massDistributionBetaProfile)
Expand All @@ -160,19 +161,19 @@ program Test_Mass_Distributions
& absTol=1.0d-6 &
& )
call Assert( &
& "Radial moment, m=1" , &
& "Radial moment, m=1, from 0 to 1" , &
& +massDistribution_%densityRadialMoment (1.0d0,0.0d0,1.0d0) , &
& +0.3465735902799727d0 , &
& absTol=1.0d-6 &
& )
call Assert( &
& "Radial moment, m=2" , &
& "Radial moment, m=2, from 0 to 1" , &
& +massDistribution_%densityRadialMoment (2.0d0,0.0d0,1.0d0) , &
& +0.2146018366025517d0 , &
& absTol=1.0d-6 &
& )
call Assert( &
& "Radial moment, m=3" , &
& "Radial moment, m=3, from 0 to 1" , &
& +massDistribution_%densityRadialMoment (3.0d0,0.0d0,1.0d0) , &
& +0.1534264097200273d0 , &
& absTol=1.0d-6 &
Expand All @@ -184,6 +185,36 @@ program Test_Mass_Distributions
& absTol=1.0d-6 &
& )
end select
call Unit_Tests_End_Group()
deallocate(massDistribution_)
call Unit_Tests_Begin_Group("β = 2")
allocate(massDistributionBetaProfile :: massDistribution_)
select type (massDistribution_)
type is (massDistributionBetaProfile)
massDistribution_=massDistributionBetaProfile(beta=2.0d0,dimensionless=.true.)
end select
select type (massDistribution_)
class is (massDistributionSpherical)
call Assert( &
& "Radial moment, m=1, unbounded" , &
& +massDistribution_%densityRadialMoment (1.0d0 ) , &
& +1.0d0/4.0d0 , &
& absTol=1.0d-6 &
& )
call Assert( &
& "Radial moment, m=2, unbounded" , &
& +massDistribution_%densityRadialMoment (2.0d0 ) , &
& +Pi/16.0d0 , &
& absTol=1.0d-6 &
& )
call Assert( &
& "Radial moment, m=3, unbounded" , &
& +massDistribution_%densityRadialMoment (3.0d0 ) , &
& +1.0d0/4.0d0 , &
& absTol=1.0d-6 &
& )
end select
call Unit_Tests_End_Group()
! Ensure that a dimensionful profile produces the correct mass inside of its outer radius.
deallocate(massDistribution_)
allocate(massDistributionBetaProfile :: massDistribution_)
Expand Down

0 comments on commit deaa21a

Please sign in to comment.