diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index c590e2883..3675a21c0 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -1220,7 +1220,7 @@ bool BaseBinaryStar::ResolveSupernova() { double aPrev_2 = aPrev * aPrev; double aPrev_3 = aPrev_2 * aPrev; - double omega = std::sqrt(G_SN * totalMassPrev / aPrev_3); // rad/s - Keplerian orbital frequency + double omega = std::sqrt(G_km_Msol_s * totalMassPrev / aPrev_3); // rad/s - Keplerian orbital frequency Vector3d separationVectorPrev = Vector3d(aPrev * (cosEccAnomaly - eccentricityPrev), aPrev * (sinEccAnomaly) * sqrt1MinusEccPrevSquared, @@ -1234,7 +1234,8 @@ bool BaseBinaryStar::ResolveSupernova() { Vector3d orbitalAngularMomentumVectorPrev = cross(separationVectorPrev, relativeVelocityVectorPrev); // km^2 s^-1 - Specific orbital angular momentum vector - Vector3d eccentricityVectorPrev = cross(relativeVelocityVectorPrev, orbitalAngularMomentumVectorPrev) / (G_SN * totalMassPrev) - separationVectorPrev.hat; // -- - Laplace-Runge-Lenz vector (magnitude = eccentricity) + Vector3d eccentricityVectorPrev = cross(relativeVelocityVectorPrev, orbitalAngularMomentumVectorPrev) / + (G_km_Msol_s * totalMassPrev) - separationVectorPrev.hat; // -- - Laplace-Runge-Lenz vector (magnitude = eccentricity) m_OrbitalVelocityPreSN = relativeVelocityVectorPrev.mag; // km/s - Set the Pre-SN orbital velocity and m_uK = m_Supernova->SN_KickMagnitude() / m_OrbitalVelocityPreSN; // -- - Dimensionless kick magnitude @@ -1268,11 +1269,12 @@ bool BaseBinaryStar::ResolveSupernova() { m_OrbitalAngularMomentumVector = orbitalAngularMomentumVector / orbitalAngularMomentum; // set unit vector here to make printing out the inclination vector easier Vector3d eccentricityVector = cross(relativeVelocityVector, orbitalAngularMomentumVector) / - (G_SN * totalMass) - separationVectorPrev / separationPrev; // PostSN Laplace-Runge-Lenz vector + (G_km_Msol_s * totalMass) - separationVectorPrev / separationPrev; // PostSN Laplace-Runge-Lenz vector m_Eccentricity = eccentricityVector.mag; // PostSN eccentricity double eccSquared = m_Eccentricity * m_Eccentricity; // useful function of eccentricity - double semiMajorAxis_km = (orbitalAngularMomentum * orbitalAngularMomentum) / (G_SN * totalMass * (1.0 - eccSquared)); // km - PostSN semi-major axis + double semiMajorAxis_km = (orbitalAngularMomentum * orbitalAngularMomentum) / + (G_km_Msol_s * totalMass * (1.0 - eccSquared)); // km - PostSN semi-major axis m_SemiMajorAxis = semiMajorAxis_km * KM_TO_AU; // AU - PostSN semi-major axis // Note: similar to above, @@ -1288,7 +1290,7 @@ bool BaseBinaryStar::ResolveSupernova() { m_Unbound = true; // Calculate the asymptotic Center of Mass velocity - double relativeVelocityAtInfinity = (G_SN*totalMass/orbitalAngularMomentum) * std::sqrt(eccSquared - 1.0); + double relativeVelocityAtInfinity = (G_km_Msol_s*totalMass/orbitalAngularMomentum) * std::sqrt(eccSquared - 1.0); Vector3d relativeVelocityVectorAtInfinity = relativeVelocityAtInfinity * (-1.0 * (eccentricityVector.hat / m_Eccentricity) + std::sqrt(1.0 - 1.0 / eccSquared) * cross(orbitalAngularMomentumVector.hat, eccentricityVector.hat)); @@ -1722,7 +1724,7 @@ double BaseBinaryStar::CalculateMassTransferOrbit(const double p double massD = p_DonorMass; // donor mass double massAtimesMassD = massA * massD; // accretor mass * donor mass double massAplusMassD = massA + massD; // accretor mass + donor mass - double jOrb = (massAtimesMassD / massAplusMassD) * std::sqrt(semiMajorAxis * G1 * massAplusMassD); // orbital angular momentum + double jOrb = (massAtimesMassD / massAplusMassD) * std::sqrt(semiMajorAxis * G_AU_Msol_yr * massAplusMassD); // orbital angular momentum double jLoss; // specific angular momentum carried away by non-conservative mass transfer if (utils::Compare(p_DeltaMassDonor, 0.0) < 0) { // mass loss from donor? @@ -2088,7 +2090,7 @@ double BaseBinaryStar::CalculateTotalEnergy(const double p_SemiMajorAxis, double Is1 = ks1 * m1 * R1 * R1 * RSOL_TO_AU_2; double Is2 = ks2 * m2 * R2 * R2 * RSOL_TO_AU_2; - return (0.5 * Is1 * w1 * w1) + (0.5 * Is2 * w2 * w2) - (0.5 * G1 * m1 * m2 / p_SemiMajorAxis); + return (0.5 * Is1 * w1 * w1) + (0.5 * Is2 * w2 * w2) - (0.5 * G_AU_Msol_yr * m1 * m2 / p_SemiMajorAxis); } @@ -2329,7 +2331,8 @@ void BaseBinaryStar::EvaluateBinary(const double p_Dt) { m_Star1->SetOmega(m_Omega); // synchronise star 1 m_Star2->SetOmega(m_Omega); // synchronise star 2 - m_SemiMajorAxis = std::cbrt(G1) * std::cbrt((m_Star1->Mass() + m_Star2->Mass())) / PPOW(m_Omega, 2.0 / 3.0); // re-calculate semi-major axis + m_SemiMajorAxis = std::cbrt(G_AU_Msol_yr) * std::cbrt((m_Star1->Mass() + m_Star2->Mass())) / + PPOW(m_Omega, 2.0 / 3.0); // re-calculate semi-major axis m_Eccentricity = 0.0; // circularise m_TotalAngularMomentum = CalculateAngularMomentum(); // re-calculate total angular momentum diff --git a/src/BaseBinaryStar.h b/src/BaseBinaryStar.h index 4868d929d..0386bfebf 100644 --- a/src/BaseBinaryStar.h +++ b/src/BaseBinaryStar.h @@ -220,7 +220,7 @@ class BaseBinaryStar { bool MergesInHubbleTime() const { return m_Flags.mergesInHubbleTime; } double Omega() const { return m_Omega; } bool OptimisticCommonEnvelope() const { return m_CEDetails.optimisticCE; } - double OrbitalAngularVelocity() const { return std::sqrt(G1 * (m_Star1->Mass() + m_Star2->Mass()) / (m_SemiMajorAxis * m_SemiMajorAxis * m_SemiMajorAxis)); } // rads/year + double OrbitalAngularVelocity() const { return std::sqrt(G_AU_Msol_yr * (m_Star1->Mass() + m_Star2->Mass()) / (m_SemiMajorAxis * m_SemiMajorAxis * m_SemiMajorAxis)); } // rads/year double OrbitalVelocityPreSN() const { return m_OrbitalVelocityPreSN; } double Periastron() const { return m_SemiMajorAxis * (1.0 - m_Eccentricity); } double PeriastronRsol() const { return Periastron() * AU_TO_RSOL; } @@ -456,11 +456,11 @@ class BaseBinaryStar { double CalculateOrbitalAngularMomentum(const double p_Star1Mass, const double p_Star2Mass, const double p_SemiMajorAxis, - const double p_Eccentricity) const { return ((p_Star1Mass * p_Star2Mass) / (p_Star1Mass + p_Star2Mass)) * std::sqrt(G1 * (p_Star1Mass + p_Star2Mass) * p_SemiMajorAxis * (1.0 - (p_Eccentricity * p_Eccentricity))); } + const double p_Eccentricity) const { return ((p_Star1Mass * p_Star2Mass) / (p_Star1Mass + p_Star2Mass)) * std::sqrt(G_AU_Msol_yr * (p_Star1Mass + p_Star2Mass) * p_SemiMajorAxis * (1.0 - (p_Eccentricity * p_Eccentricity))); } double CalculateOrbitalEnergy(const double p_Mu, const double p_Mass, - const double p_SemiMajorAxis) const { return -(G1 * p_Mu * p_Mass) / (2.0 * p_SemiMajorAxis); } + const double p_SemiMajorAxis) const { return -(G_AU_Msol_yr * p_Mu * p_Mass) / (2.0 * p_SemiMajorAxis); } double CalculateZetaRocheLobe(const double p_jLoss) const; @@ -662,7 +662,7 @@ class BaseBinaryStar { // function: (I_1 + I_2) Omega + L(Omega) - p_Ltot = 0 // where L(Omega) = b*Omega(-1/3) double a = p_I1 + p_I2; // I_1 + I_2 - double b = PPOW(G1, 2.0 / 3.0) * p_M1 * p_M2 / std::cbrt(p_M1 + p_M2); + double b = PPOW(G_AU_Msol_yr, 2.0 / 3.0) * p_M1 * p_M2 / std::cbrt(p_M1 + p_M2); double c = -p_Ltot; auto func = [a, b, c](double x) -> double { return (a * x) + (b / std::cbrt(x)) + c; }; // functor diff --git a/src/BinaryConstituentStar.cpp b/src/BinaryConstituentStar.cpp index 154e18cd3..613fb4285 100644 --- a/src/BinaryConstituentStar.cpp +++ b/src/BinaryConstituentStar.cpp @@ -358,7 +358,7 @@ double BinaryConstituentStar::CalculateCircularisationTimescale(const double p_S double rOverAPow21Over2 = rOverAPow10 * rOverA * std::sqrt(rOverA); // sqrt() is faster than pow() double secondOrderTidalCoeff = 1.592E-09 * PPOW(Mass(), 2.84); // aka E_2. - double freeFallFactor = std::sqrt(G1 * Mass() / rInAUPow3); + double freeFallFactor = std::sqrt(G_AU_Msol_yr * Mass() / rInAUPow3); timescale = 1.0 / ((21.0 / 2.0) * freeFallFactor * q2 * PPOW(1.0 + q2, 11.0/6.0) * secondOrderTidalCoeff * rOverAPow21Over2); } break; @@ -408,7 +408,7 @@ double BinaryConstituentStar::CalculateSynchronisationTimescale(const double p_S double e2 = 1.592E-9 * PPOW(Mass(), 2.84); // second order tidal coefficient (a.k.a. E_2) double rAU = Radius() * RSOL_TO_AU; double rAU_3 = rAU * rAU * rAU; - double freeFallFactor = std::sqrt(G1 * Mass() / rAU_3); + double freeFallFactor = std::sqrt(G_AU_Msol_yr * Mass() / rAU_3); timescale = 1.0 / (coeff2 * freeFallFactor * gyrationRadiusSquared_1 * q2 * q2 * PPOW(1.0 + q2, 5.0 / 6.0) * e2 * PPOW(rOverA, 17.0 / 2.0)); } break; diff --git a/src/constants.h b/src/constants.h index 59bc67515..33f275ad1 100755 --- a/src/constants.h +++ b/src/constants.h @@ -247,8 +247,8 @@ constexpr double HUBBLE_TIME = 1 / H0SI; constexpr double G = 6.67E-11; // Gravitational constant in m^3 kg^-1 s^-2 (more accurately known as G M_sol) constexpr double G_CGS = 6.6743E-8; // Gravitational constant in cm^3 g^-1 s^-2 -constexpr double G1 = 4.0 * PI_2; // Gravitational constant in AU^3 Msol^-1 yr^-2 -constexpr double G_SN = G * 1.0E-9 / KG_TO_MSOL; // Gravitational constant in km^3 Msol^-1 s^-2, for use in the ResolveSupernova() function +constexpr double G_AU_Msol_yr = 4.0 * PI_2; // Gravitational constant in AU^3 Msol^-1 yr^-2 +constexpr double G_km_Msol_s = G * 1.0E-9 / KG_TO_MSOL; // Gravitational constant in km^3 Msol^-1 s^-2, for use in the ResolveSupernova() function constexpr double G_SOLAR_YEAR = 3.14E7; // Gravitational constant in Lsol Rsol yr Msol^-2 for calculating photon tiring limit constexpr double RSOL = 6.957E8; // Solar Radius (in m) diff --git a/src/utils.cpp b/src/utils.cpp index 063f9710e..47103d267 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -233,7 +233,7 @@ namespace utils { double ConvertPeriodInDaysToSemiMajorAxisInAU(const double p_Mass1, const double p_Mass2, const double p_Period) { double a_cubed_SI_top = G * ((p_Mass1 * MSOL_TO_KG) + (p_Mass2 * MSOL_TO_KG)) * p_Period * p_Period * SECONDS_IN_DAY * SECONDS_IN_DAY; - double a_cubed_SI = a_cubed_SI_top / G1; + double a_cubed_SI = a_cubed_SI_top / G_AU_Msol_yr; double a_SI = std::cbrt(a_cubed_SI); return a_SI / AU;