From b770c99b247ff894a776401b88e57d45068ffc69 Mon Sep 17 00:00:00 2001 From: Veome Kapil Date: Wed, 10 Apr 2024 14:31:35 +1000 Subject: [PATCH] Modified CalculateImKlmTidal() to return a tuple with the appropriate moments of ImK in order --- src/BaseBinaryStar.cpp | 22 ++++++++----------- src/BaseBinaryStar.h | 43 ++++++++++++++++++------------------ src/BaseStar.h | 2 +- src/GiantBranch.cpp | 50 +++++++++++++++++++++++++++++++----------- src/GiantBranch.h | 2 +- src/Star.h | 2 +- src/constants.h | 1 + 7 files changed, 72 insertions(+), 50 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 623a8b207..730407cce 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -2342,9 +2342,10 @@ void BaseBinaryStar::EvaluateBinary(const double p_Dt) { m_Omega = OrbitalAngularVelocity(); } - - double DSemiMajorAxis1Dt = CalculateDSemiMajorAxisTidalDt(m_Star1->CalculateImKlmTidal(m_Omega, 1, 0), m_Star1->CalculateImKlmTidal(m_Omega, 1, 2), - m_Star1->CalculateImKlmTidal(m_Omega, 2, 2), m_Star1->CalculateImKlmTidal(m_Omega, 3, 2), + DBL_DBL_DBL_DBL ImKlm1 = m_Star1->CalculateImKlmTidal(m_Omega); + DBL_DBL_DBL_DBL ImKlm2 = m_Star2->CalculateImKlmTidal(m_Omega); + + double DSemiMajorAxis1Dt = CalculateDSemiMajorAxisTidalDt(ImKlm1, m_Star1->Mass(), m_Star1->Radius(), m_Star2->Mass(), @@ -2352,8 +2353,7 @@ void BaseBinaryStar::EvaluateBinary(const double p_Dt) { m_SemiMajorAxis, m_Eccentricity); - double DSemiMajorAxis2Dt = CalculateDSemiMajorAxisTidalDt(m_Star1->CalculateImKlmTidal(m_Omega, 1, 0), m_Star1->CalculateImKlmTidal(m_Omega, 1, 2), - m_Star1->CalculateImKlmTidal(m_Omega, 2, 2), m_Star1->CalculateImKlmTidal(m_Omega, 3, 2), + double DSemiMajorAxis2Dt = CalculateDSemiMajorAxisTidalDt(ImKlm2, m_Star2->Mass(), m_Star2->Radius(), m_Star1->Mass(), @@ -2361,8 +2361,7 @@ void BaseBinaryStar::EvaluateBinary(const double p_Dt) { m_SemiMajorAxis, m_Eccentricity); - double DEccentricity1Dt = CalculateDEccentricityTidalDt(m_Star1->CalculateImKlmTidal(m_Omega, 1, 0), m_Star1->CalculateImKlmTidal(m_Omega, 1, 2), - m_Star1->CalculateImKlmTidal(m_Omega, 2, 2), m_Star1->CalculateImKlmTidal(m_Omega, 3, 2), + double DEccentricity1Dt = CalculateDEccentricityTidalDt(ImKlm1, m_Star1->Mass(), m_Star1->Radius(), m_Star2->Mass(), @@ -2370,8 +2369,7 @@ void BaseBinaryStar::EvaluateBinary(const double p_Dt) { m_SemiMajorAxis, m_Eccentricity); - double DEccentricity2Dt = CalculateDEccentricityTidalDt(m_Star1->CalculateImKlmTidal(m_Omega, 1, 0), m_Star1->CalculateImKlmTidal(m_Omega, 1, 2), - m_Star1->CalculateImKlmTidal(m_Omega, 2, 2), m_Star1->CalculateImKlmTidal(m_Omega, 3, 2), + double DEccentricity2Dt = CalculateDEccentricityTidalDt(ImKlm2, m_Star2->Mass(), m_Star2->Radius(), m_Star1->Mass(), @@ -2379,8 +2377,7 @@ void BaseBinaryStar::EvaluateBinary(const double p_Dt) { m_SemiMajorAxis, m_Eccentricity); - double DOmega1Dt = CalculateDOmegaTidalDt(m_Star1->CalculateImKlmTidal(m_Omega, 1, 0), m_Star1->CalculateImKlmTidal(m_Omega, 1, 2), - m_Star1->CalculateImKlmTidal(m_Omega, 2, 2), m_Star1->CalculateImKlmTidal(m_Omega, 3, 2), + double DOmega1Dt = CalculateDOmegaTidalDt(ImKlm1, m_Star1->Mass(), m_Star1->Radius(), m_Star1->CalculateMomentOfInertiaAU(), @@ -2389,8 +2386,7 @@ void BaseBinaryStar::EvaluateBinary(const double p_Dt) { m_SemiMajorAxis, m_Eccentricity); - double DOmega2Dt = CalculateDOmegaTidalDt(m_Star1->CalculateImKlmTidal(m_Omega, 1, 0), m_Star1->CalculateImKlmTidal(m_Omega, 1, 2), - m_Star1->CalculateImKlmTidal(m_Omega, 2, 2), m_Star1->CalculateImKlmTidal(m_Omega, 3, 2), + double DOmega2Dt = CalculateDOmegaTidalDt(ImKlm2, m_Star2->Mass(), m_Star2->Radius(), m_Star2->CalculateMomentOfInertiaAU(), diff --git a/src/BaseBinaryStar.h b/src/BaseBinaryStar.h index 7ea37d3d4..281ef84db 100644 --- a/src/BaseBinaryStar.h +++ b/src/BaseBinaryStar.h @@ -774,12 +774,9 @@ class BaseBinaryStar { * Zahn, 1977, Eq. (3.6) * * - * double CalculateDSemiMajorAxisTidalDt(const double p_ImK10, const double p_ImK12, const double p_ImK22, const double p_ImK32, const double p_M1, const double p_R1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity) + * double CalculateDSemiMajorAxisTidalDt(const std::tuple p_ImKlm, const double p_M1, const double p_R1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity) * - * @param [IN] p_ImK10 Imaginary (1,0) component of potential tidal love number of star (unitless) - * @param [IN] p_ImK12 Imaginary (1,2) component of potential tidal love number of star (unitless) - * @param [IN] p_ImK22 Imaginary (2,2) component of potential tidal love number of star (unitless) - * @param [IN] p_ImK32 Imaginary (3,2) component of potential tidal love number of star (unitless) + * @param [IN] p_ImKlm Imaginary [(1,0), (1,2), (2,2), (3,2)] components of the potential tidal love number of star (unitless) * @param [IN] p_M1 Mass of star (Msol) * @param [IN] p_R1 Radius of star (Rsol) * @param [IN] p_M2 Mass of companion star (Msol) @@ -788,12 +785,16 @@ class BaseBinaryStar { * @param [IN] p_Eccentricity Eccentricity for binary * @return Change in semi-major axis for binary (AU/yr) */ - double CalculateDSemiMajorAxisTidalDt(const double p_ImK10, const double p_ImK12, const double p_ImK22, const double p_ImK32, const double p_M1, const double p_R1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity) { + double CalculateDSemiMajorAxisTidalDt(const std::tuple p_ImKlm, const double p_M1, const double p_R1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity) { + + double ImK10, ImK12, ImK22, ImK32; + std::tie(ImK10, ImK12, ImK22, ImK32) = p_ImKlm; + double R1_AU = p_R1 * RSOL_TO_AU; double R1_over_a = R1_AU / p_SemiMajorAxis; double R1_over_a_7 = R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a; - return - (3.0 / p_Omega) * (1.0 + (p_M2 / p_M1)) * (G_AU_Msol_yr * p_M2/ R1_AU / R1_AU) * R1_over_a_7 * (p_ImK22 + ((p_Eccentricity*p_Eccentricity) * ((3.0 * p_ImK10 / 4.0) + (p_ImK12 / 8.0) - (5.0 * p_ImK22) + (147.0 * p_ImK32 / 8.0)))); + return - (3.0 / p_Omega) * (1.0 + (p_M2 / p_M1)) * (G_AU_Msol_yr * p_M2/ R1_AU / R1_AU) * R1_over_a_7 * (ImK22 + ((p_Eccentricity*p_Eccentricity) * ((3.0 * ImK10 / 4.0) + (ImK12 / 8.0) - (5.0 * ImK22) + (147.0 * ImK32 / 8.0)))); } /* @@ -801,12 +802,9 @@ class BaseBinaryStar { * Zahn, 1977, Eq. (3.7) * * - * double CalculateDEccentricityTidalDt(const double p_ImK10, const double p_ImK12, const double p_ImK22, const double p_ImK32, const double p_M1, const double p_R1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity) + * double CalculateDEccentricityTidalDt(const std::tuple p_ImKlm, const double p_M1, const double p_R1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity) * - * @param [IN] p_ImK10 Imaginary (1,0) component of potential tidal love number of star (unitless) - * @param [IN] p_ImK12 Imaginary (1,2) component of potential tidal love number of star (unitless) - * @param [IN] p_ImK22 Imaginary (2,2) component of potential tidal love number of star (unitless) - * @param [IN] p_ImK32 Imaginary (3,2) component of potential tidal love number of star (unitless) + * @param [IN] p_ImKlm Imaginary [(1,0), (1,2), (2,2), (3,2)] components of the potential tidal love number of star (unitless) * @param [IN] p_M1 Mass of star (Msol) * @param [IN] p_R1 Radius of star (Rsol) * @param [IN] p_M2 Mass of companion star (Msol) @@ -815,13 +813,16 @@ class BaseBinaryStar { * @param [IN] p_Eccentricity Eccentricity for binary * @return Change in Eccentricity for binary (1/yr) */ - double CalculateDEccentricityTidalDt(const double p_ImK10, const double p_ImK12, const double p_ImK22, const double p_ImK32, const double p_M1, const double p_R1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity) { + double CalculateDEccentricityTidalDt(const std::tuple p_ImKlm, const double p_M1, const double p_R1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity) { + double ImK10, ImK12, ImK22, ImK32; + std::tie(ImK10, ImK12, ImK22, ImK32) = p_ImKlm; + double R1_AU = p_R1 * RSOL_TO_AU; double R1_over_a = R1_AU / p_SemiMajorAxis; double R1_over_a_8 = R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a; - return -(3.0 / 4.0) * (p_Eccentricity/p_Omega) * (1.0 + (p_M2 / p_M1)) * (G_AU_Msol_yr * p_M2 / R1_AU / R1_AU / R1_AU) * R1_over_a_8 * ((3.0 * p_ImK10 / 2.0) - (p_ImK12 / 4.0) - p_ImK22 + (49.0 * p_ImK32 / 4.0)); + return -(3.0 / 4.0) * (p_Eccentricity/p_Omega) * (1.0 + (p_M2 / p_M1)) * (G_AU_Msol_yr * p_M2 / R1_AU / R1_AU / R1_AU) * R1_over_a_8 * ((3.0 * ImK10 / 2.0) - (ImK12 / 4.0) - ImK22 + (49.0 * ImK32 / 4.0)); } /* @@ -829,12 +830,9 @@ class BaseBinaryStar { * Zahn, 1977, Eq. (3.8) * * - * double CalculateDOmegaTidalDt(const double p_ImK10, const double p_ImK12, const double p_ImK22, const double p_ImK32, const double p_M1, const double p_R1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity) + * double CalculateDOmegaTidalDt(const std::tuple p_ImKlm, const double p_M1, const double p_R1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity) * - * @param [IN] p_ImK10 Imaginary (1,0) component of potential tidal love number of star (unitless) - * @param [IN] p_ImK12 Imaginary (1,2) component of potential tidal love number of star (unitless) - * @param [IN] p_ImK22 Imaginary (2,2) component of potential tidal love number of star (unitless) - * @param [IN] p_ImK32 Imaginary (3,2) component of potential tidal love number of star (unitless) + * @param [IN] p_ImKlm Imaginary [(1,0), (1,2), (2,2), (3,2)] components of the potential tidal love number of star (unitless) * @param [IN] p_M1 Mass of star (Msol) * @param [IN] p_R1 Radius of star (Rsol) * @param [IN] p_I1 Moment of Inertia of star (Msol * AU^2) @@ -843,13 +841,16 @@ class BaseBinaryStar { * @param [IN] p_Eccentricity Eccentricity for binary * @return Change in Omega for star (1/yr/yr) */ - double CalculateDOmegaTidalDt(const double p_ImK10, const double p_ImK12, const double p_ImK22, const double p_ImK32, const double p_M1, const double p_R1, const double p_I1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity) { + double CalculateDOmegaTidalDt(const std::tuple p_ImKlm, const double p_M1, const double p_R1, const double p_I1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity) { + double ImK10, ImK12, ImK22, ImK32; + std::tie(ImK10, ImK12, ImK22, ImK32) = p_ImKlm; + double R1_AU = p_R1 * RSOL_TO_AU; double R1_over_a = R1_AU / p_SemiMajorAxis; double R1_over_a_6 = R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a; - return (3.0 / 2.0) * (1/p_I1) * (G_AU_Msol_yr * p_M2 * p_M2 / R1_AU) * R1_over_a_6 * (p_ImK22 + ((p_Eccentricity*p_Eccentricity) * ((p_ImK12 / 4.0) - (5.0 * p_ImK22) + (49.0 * p_ImK32 / 4.0)))); + return (3.0 / 2.0) * (1/p_I1) * (G_AU_Msol_yr * p_M2 * p_M2 / R1_AU) * R1_over_a_6 * (ImK22 + ((p_Eccentricity*p_Eccentricity) * ((ImK12 / 4.0) - (5.0 * ImK22) + (49.0 * ImK32 / 4.0)))); } }; diff --git a/src/BaseStar.h b/src/BaseStar.h index f7e0bc3f8..6553a29d1 100644 --- a/src/BaseStar.h +++ b/src/BaseStar.h @@ -179,7 +179,7 @@ class BaseStar { virtual void CalculateGBParams(const double p_Mass, DBL_VECTOR &p_GBParams) { } // Default is NO-OP virtual void CalculateGBParams() { CalculateGBParams(m_Mass0, m_GBParams); } // Use class member variables - virtual double CalculateImKlmTidal(const double p_Omega, const int p_l, const int p_m) { return 0.0; } // Default is 0.0 + virtual DBL_DBL_DBL_DBL CalculateImKlmTidal(const double p_Omega) { return std::make_tuple(0.0, 0.0, 0.0, 0.0); } // Default is 0.0 void CalculateLambdas() { CalculateLambdas(m_Mass - m_CoreMass); } // Use class member variables void CalculateLambdas(const double p_EnvMass); diff --git a/src/GiantBranch.cpp b/src/GiantBranch.cpp index f263339a1..21b0b9a04 100644 --- a/src/GiantBranch.cpp +++ b/src/GiantBranch.cpp @@ -1131,18 +1131,18 @@ double GiantBranch::CalculateMomentOfInertia() const { } /* - * Calculate (n=2, l, m) imaginary component of potential tidal Love number (Dynamical Tides ONLY for now) + * Calculate the (l,m) = [(1,0), (1,2), (2,2), (3,2)] imaginary components of the potential tidal Love number + * (Dynamical Tides only for now) * * Zahn, 1977, Eq. (5.5) , with the value of E_2 coming from Kushnir et al., 2017, by comparing Eq. (8) to Eq. (1) * - * double CalculateImKlmTidal(const double p_Omega) + * std::tuple CalculateImKlmTidal(const double p_Omega) * * @param [IN] p_Omega Orbital angular frequency (1/yr) - * @param [IN] l l mode - * @param [IN] m m mode - * @return Imaginary component of potential tidal love number (unitless) + * @return [(1,0), (1,2), (2,2), (3,2)] Imaginary components of the + * potential tidal love number (unitless) */ -double GiantBranch::CalculateImKlmTidal(const double p_Omega, const int p_l, const int p_m) { +std::tuple GiantBranch::CalculateImKlmTidal(const double p_Omega) { double beta2Dynamical = 1.0; double rhoFactorDynamcial = 0.1; double radiusAU = m_Radius * RSOL_TO_AU; @@ -1152,15 +1152,39 @@ double GiantBranch::CalculateImKlmTidal(const double p_Omega, const int p_l, con double coreRadius_over_radius_9 = coreRadius_over_radius_3 * coreRadius_over_radius_3 * coreRadius_over_radius_3; double mass_over_coreMass = m_Mass / m_CoreMass; - double E2Dynamical = (2.0 / 3.0) * coreRadius_over_radius_9 * mass_over_coreMass * std::cbrt(mass_over_coreMass) * beta2Dynamical * rhoFactorDynamcial; - double slm = ((p_l*p_Omega) - (p_m*Omega())) * std::sqrt(radiusAU * radiusAU * radiusAU / G_AU_Msol_yr / m_Mass); - double slm_4_3 = slm * std::cbrt(slm); - double slm_8_3 = slm_4_3 * slm_4_3; - double klmDynamical = E2Dynamical * slm_8_3; - - return klmDynamical; + double sqrt_R3_over_G_M = std::sqrt(radiusAU * radiusAU * radiusAU / G_AU_Msol_yr / m_Mass); + // (l=1, m=0), Dynamical Tide + double s10 = (p_Omega) * sqrt_R3_over_G_M; + double s10_4_3 = s10 * std::cbrt(s10); + double s10_8_3 = s10_4_3 * s10_4_3; + if (s10 < 0) s10_8_3 = - std::abs(s10_8_3); + double k10Dynamical = E2Dynamical * s10_8_3; + + // (l=1, m=2), Dynamical Tide + double s12 = (p_Omega - (2*Omega())) * sqrt_R3_over_G_M; + double s12_4_3 = s12 * std::cbrt(s12); + double s12_8_3 = s12_4_3 * s12_4_3; + if (s12 < 0) s12_8_3 = - std::abs(s12_8_3); + double k12Dynamical = E2Dynamical * s12_8_3; + + // (l=2, m=2), Dynamical Tide + double s22 = ((2*p_Omega) - (2*Omega())) * sqrt_R3_over_G_M; + double s22_4_3 = s22 * std::cbrt(s22); + double s22_8_3 = s22_4_3 * s22_4_3; + if (s22 < 0) s22_8_3 = - std::abs(s22_8_3); + double k22Dynamical = E2Dynamical * s22_8_3; + + // (l=3, m=2), Dynamical Tide + double s32 = ((3*p_Omega) - (2*Omega())) * sqrt_R3_over_G_M; + double s32_4_3 = s32 * std::cbrt(s32); + double s32_8_3 = s32_4_3 * s32_4_3; + if (s32 < 0) s32_8_3 = - std::abs(s32_8_3); + double k32Dynamical = E2Dynamical * s32_8_3; + + // return klmDynamical; + return std::make_tuple(k10Dynamical, k12Dynamical, k22Dynamical, k32Dynamical); } diff --git a/src/GiantBranch.h b/src/GiantBranch.h index ce50e5215..f244994a6 100755 --- a/src/GiantBranch.h +++ b/src/GiantBranch.h @@ -52,7 +52,7 @@ class GiantBranch: virtual public BaseStar, public MainSequence { void CalculateGBParams() { CalculateGBParams(m_Mass0, m_GBParams); } // Use class member variables static double CalculateHRateConstant_Static(const double p_Mass); - double CalculateImKlmTidal(const double p_Omega, const int p_l, const int p_m); + DBL_DBL_DBL_DBL CalculateImKlmTidal(const double p_Omega); virtual double CalculateInitialSupernovaMass() const { return m_Mass; } // Use class member variables double CalculateLifetimeToHeIgnition(const double p_Mass, const double p_Tinf1_FGB, const double p_Tinf2_FGB); diff --git a/src/Star.h b/src/Star.h index 333adf533..224e69e45 100755 --- a/src/Star.h +++ b/src/Star.h @@ -162,7 +162,7 @@ class Star { double CalculateConvectiveEnvelopeMass() { return m_Star->CalculateConvectiveEnvelopeMass(); } double CalculateEddyTurnoverTimescale() { return m_Star->CalculateEddyTurnoverTimescale(); } - double CalculateImKlmTidal(const double p_Omega, const int p_l, const int p_m) { return m_Star->CalculateImKlmTidal(p_Omega, p_l, p_m); } + DBL_DBL_DBL_DBL CalculateImKlmTidal(const double p_Omega) { return m_Star->CalculateImKlmTidal(p_Omega); } void CalculateLambdas() { m_Star->CalculateLambdas(); } void CalculateLambdas(const double p_EnvMass) { m_Star->CalculateLambdas(p_EnvMass); } diff --git a/src/constants.h b/src/constants.h index 77c1ca1b0..3d8afd32a 100755 --- a/src/constants.h +++ b/src/constants.h @@ -38,6 +38,7 @@ typedef unsigned long int OBJECT_I typedef std::vector DBL_VECTOR; typedef std::tuple DBL_DBL; typedef std::tuple DBL_DBL_DBL; +typedef std::tuple DBL_DBL_DBL_DBL; typedef std::tuple STR_STR_STR_STR; // Hash for Enum Class