Skip to content

Commit

Permalink
Modified CalculateImKlmTidal() to return a tuple with the appropriate…
Browse files Browse the repository at this point in the history 
… moments of ImK in order
  • Loading branch information
@veome22
veome22 committed Apr 10, 2024
1 parent cc1be8c commit b770c99
Show file tree
Hide file tree
Showing 7 changed files with 72 additions and 50 deletions.
22 changes: 9 additions & 13 deletions src/BaseBinaryStar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2342,45 +2342,42 @@ 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(),
m_Omega,
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(),
m_Omega,
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(),
m_Omega,
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(),
m_Omega,
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(),
Expand All @@ -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(),
Expand Down
43 changes: 22 additions & 21 deletions src/BaseBinaryStar.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 <double, double, double, double> 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)
Expand All @@ -788,25 +785,26 @@ 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 <double, double, double, double> 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))));
}

/*
* Change in eccentricity based on secular equations for tidal evolution given the tidal love number
* 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 <double, double, double, double> 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)
Expand All @@ -815,26 +813,26 @@ 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 <double, double, double, double> 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));
}

/*
* Change in spin based on secular equations for tidal evolution given the tidal love number
* 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 <double, double, double, double> 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)
Expand All @@ -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 <double, double, double, double> 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))));
}

};
Expand Down
2 changes: 1 addition & 1 deletion src/BaseStar.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
50 changes: 37 additions & 13 deletions src/GiantBranch.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 <double, double, double, double> 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 <double, double, double, double> GiantBranch::CalculateImKlmTidal(const double p_Omega) {
double beta2Dynamical = 1.0;
double rhoFactorDynamcial = 0.1;
double radiusAU = m_Radius * RSOL_TO_AU;
Expand All @@ -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);
}


Expand Down
2 changes: 1 addition & 1 deletion src/GiantBranch.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
Loading

0 comments on commit b770c99

Please sign in to comment.