Skip to content

Commit

Permalink
renamed G_SN to G_km_Msol_s and G1 to G_AU_Msol_yr
Browse files Browse the repository at this point in the history
  • Loading branch information
@reinhold-willcox
reinhold-willcox committed Nov 15, 2023
1 parent 9cd502d commit 8076e3f
Show file tree
Hide file tree
Showing 5 changed files with 20 additions and 17 deletions.
19 changes: 11 additions & 8 deletions src/BaseBinaryStar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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
Expand Down Expand Up @@ -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,
Expand All @@ -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));
Expand Down Expand Up @@ -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?
Expand Down Expand Up @@ -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);
}


Expand Down Expand Up @@ -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

Expand Down
8 changes: 4 additions & 4 deletions src/BaseBinaryStar.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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; }
Expand Down Expand Up @@ -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;

Expand Down Expand Up @@ -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
Expand Down
4 changes: 2 additions & 2 deletions src/BinaryConstituentStar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand Down
4 changes: 2 additions & 2 deletions src/constants.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion src/utils.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down

0 comments on commit 8076e3f

Please sign in to comment.