Fixes an issue that manifests for low metallicities. This is a continuation of the fixes for issue #978
The problem here was that Hurley et al. 2000 has this on p29:
b47 = 1.127733ρ + 0.2344416ρ2 − 0.3793726ρ3
but the sse code has this (in zcnsts.f:
b47 = 1.127733ρ + 0.2344416ρ2 − 0.3793726ρ3
(actually, in the code it is:
gbp(66) = 1.d0 - lzd(xh(77) + lzd(xh(78) + lzd*xh(79)))
)
At first I thought that was just a mistake in the paper, but b[47] is only used once in the paper (eq 58, to calculate Tbl on p12), and there the equation 58 uses (1 - b[47]), so Jarrod has just put the "1.0 -" into the constant when it is calculated.
However, in the sse code, b[47] is also used to calculate the radius on the HG phase. I recently (in the fix for issue #978) added that code to COMPAS (see HG::CalculateRadiusOnPhase()), but I wasn't aware that the constant had the "1.0 -" prefix already added...
So, to keep the code looking like the paper, COMPAS uses the equation on p29 of the paper to calculate b[47], and adds the "1.0 -" prefix in the code (in CHeB::CalculateLifetimeOnBluePhase(), and now in HG::CalculateRadiusOnPhase()).
Fixes an issue that manifests for low metallicities. This is a continuation of the fixes for issue #978
The problem here was that Hurley et al. 2000 has this on p29:
but the sse code has this (in
zcnsts.f:(actually, in the code it is: gbp(66) = 1.d0 - lzd(xh(77) + lzd(xh(78) + lzd*xh(79))) )
At first I thought that was just a mistake in the paper, but b[47] is only used once in the paper (eq 58, to calculate Tbl on p12), and there the equation 58 uses (1 - b[47]), so Jarrod has just put the "1.0 -" into the constant when it is calculated.
However, in the sse code, b[47] is also used to calculate the radius on the HG phase. I recently (in the fix for issue #978) added that code to COMPAS (see HG::CalculateRadiusOnPhase()), but I wasn't aware that the constant had the "1.0 -" prefix already added...
So, to keep the code looking like the paper, COMPAS uses the equation on p29 of the paper to calculate b[47], and adds the "1.0 -" prefix in the code (in CHeB::CalculateLifetimeOnBluePhase(), and now in HG::CalculateRadiusOnPhase()).