Skip to content

Commit 9bed21c

Browse files
jeffrileyjeffriley
authored
Merge pull request #731 from reinhold-willcox/nan_radius_fix
Nan radius fix
2 parents 869ad48 + 301ecaa commit 9bed21c

17 files changed

+65
-60
lines changed

‎src/BaseBinaryStar.cpp‎

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -1018,7 +1018,7 @@ double BaseBinaryStar::CalculateTimeToCoalescence(const double p_SemiMajorAxis,
10181018
double f = 1.0 - (e0 * e0);
10191019
double f_3 = f * f * f;
10201020

1021-
tC = f <= 0.0 ? 0.0 : tC * (1.0 + 0.27 * e0_10 + 0.33 * e0_20 + 0.2 * e0_100) * f_3 * sqrt(f); // check f <= 0.0 just in case a rounding error hurts us
1021+
tC = f <= 0.0 ? 0.0 : tC * (1.0 + 0.27 * e0_10 + 0.33 * e0_20 + 0.2 * e0_100) * f_3 * std::sqrt(f); // check f <= 0.0 just in case a rounding error hurts us
10221022
}
10231023

10241024
return tC;
@@ -1156,7 +1156,7 @@ bool BaseBinaryStar::ResolveSupernova() {
11561156
// Pre-SN parameters
11571157
double semiMajorAxisPrev_km = m_SemiMajorAxis * AU_TO_KM; // km - Semi-Major axis
11581158
double eccentricityPrev = m_Eccentricity; // -- - Eccentricity, written with a prev to distinguish from later use
1159-
double sqrt1MinusEccPrevSquared = sqrt(1 - eccentricityPrev * eccentricityPrev); // useful function of eccentricity
1159+
double sqrt1MinusEccPrevSquared = std::sqrt(1 - eccentricityPrev * eccentricityPrev); // useful function of eccentricity
11601160

11611161
double m1Prev = m_Supernova->SN_TotalMassAtCOFormation(); // Mo - SN star pre-SN mass
11621162
double m2Prev = m_Companion->Mass(); // Mo - CP star pre-SN mass
@@ -1168,7 +1168,7 @@ bool BaseBinaryStar::ResolveSupernova() {
11681168
double sinEccAnomaly = sin(m_Supernova->SN_EccentricAnomaly());
11691169

11701170
// Derived quantities
1171-
double omega = sqrt(G_SN*totalMassPrev / (semiMajorAxisPrev_km * semiMajorAxisPrev_km*semiMajorAxisPrev_km)); // rad/s - Keplerian orbital frequency
1171+
double omega = std::sqrt(G_SN*totalMassPrev / (semiMajorAxisPrev_km * semiMajorAxisPrev_km*semiMajorAxisPrev_km)); // rad/s - Keplerian orbital frequency
11721172

11731173
Vector3d separationVectorPrev = Vector3d( semiMajorAxisPrev_km * (cosEccAnomaly - eccentricityPrev),
11741174
semiMajorAxisPrev_km * (sinEccAnomaly) * sqrt1MinusEccPrevSquared,
@@ -1249,10 +1249,10 @@ bool BaseBinaryStar::ResolveSupernova() {
12491249
m_Unbound = true;
12501250

12511251
// Calculate the asymptotic Center of Mass velocity
1252-
double relativeVelocityAtInfinity = (G_SN*totalMass/orbitalAngularMomentum) * sqrt(eccSquared - 1);
1252+
double relativeVelocityAtInfinity = (G_SN*totalMass/orbitalAngularMomentum) * std::sqrt(eccSquared - 1);
12531253
Vector3d relativeVelocityVectorAtInfinity = relativeVelocityAtInfinity
12541254
* (-1 * (eccentricityVector.hat / m_Eccentricity)
1255-
+ sqrt(1 - 1.0 / eccSquared) * cross(orbitalAngularMomentumVector.hat, eccentricityVector.hat));
1255+
+ std::sqrt(1 - 1.0 / eccSquared) * cross(orbitalAngularMomentumVector.hat, eccentricityVector.hat));
12561256

12571257
// Calculate the asymptotic velocities of Star1 (SN) and Star2 (CP)
12581258
Vector3d component1VelocityVectorAtInfinity = (m2 / totalMass) * relativeVelocityVectorAtInfinity + centerOfMassVelocity;
@@ -1659,7 +1659,7 @@ double BaseBinaryStar::CalculateMassTransferOrbit(const double p
16591659
double massD = p_DonorMass; // donor mass
16601660
double massAtimesMassD = massA * massD; // accretor mass * donor mass
16611661
double massAplusMassD = massA + massD; // accretor mass + donor mass
1662-
double jOrb = (massAtimesMassD / massAplusMassD) * sqrt(semiMajorAxis * G1 * massAplusMassD); // orbital angular momentum
1662+
double jOrb = (massAtimesMassD / massAplusMassD) * std::sqrt(semiMajorAxis * G1 * massAplusMassD); // orbital angular momentum
16631663
double jLoss; // specific angular momentum carried away by non-conservative mass transfer
16641664

16651665
int numberIterations = fmax( floor (fabs(p_DeltaMassDonor/(MAXIMUM_MASS_TRANSFER_FRACTION_PER_STEP*massD))), 1); // number of iterations
@@ -2079,7 +2079,7 @@ double BaseBinaryStar::CalculateAngularMomentum(const double p_SemiMajorAxis,
20792079

20802080
double Is1 = ks1 * m1 * R1 * R1;
20812081
double Is2 = ks2 * m2 * R2 * R2;
2082-
double Jorb = ((m1 * m2) / (m1 + m2)) * sqrt(G1 * (m1 + m2) * p_SemiMajorAxis * (1.0 - (p_Eccentricity * p_Eccentricity)));
2082+
double Jorb = ((m1 * m2) / (m1 + m2)) * std::sqrt(G1 * (m1 + m2) * p_SemiMajorAxis * (1.0 - (p_Eccentricity * p_Eccentricity)));
20832083

20842084
return (Is1 * w1) + (Is2 * w2) + Jorb;
20852085
}

‎src/BaseBinaryStar.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -214,7 +214,7 @@ class BaseBinaryStar {
214214
MT_TRACKING MassTransferTrackerHistory() const { return m_MassTransferTrackerHistory; }
215215
bool MergesInHubbleTime() const { return m_Flags.mergesInHubbleTime; }
216216
bool OptimisticCommonEnvelope() const { return m_CEDetails.optimisticCE; }
217-
double OrbitalAngularVelocity() const { return sqrt(G1 * (m_Star1->Mass() + m_Star2->Mass()) / (m_SemiMajorAxis * m_SemiMajorAxis * m_SemiMajorAxis)); } // rads/year
217+
double OrbitalAngularVelocity() const { return std::sqrt(G1 * (m_Star1->Mass() + m_Star2->Mass()) / (m_SemiMajorAxis * m_SemiMajorAxis * m_SemiMajorAxis)); } // rads/year
218218
double OrbitalVelocityPreSN() const { return m_OrbitalVelocityPreSN; }
219219
double Periastron() const { return m_SemiMajorAxis * (1.0-m_Eccentricity); }
220220
double PeriastronRsol() const { return Periastron() * AU_TO_RSOL; }
@@ -444,7 +444,7 @@ class BaseBinaryStar {
444444

445445
double CalculateOrbitalAngularMomentum(const double p_Mu,
446446
const double p_Mass,
447-
const double p_SemiMajorAxis) const { return p_Mu * sqrt(G1 * p_Mass * p_SemiMajorAxis); }
447+
const double p_SemiMajorAxis) const { return p_Mu * std::sqrt(G1 * p_Mass * p_SemiMajorAxis); }
448448

449449
double CalculateOrbitalEnergy(const double p_Mu,
450450
const double p_Mass,

‎src/BaseStar.cpp‎

Lines changed: 12 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -1212,7 +1212,7 @@ double BaseStar::CalculateLuminosityAtZAMS(const double p_MZAMS) {
12121212

12131213
// pow() is slow - use multiplication where it makes sense
12141214
// sqrt() is much faster than pow()
1215-
double m_0_5 = sqrt(p_MZAMS);
1215+
double m_0_5 = std::sqrt(p_MZAMS);
12161216
double m_2 = p_MZAMS * p_MZAMS;
12171217
double m_3 = m_2 * p_MZAMS;
12181218
double m_5 = m_3 * m_2;
@@ -1295,8 +1295,8 @@ double BaseStar::CalculateRadiusAtZAMS(const double p_MZAMS) const {
12951295
#define coeff(x) m_RCoefficients[static_cast<int>(R_Coeff::x)] // for convenience and readability - undefined at end of function
12961296

12971297
// pow() is slow - use multiplication where it makes sense
1298-
// sqrt() is much faster than pow()
1299-
double m_0_5 = sqrt(p_MZAMS);
1298+
// std::sqrt() is much faster than pow()
1299+
double m_0_5 = std::sqrt(p_MZAMS);
13001300
double m_2 = p_MZAMS * p_MZAMS;
13011301
double m_2_5 = m_2 * m_0_5;
13021302
double m_6 = m_2 * m_2 * m_2;
@@ -1458,7 +1458,7 @@ double BaseStar::CalculateMassLossRateNieuwenhuijzenDeJager() const {
14581458
double rate = 0.0;
14591459
if (utils::Compare(m_Luminosity, NJ_MINIMUM_LUMINOSITY) > 0) { // check for minimum luminosity
14601460
double smoothTaper = min(1.0, (m_Luminosity - 4000.0) / 500.0); // Smooth taper between no mass loss and mass loss
1461-
rate = sqrt((m_Metallicity / ZSOL)) * smoothTaper * 9.6E-15 * PPOW(m_Radius, 0.81) * PPOW(m_Luminosity, 1.24) * PPOW(m_Mass, 0.16);
1461+
rate = std::sqrt((m_Metallicity / ZSOL)) * smoothTaper * 9.6E-15 * PPOW(m_Radius, 0.81) * PPOW(m_Luminosity, 1.24) * PPOW(m_Mass, 0.16);
14621462
} else {
14631463
rate = 0.0;
14641464
}
@@ -1477,7 +1477,7 @@ double BaseStar::CalculateMassLossRateNieuwenhuijzenDeJager() const {
14771477
*/
14781478
double BaseStar::CalculateMassLossRateLBV(const LBV_PRESCRIPTION p_LBV_prescription) {
14791479
double rate = 0.0;
1480-
double HD_limit_factor = m_Radius * sqrt(m_Luminosity) * 1.0E-5; // calculate factor by which you are above the HD limit
1480+
double HD_limit_factor = m_Radius * std::sqrt(m_Luminosity) * 1.0E-5; // calculate factor by which you are above the HD limit
14811481
if ((utils::Compare(m_Luminosity, LBV_LUMINOSITY_LIMIT_STARTRACK) > 0) && (utils::Compare(HD_limit_factor, 1.0) > 0)) { // check if luminous blue variable
14821482
m_LBVphaseFlag = true; // mark the star as LBV
14831483
m_DominantMassLossRate = MASS_LOSS_TYPE::LUMINOUS_BLUE_VARIABLE;
@@ -1959,7 +1959,7 @@ double BaseStar::CalculateThermalMassAcceptanceRate(const double p_Radius) const
19591959
* @return Effective temperature of the star (Tsol)
19601960
*/
19611961
double BaseStar::CalculateTemperatureOnPhase_Static(const double p_Luminosity, const double p_Radius) {
1962-
return sqrt(sqrt(p_Luminosity)) / sqrt(p_Radius); // sqrt() is much faster than pow()
1962+
return std::sqrt(std::sqrt(p_Luminosity)) / std::sqrt(p_Radius); // sqrt() is much faster than pow()
19631963
}
19641964

19651965

@@ -2188,7 +2188,7 @@ double BaseStar::CalculateZAMSAngularFrequency(const double p_MZAMS, const doubl
21882188
double BaseStar::CalculateOmegaBreak() const {
21892189
constexpr double RSOL_TO_AU_3 = RSOL_TO_AU * RSOL_TO_AU * RSOL_TO_AU;
21902190

2191-
return _2_PI * sqrt(m_Mass / (RSOL_TO_AU_3 * m_Radius * m_Radius * m_Radius));
2191+
return _2_PI * std::sqrt(m_Mass / (RSOL_TO_AU_3 * m_Radius * m_Radius * m_Radius));
21922192
}
21932193

21942194

@@ -2255,7 +2255,7 @@ double BaseStar::CalculateLifetimeToBGB(const double p_Mass) const {
22552255
// pow() is slow - use multiplication (sqrt() is much faster than pow())
22562256
double m_2 = p_Mass * p_Mass;
22572257
double m_4 = m_2 * m_2;
2258-
double m_5_5 = m_4 * p_Mass * sqrt(p_Mass);
2258+
double m_5_5 = m_4 * p_Mass * std::sqrt(p_Mass);
22592259
double m_7 = m_4 * m_2 * p_Mass;
22602260

22612261
return (a[1] + (a[2] * m_4) + (a[3] * m_5_5) + m_7) / ((a[4] * m_2) + (a[5] * m_7));
@@ -2292,7 +2292,7 @@ double BaseStar::CalculateLifetimeToBAGB(const double p_tHeI, const double p_tHe
22922292
* @return Dynamical timescale in Myr
22932293
*/
22942294
double BaseStar::CalculateDynamicalTimescale_Static(const double p_Mass, const double p_Radius) {
2295-
return 5.0 * 1.0E-5 * p_Radius * sqrt(p_Radius) * YEAR_TO_MYR / sqrt(p_Mass); // sqrt() is much faster than pow()
2295+
return 5.0 * 1.0E-5 * p_Radius * std::sqrt(p_Radius) * YEAR_TO_MYR / std::sqrt(p_Mass); // sqrt() is much faster than pow()
22962296
}
22972297

22982298

@@ -2378,7 +2378,7 @@ double BaseStar::CalculateEddyTurnoverTimescale() {
23782378

23792379
double rEnv = CalculateRadialExtentConvectiveEnvelope();
23802380

2381-
return 0.4311 * PPOW((m_Mass * rEnv * (m_Radius - (0.5 * rEnv))) / (3.0 * m_Luminosity), 1.0 / 3.0);
2381+
return 0.4311 * cbrt((m_Mass * rEnv * (m_Radius - (0.5 * rEnv))) / (3.0 * m_Luminosity));
23822382
}
23832383

23842384

@@ -2451,7 +2451,7 @@ double BaseStar::CalculateEddyTurnoverTimescale() {
24512451
* @return Drawn kick magnitude (km s^-1)
24522452
*/
24532453
double BaseStar::DrawKickMagnitudeDistributionMaxwell(const double p_Sigma, const double p_Rand) const {
2454-
return p_Sigma * sqrt(gsl_cdf_chisq_Pinv(p_Rand, 3)); // a Maxwellian is a chi distribution with three degrees of freedom
2454+
return p_Sigma * std::sqrt(gsl_cdf_chisq_Pinv(p_Rand, 3)); // a Maxwellian is a chi distribution with three degrees of freedom
24552455
}
24562456

24572457

@@ -2615,7 +2615,7 @@ double BaseStar::DrawSNKickMagnitude(const double p_Sigma,
26152615

26162616
case KICK_MAGNITUDE_DISTRIBUTION::MULLER2016MAXWELLIAN: { // MULLER2016-MAXWELLIAN
26172617

2618-
double mullerSigma = DrawRemnantKickMuller(p_COCoreMass) / sqrt(3.0);
2618+
double mullerSigma = DrawRemnantKickMuller(p_COCoreMass) / std::sqrt(3.0);
26192619

26202620
kickMagnitude = DrawKickMagnitudeDistributionMaxwell(mullerSigma, p_Rand);
26212621
} break;

‎src/BinaryConstituentStar.cpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -334,10 +334,10 @@ double BinaryConstituentStar::CalculateCircularisationTimescale(const double p_S
334334
double rInAU = Radius() * RSOL_TO_AU;
335335
double rInAUPow3 = rInAU * rInAU * rInAU; // use multiplication - pow() is slow
336336
double rOverAPow10 = rOverA * rOverA * rOverA * rOverA * rOverA * rOverA * rOverA * rOverA * rOverA * rOverA; // use multiplication - pow() is slow
337-
double rOverAPow21Over2 = rOverAPow10 * rOverA * sqrt(rOverA); // srqt() is faster than pow()
337+
double rOverAPow21Over2 = rOverAPow10 * rOverA * std::sqrt(rOverA); // sqrt() is faster than pow()
338338

339339
double secondOrderTidalCoeff = 1.592E-09 * PPOW(Mass(), 2.84); // aka E_2.
340-
double freeFallFactor = sqrt(G1 * Mass() / rInAUPow3);
340+
double freeFallFactor = std::sqrt(G1 * Mass() / rInAUPow3);
341341

342342
timescale = 1.0 / ((21.0 / 2.0) * freeFallFactor * q2 * PPOW(1.0 + q2, 11.0/6.0) * secondOrderTidalCoeff * rOverAPow21Over2);
343343
} break;
@@ -387,7 +387,7 @@ double BinaryConstituentStar::CalculateSynchronisationTimescale(const double p_S
387387
double e2 = 1.592E-9 * PPOW(Mass(), 2.84); // second order tidal coefficient (a.k.a. E_2)
388388
double rAU = Radius() * RSOL_TO_AU;
389389
double rAU_3 = rAU * rAU * rAU;
390-
double freeFallFactor = sqrt(G1 * Mass() / rAU_3);
390+
double freeFallFactor = std::sqrt(G1 * Mass() / rAU_3);
391391

392392
timescale = 1.0 / (coeff2 * freeFallFactor * gyrationRadiusSquared_1 * q2 * q2 * PPOW(1.0 + q2, 5.0 / 6.0) * e2 * PPOW(rOverA, 17.0 / 2.0));
393393
} break;

‎src/CHeB.cpp‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -94,13 +94,13 @@ double CHeB::CalculateLambdaDewi() const {
9494

9595
double lambda3 = std::min(-0.9, 0.58 + (0.75 * log10(m_Mass))) - (0.08 * log10(m_Luminosity)); // (A.4) Claeys+2014
9696
double lambda1 = std::min(lambda3, std::min(0.8, 1.25 - (0.15 * log10(m_Luminosity)))); // (A.5) Top, Claeys+2014
97-
double lambda2 = 0.42 * PPOW(m_RZAMS / m_Radius, 0.4); // (A.2) Claeys+2014
97+
double lambda2 = 0.42 * PPOW(m_RZAMS / m_Radius, 0.4); // (A.2) Claeys+2014 // RTW - Consider replacing this with a 2/5 root function (somehow) to avoid NaNs if the base is negative
9898
double envMass = utils::Compare(m_CoreMass, 0.0) > 0 && utils::Compare(m_Mass, m_CoreMass) > 0 ? m_Mass - m_CoreMass : 0.0;
9999

100100
double lambdaCE;
101101

102102
if (utils::Compare(envMass, 1.0) >= 0) lambdaCE = 2.0 * lambda1; // (A.1) Bottom, Claeys+2014
103-
else if (utils::Compare(envMass, 0.0) > 0) lambdaCE = 2.0 * (lambda2 + (sqrt(envMass) * (lambda1 - lambda2))); // (A.1) Mid, Claeys+2014
103+
else if (utils::Compare(envMass, 0.0) > 0) lambdaCE = 2.0 * (lambda2 + (std::sqrt(envMass) * (lambda1 - lambda2))); // (A.1) Mid, Claeys+2014
104104
else lambdaCE = 2.0 * lambda2; // (A.1) Top, Claeys+2014
105105

106106
return lambdaCE;
@@ -778,9 +778,9 @@ double CHeB::CalculateRadiusRho(const double p_Mass, const double p_Tau) const {
778778

779779
double ty_tx = ty - tx;
780780

781-
double one = PPOW(log((Ry / Rmin)), 1.0 / 3.0);
781+
double one = std::cbrt(log((Ry / Rmin)));
782782
double two = (p_Tau - tx) / ty_tx;
783-
double three = PPOW(log((Rx / Rmin)), 1.0 / 3.0);
783+
double three = std::cbrt(log((Rx / Rmin)));
784784
double four = (ty - p_Tau) / ty_tx;
785785

786786
return (one * two) - (three * four);

‎src/GiantBranch.cpp‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -87,7 +87,7 @@ void GiantBranch::CalculateTimescales(const double p_Mass, DBL_VECTOR &p_Timesca
8787
* @return Core mass - Luminosity relation parameter B
8888
*/
8989
double GiantBranch::CalculateCoreMass_Luminosity_B_Static(const double p_Mass) {
90-
return std::max(3.0E4, (500.0 + (1.75E4 * PPOW(p_Mass, 0.6))));
90+
return std::max(3.0E4, (500.0 + (1.75E4 * PPOW(p_Mass, 0.6)))); // RTW - Consider replacing this with a 3/5 root function (somehow) to avoid NaNs if the base is negative
9191
}
9292

9393

@@ -720,7 +720,7 @@ double GiantBranch::CalculateRadialExtentConvectiveEnvelope() const{
720720
double GiantBranch::CalculateCoreMassAtBAGB(const double p_Mass) const {
721721
#define b m_BnCoefficients // for convenience and readability - undefined at end of function
722722

723-
return sqrt(sqrt((b[36] * PPOW(p_Mass, b[37])) + b[38])); // sqrt() is much faster than PPOW()
723+
return std::sqrt(std::sqrt((b[36] * PPOW(p_Mass, b[37])) + b[38])); // sqrt() is much faster than PPOW()
724724

725725
#undef b
726726
}
@@ -743,7 +743,7 @@ double GiantBranch::CalculateCoreMassAtBAGB(const double p_Mass) const {
743743
double GiantBranch::CalculateCoreMassAtBAGB_Static(const double p_Mass, const DBL_VECTOR &p_BnCoefficients) {
744744
#define b p_BnCoefficients // for convenience and readability - undefined at end of function
745745

746-
return sqrt(sqrt((b[36] * PPOW(p_Mass, b[37])) + b[38])); // sqrt() is much faster than PPOW()
746+
return std::sqrt(std::sqrt((b[36] * PPOW(p_Mass, b[37])) + b[38])); // sqrt() is much faster than PPOW()
747747

748748
#undef b
749749
}
@@ -771,7 +771,7 @@ double GiantBranch::CalculateCoreMassAtBGB(const double p_Mass, const DBL_VECTOR
771771
double Mc_MHeF = BaseStar::CalculateCoreMassGivenLuminosity_Static(luminosity, p_GBParams);
772772
double c = (Mc_MHeF * Mc_MHeF * Mc_MHeF * Mc_MHeF) - (MC_L_C1 * PPOW(massCutoffs(MHeF), MC_L_C2)); // pow() is slow - use multiplication
773773

774-
return std::min((0.95 * gbParams(McBAGB)), sqrt(sqrt(c + (MC_L_C1 * PPOW(p_Mass, MC_L_C2))))); // sqrt is much faster than PPOW()
774+
return std::min((0.95 * gbParams(McBAGB)), std::sqrt(std::sqrt(c + (MC_L_C1 * PPOW(p_Mass, MC_L_C2))))); // sqrt is much faster than PPOW()
775775

776776
#undef massCutoffs
777777
#undef gbParams
@@ -807,7 +807,7 @@ double GiantBranch::CalculateCoreMassAtBGB_Static(const double p_Mass,
807807
double Mc_MHeF = BaseStar::CalculateCoreMassGivenLuminosity_Static(luminosity, p_GBParams);
808808
double c = (Mc_MHeF * Mc_MHeF * Mc_MHeF * Mc_MHeF) - (MC_L_C1 * PPOW(massCutoffs(MHeF), MC_L_C2)); // pow() is slow - use multiplication
809809

810-
return std::min((0.95 * gbParams(McBAGB)), sqrt(sqrt(c + (MC_L_C1 * PPOW(p_Mass, MC_L_C2))))); // sqrt is much faster than PPOW()
810+
return std::min((0.95 * gbParams(McBAGB)), std::sqrt(std::sqrt(c + (MC_L_C1 * PPOW(p_Mass, MC_L_C2))))); // sqrt is much faster than PPOW()
811811

812812
#undef massCutoffs
813813
#undef gbParams
@@ -860,7 +860,7 @@ double GiantBranch::CalculateCoreMassAtHeIgnition(const double p_Mass) const {
860860
double McBAGB = CalculateCoreMassAtBAGB(p_Mass);
861861
double c = (Mc_MHeF * Mc_MHeF * Mc_MHeF * Mc_MHeF) - (MC_L_C1 * PPOW(massCutoffs(MHeF), MC_L_C2)); // pow() is slow - use multiplication
862862

863-
coreMass = std::min((0.95 * McBAGB), sqrt(sqrt(c + (MC_L_C1 * PPOW(p_Mass, MC_L_C2))))); // sqrt() is much faster than PPOW()
863+
coreMass = std::min((0.95 * McBAGB), std::sqrt(std::sqrt(c + (MC_L_C1 * PPOW(p_Mass, MC_L_C2))))); // sqrt() is much faster than PPOW()
864864
}
865865

866866
return coreMass;

0 commit comments

Comments
 (0)