Skip to content

Commit 56ba593

Browse files
veome22veome22
committed
Updated dynamical tides to match paper, set m=2 everywhere
1 parent 76b6fce commit 56ba593

File tree

1 file changed

+19
-13
lines changed

1 file changed

+19
-13
lines changed

‎src/BaseStar.cpp‎

Lines changed: 19 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -3234,7 +3234,7 @@ DBL_DBL_DBL_DBL BaseStar::CalculateImKnmDynamical(const double p_Omega, const do
32343234
double radiusAU = m_Radius * RSOL_TO_AU;
32353235
double coreRadiusAU = CalculateConvectiveCoreRadius() * RSOL_TO_AU;
32363236
double convectiveEnvRadiusAU = CalculateRadialExtentConvectiveEnvelope() * RSOL_TO_AU;
3237-
double radiusIntershellAU = radiusAU - convectiveEnvRadiusAU; // Outer radial coordinate of radiative intershell
3237+
double radiusIntershellAU = radiusAU - convectiveEnvRadiusAU; // Outer radial coordinate of radiative intershell, or inner radial coordinate of convective envelope
32383238

32393239
// There should be no Dynamical tides if the entire star is convective, i.e. if there are no convective-radiative boundaries.
32403240
// If so, return 0.0 for all dynamical components of ImKnm.
@@ -3275,34 +3275,34 @@ DBL_DBL_DBL_DBL BaseStar::CalculateImKnmDynamical(const double p_Omega, const do
32753275
double coreRadiusOverRadius_3 = coreRadiusOverRadius * coreRadiusOverRadius * coreRadiusOverRadius;
32763276
double coreRadiusOverRadius_9 = coreRadiusOverRadius_3 * coreRadiusOverRadius_3 * coreRadiusOverRadius_3;
32773277
double massOverCoreMass = m_Mass / coreMass;
3278-
double E2Dynamical = (2.0 / 3.0) * coreRadiusOverRadius_9 * massOverCoreMass * std::cbrt(massOverCoreMass) * beta2Dynamical * rhoFactorDynamcial;
3278+
double E2Core = (2.0 / 3.0) * coreRadiusOverRadius_9 * massOverCoreMass * std::cbrt(massOverCoreMass) * beta2Dynamical * rhoFactorDynamcial;
32793279

32803280
// (l=2, n=1, m=0), Gravity Wave dissipation from core boundary
32813281
double s10 = w10 * sqrtR3OverG_M;
32823282
double s10_4_3 = s10 * std::cbrt(s10);
32833283
double s10_8_3 = s10_4_3 * s10_4_3;
3284-
k10GravityCore = E2Dynamical * std::copysign(s10_8_3, w10);
3284+
k10GravityCore = E2Core * std::copysign(s10_8_3, w10);
32853285
if (std::isnan(k10GravityCore)) k10GravityCore = 0.0;
32863286

32873287
// (l=2, n=1, m=2), Gravity Wave dissipation from core boundary
32883288
double s12 = w12 * sqrtR3OverG_M;
32893289
double s12_4_3 = s12 * std::cbrt(s12);
32903290
double s12_8_3 = s12_4_3 * s12_4_3;
3291-
k12GravityCore = E2Dynamical * std::copysign(s12_8_3, w12);
3291+
k12GravityCore = E2Core * std::copysign(s12_8_3, w12);
32923292
if (std::isnan(k12GravityCore)) k12GravityCore = 0.0;
32933293

32943294
// (l=2, n=2, m=2), Gravity Wave dissipation from core boundary
32953295
double s22 = w22 * sqrtR3OverG_M;
32963296
double s22_4_3 = s22 * std::cbrt(s22);
32973297
double s22_8_3 = s22_4_3 * s22_4_3;
3298-
k22GravityCore = E2Dynamical * std::copysign(s22_8_3, w22);
3298+
k22GravityCore = E2Core * std::copysign(s22_8_3, w22);
32993299
if (std::isnan(k22GravityCore)) k22GravityCore = 0.0;
33003300

33013301
// (l=2, n=3, m=2), Gravity Wave dissipation from core boundary
33023302
double s32 = w32 * sqrtR3OverG_M;
33033303
double s32_4_3 = s32 * std::cbrt(s32);
33043304
double s32_8_3 = s32_4_3 * s32_4_3;
3305-
k32GravityCore = E2Dynamical * std::copysign(s32_8_3, w32);
3305+
k32GravityCore = E2Core * std::copysign(s32_8_3, w32);
33063306
if (std::isnan(k32GravityCore)) k32GravityCore = 0.0;
33073307
}
33083308

@@ -3314,8 +3314,8 @@ DBL_DBL_DBL_DBL BaseStar::CalculateImKnmDynamical(const double p_Omega, const do
33143314
if ((utils::Compare(convectiveEnvRadiusAU / radiusAU, TIDES_MINIMUM_FRACTIONAL_EXTENT) > 0) || (utils::Compare(envMass / m_Mass, TIDES_MINIMUM_FRACTIONAL_EXTENT) > 0)) {
33153315

33163316
constexpr double dynPrefactor = 3.207452512782476; // 3^(11/3) * Gamma(1/3)^2 / 40 PI
3317-
constexpr double m_l_factor_22 = 0.183440402716368; // m * (l(l+1))^{-4/3}
3318-
double cbrtdNdlnr = std::cbrt(G_AU_Msol_yr * radIntershellMass / radiusIntershellAU / (radiusAU - radiusIntershellAU) / (radiusAU - radiusIntershellAU));
3317+
constexpr double m_l_factor_22 = 0.183440402716368; // m * (l(l+1))^{-4/3}, assuming l=2, m=2
3318+
double cbrtdNdlnr = std::cbrt(G_AU_Msol_yr * radIntershellMass / radiusIntershellAU / radiusIntershellAU / (radiusAU - radiusIntershellAU));
33193319

33203320
double alpha = radiusIntershellAU / radiusAU;
33213321
double oneMinusAlpha = 1.0 - alpha;
@@ -3335,27 +3335,33 @@ DBL_DBL_DBL_DBL BaseStar::CalculateImKnmDynamical(const double p_Omega, const do
33353335
double alpha_2_3Minus_1 = (alpha * 2.0 / 3.0) - 1.0;
33363336

33373337
// Assume GW dissipation from the envelope boundary only acts if the radiative zone extends to the core, i.e. if there is no convective core.
3338-
if (utils::Compare(coreRadiusAU/radiusAU, TIDES_MINIMUM_FRACTIONAL_EXTENT) < 0) {
3338+
if (utils::Compare(coreRadiusAU/radiusAU, TIDES_MINIMUM_FRACTIONAL_EXTENT) < 0) {
33393339
double Epsilon = alpha_11 * envMass / m_Mass * oneMinusGamma_2 * alpha_2_3Minus_1 * alpha_2_3Minus_1 / beta_2 / oneMinusAlpha_3 / oneMinusAlpha_2;
3340+
double R3_Epsilon_dNdlnr_factor = R3OverG_M * Epsilon / cbrtdNdlnr;
3341+
double E2Envelope = dynPrefactor * m_l_factor_22 * R3_Epsilon_dNdlnr_factor;
33403342

3341-
// (l=2, n=1, m=0), Gravity Wave dissipation from envelope boundary is always 0.0 since m * (l(l+1))^{-4/3} = 0
3343+
// (l=2, n=1, m=0), Gravity Wave dissipation from envelope boundary
3344+
double w10_4_3 = w10 * std::cbrt(w10);
3345+
double w10_8_3 = w10_4_3 * w10_4_3;
3346+
k10GravityEnv = E2Envelope * std::copysign(w10_8_3, w10);
3347+
if (std::isnan(k10GravityEnv)) k10GravityEnv = 0.0;
33423348

33433349
// (l=2, n=1, m=2), Gravity Wave dissipation from envelope boundary
33443350
double w12_4_3 = w12 * std::cbrt(w12);
33453351
double w12_8_3 = w12_4_3 * w12_4_3;
3346-
k12GravityEnv = dynPrefactor * m_l_factor_22 * std::copysign(w12_8_3, w12) * R3OverG_M * Epsilon / cbrtdNdlnr;
3352+
k12GravityEnv = E2Envelope * std::copysign(w12_8_3, w12);
33473353
if (std::isnan(k12GravityEnv)) k12GravityEnv = 0.0;
33483354

33493355
// (l=2, n=2, m=2), Gravity Wave dissipation from envelope boundary
33503356
double w22_4_3 = w22 * std::cbrt(w22);
33513357
double w22_8_3 = w22_4_3 * w22_4_3;
3352-
k22GravityEnv = dynPrefactor * m_l_factor_22 * std::copysign(w22_8_3, w22)* R3OverG_M * Epsilon / cbrtdNdlnr;
3358+
k22GravityEnv = E2Envelope * std::copysign(w22_8_3, w22);
33533359
if (std::isnan(k22GravityEnv)) k22GravityEnv = 0.0;
33543360

33553361
// (l=2, n=3, m=2), Gravity Wave dissipation from envelope boundary
33563362
double w32_4_3 = w32 * std::cbrt(w32);
33573363
double w32_8_3 = w32_4_3 * w32_4_3;
3358-
k32GravityEnv = dynPrefactor * m_l_factor_22 * std::copysign(w32_8_3, w32) * R3OverG_M * Epsilon / cbrtdNdlnr;
3364+
k32GravityEnv = E2Envelope * std::copysign(w32_8_3, w32);
33593365
if (std::isnan(k32GravityEnv)) k32GravityEnv = 0.0;
33603366
}
33613367

0 commit comments

Comments
 (0)