Skip to content

Commit ad8c2ed

Browse files
veome22veome22
committed
Promoted tidal ImKnm to global variables to compute timescales as outputs, added all tidal quantities to outputs
1 parent b10b3e8 commit ad8c2ed

File tree

3 files changed

+170
-40
lines changed

3 files changed

+170
-40
lines changed

‎src/BaseBinaryStar.cpp‎

Lines changed: 65 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -398,8 +398,11 @@ void BaseBinaryStar::SetRemainingValues() {
398398
m_PhiE = DEFAULT_INITIAL_DOUBLE_VALUE;
399399
m_PsiE = DEFAULT_INITIAL_DOUBLE_VALUE;
400400

401-
m_SynchronizationTimescale = DEFAULT_INITIAL_DOUBLE_VALUE;
401+
m_SynchronizationTimescale1 = DEFAULT_INITIAL_DOUBLE_VALUE;
402+
m_SynchronizationTimescale2 = DEFAULT_INITIAL_DOUBLE_VALUE;
402403
m_CircularizationTimescale = DEFAULT_INITIAL_DOUBLE_VALUE;
404+
m_ImKnm1_tidal = std::tuple<double, double, double, double> (DEFAULT_INITIAL_DOUBLE_VALUE, DEFAULT_INITIAL_DOUBLE_VALUE, DEFAULT_INITIAL_DOUBLE_VALUE, DEFAULT_INITIAL_DOUBLE_VALUE);;
405+
m_ImKnm2_tidal = std::tuple<double, double, double, double> (DEFAULT_INITIAL_DOUBLE_VALUE, DEFAULT_INITIAL_DOUBLE_VALUE, DEFAULT_INITIAL_DOUBLE_VALUE, DEFAULT_INITIAL_DOUBLE_VALUE);;
403406

404407
// RLOF details
405408
m_RLOFDetails.experiencedRLOF = false;
@@ -625,11 +628,36 @@ COMPAS_VARIABLE BaseBinaryStar::BinaryPropertyValue(const T_ANY_PROPERTY p_Prope
625628
case BINARY_PROPERTY::STELLAR_TYPE_NAME_2_PRE_COMMON_ENVELOPE: value = STELLAR_TYPE_LABEL.at(StellarType2PreCEE()); break;
626629
case BINARY_PROPERTY::SUPERNOVA_ORBIT_INCLINATION_ANGLE: value = SN_OrbitInclinationAngle(); break;
627630
case BINARY_PROPERTY::SUPERNOVA_STATE: value = SN_State(); break;
628-
case BINARY_PROPERTY::SYNCHRONIZATION_TIMESCALE: value = SynchronizationTimescale(); break;
631+
case BINARY_PROPERTY::SYNCHRONIZATION_TIMESCALE_1: value = SynchronizationTimescale1(); break;
632+
case BINARY_PROPERTY::SYNCHRONIZATION_TIMESCALE_2: value = SynchronizationTimescale2(); break;
629633
case BINARY_PROPERTY::SYSTEMIC_SPEED: value = SystemicSpeed(); break;
630634
case BINARY_PROPERTY::SYSTEMIC_VELOCITY_X: value = SystemicVelocityX(); break;
631635
case BINARY_PROPERTY::SYSTEMIC_VELOCITY_Y: value = SystemicVelocityY(); break;
632636
case BINARY_PROPERTY::SYSTEMIC_VELOCITY_Z: value = SystemicVelocityZ(); break;
637+
case BINARY_PROPERTY::TIDAL_POTENTIAL_LOVE_NUMBER_1_10: std::tie(value, std::ignore, std::ignore, std::ignore) = ImKnm1_tidal(); break;
638+
case BINARY_PROPERTY::TIDAL_POTENTIAL_LOVE_NUMBER_1_12: std::tie(std::ignore, value, std::ignore, std::ignore) = ImKnm1_tidal(); break;
639+
case BINARY_PROPERTY::TIDAL_POTENTIAL_LOVE_NUMBER_1_22: std::tie(std::ignore, std::ignore, value, std::ignore) = ImKnm1_tidal(); break;
640+
case BINARY_PROPERTY::TIDAL_POTENTIAL_LOVE_NUMBER_1_32: std::tie(std::ignore, std::ignore, std::ignore, value) = ImKnm1_tidal(); break;
641+
case BINARY_PROPERTY::TIDAL_POTENTIAL_LOVE_NUMBER_2_10: std::tie(value, std::ignore, std::ignore, std::ignore) = ImKnm2_tidal(); break;
642+
case BINARY_PROPERTY::TIDAL_POTENTIAL_LOVE_NUMBER_2_12: std::tie(std::ignore, value, std::ignore, std::ignore) = ImKnm2_tidal(); break;
643+
case BINARY_PROPERTY::TIDAL_POTENTIAL_LOVE_NUMBER_2_22: std::tie(std::ignore, std::ignore, value, std::ignore) = ImKnm2_tidal(); break;
644+
case BINARY_PROPERTY::TIDAL_POTENTIAL_LOVE_NUMBER_2_32: std::tie(std::ignore, std::ignore, std::ignore, value) = ImKnm2_tidal(); break;
645+
case BINARY_PROPERTY::TIDAL_POTENTIAL_LOVE_NUMBER_1_10_EQ: std::tie(value, std::ignore, std::ignore, std::ignore) = ImKnm1_tidal_eq(); break;
646+
case BINARY_PROPERTY::TIDAL_POTENTIAL_LOVE_NUMBER_1_12_EQ: std::tie(std::ignore, value, std::ignore, std::ignore) = ImKnm1_tidal_eq(); break;
647+
case BINARY_PROPERTY::TIDAL_POTENTIAL_LOVE_NUMBER_1_22_EQ: std::tie(std::ignore, std::ignore, value, std::ignore) = ImKnm1_tidal_eq(); break;
648+
case BINARY_PROPERTY::TIDAL_POTENTIAL_LOVE_NUMBER_1_32_EQ: std::tie(std::ignore, std::ignore, std::ignore, value) = ImKnm1_tidal_eq(); break;
649+
case BINARY_PROPERTY::TIDAL_POTENTIAL_LOVE_NUMBER_2_10_EQ: std::tie(value, std::ignore, std::ignore, std::ignore) = ImKnm2_tidal_eq(); break;
650+
case BINARY_PROPERTY::TIDAL_POTENTIAL_LOVE_NUMBER_2_12_EQ: std::tie(std::ignore, value, std::ignore, std::ignore) = ImKnm2_tidal_eq(); break;
651+
case BINARY_PROPERTY::TIDAL_POTENTIAL_LOVE_NUMBER_2_22_EQ: std::tie(std::ignore, std::ignore, value, std::ignore) = ImKnm2_tidal_eq(); break;
652+
case BINARY_PROPERTY::TIDAL_POTENTIAL_LOVE_NUMBER_2_32_EQ: std::tie(std::ignore, std::ignore, std::ignore, value) = ImKnm2_tidal_eq(); break;
653+
case BINARY_PROPERTY::TIDAL_POTENTIAL_LOVE_NUMBER_1_10_DYN: std::tie(value, std::ignore, std::ignore, std::ignore) = ImKnm1_tidal_dyn();break;
654+
case BINARY_PROPERTY::TIDAL_POTENTIAL_LOVE_NUMBER_1_12_DYN: std::tie(std::ignore, value, std::ignore, std::ignore) = ImKnm1_tidal_dyn();break;
655+
case BINARY_PROPERTY::TIDAL_POTENTIAL_LOVE_NUMBER_1_22_DYN: std::tie(std::ignore, std::ignore, value, std::ignore) = ImKnm1_tidal_dyn();break;
656+
case BINARY_PROPERTY::TIDAL_POTENTIAL_LOVE_NUMBER_1_32_DYN: std::tie(std::ignore, std::ignore, std::ignore, value) = ImKnm1_tidal_dyn();break;
657+
case BINARY_PROPERTY::TIDAL_POTENTIAL_LOVE_NUMBER_2_10_DYN: std::tie(value, std::ignore, std::ignore, std::ignore) = ImKnm2_tidal_dyn();break;
658+
case BINARY_PROPERTY::TIDAL_POTENTIAL_LOVE_NUMBER_2_12_DYN: std::tie(std::ignore, value, std::ignore, std::ignore) = ImKnm2_tidal_dyn();break;
659+
case BINARY_PROPERTY::TIDAL_POTENTIAL_LOVE_NUMBER_2_22_DYN: std::tie(std::ignore, std::ignore, value, std::ignore) = ImKnm2_tidal_dyn();break;
660+
case BINARY_PROPERTY::TIDAL_POTENTIAL_LOVE_NUMBER_2_32_DYN: std::tie(std::ignore, std::ignore, std::ignore, value) = ImKnm2_tidal_dyn();break;
633661
case BINARY_PROPERTY::TIME: value = Time(); break;
634662
case BINARY_PROPERTY::TIME_TO_COALESCENCE: value = TimeToCoalescence(); break;
635663
case BINARY_PROPERTY::TOTAL_ANGULAR_MOMENTUM: value = TotalAngularMomentum(); break;
@@ -2705,31 +2733,35 @@ void BaseBinaryStar::ProcessTides(const double p_Dt) {
27052733

27062734
// Evolve binary semi-major axis, eccentricity, and spin of each star based on Kapil et al., 2025
27072735

2708-
DBL_DBL_DBL_DBL ImKnm1 = m_Star1->CalculateImKnmTidal(omega, m_SemiMajorAxis, m_Star2->Mass());
2709-
DBL_DBL_DBL_DBL ImKnm2 = m_Star2->CalculateImKnmTidal(omega, m_SemiMajorAxis, m_Star1->Mass());
2736+
m_ImKnm1_tidal = m_Star1->CalculateImKnmTidal(omega, m_SemiMajorAxis, m_Star2->Mass());
2737+
m_ImKnm2_tidal = m_Star2->CalculateImKnmTidal(omega, m_SemiMajorAxis, m_Star1->Mass());
27102738

2711-
double DSemiMajorAxis1Dt = CalculateDSemiMajorAxisTidalDt(ImKnm1, m_Star1); // change in semi-major axis from star1
2712-
double DSemiMajorAxis2Dt = CalculateDSemiMajorAxisTidalDt(ImKnm2, m_Star2); // change in semi-major axis from star2
2739+
m_DSemiMajorAxis1Dt_tidal = CalculateDSemiMajorAxisTidalDt(m_ImKnm1_tidal, m_Star1); // change in semi-major axis from star1
2740+
m_DSemiMajorAxis2Dt_tidal = CalculateDSemiMajorAxisTidalDt(m_ImKnm2_tidal, m_Star2); // change in semi-major axis from star2
27132741

2714-
double DEccentricity1Dt = CalculateDEccentricityTidalDt(ImKnm1, m_Star1); // change in eccentricity from star1
2715-
double DEccentricity2Dt = CalculateDEccentricityTidalDt(ImKnm2, m_Star2); // change in eccentricity from star2
2742+
m_DEccentricity1Dt_tidal = CalculateDEccentricityTidalDt(m_ImKnm1_tidal, m_Star1); // change in eccentricity from star1
2743+
m_DEccentricity2Dt_tidal = CalculateDEccentricityTidalDt(m_ImKnm2_tidal, m_Star2); // change in eccentricity from star2
27162744

2717-
double DOmega1Dt = CalculateDOmegaTidalDt(ImKnm1, m_Star1); // change in spin from star1
2718-
double DOmega2Dt = CalculateDOmegaTidalDt(ImKnm2, m_Star2); // change in spin from star2
2745+
m_DOmega1Dt_tidal = CalculateDOmegaTidalDt(m_ImKnm1_tidal, m_Star1); // change in spin from star1
2746+
m_DOmega2Dt_tidal = CalculateDOmegaTidalDt(m_ImKnm2_tidal, m_Star2); // change in spin from star2
27192747

27202748
// limit change in stellar and orbital properties from tides to a maximum fraction of the current value
27212749
double fraction_tidal_change = 1.0;
2722-
fraction_tidal_change = std::min(fraction_tidal_change, std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * OrbitalAngularVelocity() / (DOmega1Dt * p_Dt * MYR_TO_YEAR)));
2723-
fraction_tidal_change = std::min(fraction_tidal_change, std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * OrbitalAngularVelocity() / (DOmega2Dt * p_Dt * MYR_TO_YEAR)));
2724-
fraction_tidal_change = std::min(fraction_tidal_change, std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * m_SemiMajorAxis / ((DSemiMajorAxis1Dt + DSemiMajorAxis2Dt) * p_Dt * MYR_TO_YEAR)));
2725-
fraction_tidal_change = std::min(fraction_tidal_change, std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * m_Eccentricity / ((DEccentricity1Dt + DEccentricity2Dt) * p_Dt * MYR_TO_YEAR)));
2750+
fraction_tidal_change = std::min(fraction_tidal_change, std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * OrbitalAngularVelocity() / (m_DOmega1Dt_tidal * p_Dt * MYR_TO_YEAR)));
2751+
fraction_tidal_change = std::min(fraction_tidal_change, std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * OrbitalAngularVelocity() / (m_DOmega2Dt_tidal * p_Dt * MYR_TO_YEAR)));
2752+
fraction_tidal_change = std::min(fraction_tidal_change, std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * m_SemiMajorAxis / ((m_DSemiMajorAxis1Dt_tidal + m_DSemiMajorAxis2Dt_tidal) * p_Dt * MYR_TO_YEAR)));
2753+
fraction_tidal_change = std::min(fraction_tidal_change, std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * m_Eccentricity / ((m_DEccentricity1Dt_tidal + m_DEccentricity2Dt_tidal) * p_Dt * MYR_TO_YEAR)));
27262754

2727-
m_Star1->SetOmega(m_Star1->Omega() + fraction_tidal_change * (DOmega1Dt * p_Dt * MYR_TO_YEAR)); // evolve star 1 spin
2728-
m_Star2->SetOmega(m_Star2->Omega() + fraction_tidal_change * (DOmega2Dt * p_Dt * MYR_TO_YEAR)); // evolve star 2 spin
2729-
m_SemiMajorAxis = m_SemiMajorAxis + fraction_tidal_change * ((DSemiMajorAxis1Dt + DSemiMajorAxis2Dt) * p_Dt * MYR_TO_YEAR); // evolve separation
2730-
m_Eccentricity = m_Eccentricity + fraction_tidal_change * ((DEccentricity1Dt + DEccentricity2Dt) * p_Dt * MYR_TO_YEAR); // evolve eccentricity
2755+
m_Star1->SetOmega(m_Star1->Omega() + fraction_tidal_change * (m_DOmega1Dt_tidal * p_Dt * MYR_TO_YEAR)); // evolve star 1 spin
2756+
m_Star2->SetOmega(m_Star2->Omega() + fraction_tidal_change * (m_DOmega2Dt_tidal * p_Dt * MYR_TO_YEAR)); // evolve star 2 spin
2757+
m_SemiMajorAxis = m_SemiMajorAxis + fraction_tidal_change * ((m_DSemiMajorAxis1Dt_tidal + m_DSemiMajorAxis2Dt_tidal) * p_Dt * MYR_TO_YEAR); // evolve separation
2758+
m_Eccentricity = m_Eccentricity + fraction_tidal_change * ((m_DEccentricity1Dt_tidal + m_DEccentricity2Dt_tidal) * p_Dt * MYR_TO_YEAR); // evolve eccentricity
27312759

2732-
m_TotalAngularMomentum = CalculateAngularMomentum(); // re-calculate angular momenta
2760+
m_CircularizationTimescale = - m_Eccentricity / (m_DEccentricity1Dt_tidal + m_DEccentricity2Dt_tidal) * YEAR_TO_MYR; // Circularization timescale in Myr (for output files)
2761+
m_SynchronizationTimescale1 = - (m_Star1->Omega() - omega) / m_DOmega1Dt_tidal * YEAR_TO_MYR; // Synchronization timescale for Star1 in Myr (for output files)
2762+
m_SynchronizationTimescale2 = - (m_Star2->Omega() - omega) / m_DOmega2Dt_tidal * YEAR_TO_MYR; // Synchronization timescale for Star2 in Myr (for output files)
2763+
2764+
m_TotalAngularMomentum = CalculateAngularMomentum(); // re-calculate angular momenta
27332765
m_OrbitalAngularMomentum = CalculateOrbitalAngularMomentum(m_Star1->Mass(), m_Star2->Mass(), m_SemiMajorAxis, m_Eccentricity);
27342766

27352767
} break;
@@ -2999,30 +3031,30 @@ double BaseBinaryStar::ChooseTimestep(const double p_Factor) {
29993031
// yes - need to adjust dt
30003032
double omega = OrbitalAngularVelocity();
30013033

3002-
DBL_DBL_DBL_DBL ImKnm1 = m_Star1->CalculateImKnmTidal(omega, m_SemiMajorAxis, m_Star2->Mass());
3003-
DBL_DBL_DBL_DBL ImKnm2 = m_Star2->CalculateImKnmTidal(omega, m_SemiMajorAxis, m_Star1->Mass());
3034+
m_ImKnm1_tidal = m_Star1->CalculateImKnmTidal(omega, m_SemiMajorAxis, m_Star2->Mass());
3035+
m_ImKnm2_tidal = m_Star2->CalculateImKnmTidal(omega, m_SemiMajorAxis, m_Star1->Mass());
30043036

3005-
double DSemiMajorAxis1DtTidal = CalculateDSemiMajorAxisTidalDt(ImKnm1, m_Star1);
3006-
double DSemiMajorAxis2DtTidal = CalculateDSemiMajorAxisTidalDt(ImKnm2, m_Star2);
3037+
double m_DSemiMajorAxis1Dt_tidal = CalculateDSemiMajorAxisTidalDt(m_ImKnm1_tidal, m_Star1);
3038+
double m_DSemiMajorAxis2Dt_tidal = CalculateDSemiMajorAxisTidalDt(m_ImKnm2_tidal, m_Star2);
30073039

3008-
double DEccentricity1DtTidal = CalculateDEccentricityTidalDt(ImKnm1, m_Star1);
3009-
double DEccentricity2DtTidal = CalculateDEccentricityTidalDt(ImKnm2, m_Star2);
3040+
double m_DEccentricity1Dt_tidal = CalculateDEccentricityTidalDt(m_ImKnm1_tidal, m_Star1);
3041+
double m_DEccentricity2Dt_tidal = CalculateDEccentricityTidalDt(m_ImKnm2_tidal, m_Star2);
30103042

3011-
double DOmega1Dt_tidal = CalculateDOmegaTidalDt(ImKnm1, m_Star1);
3012-
double DOmega2Dt_tidal = CalculateDOmegaTidalDt(ImKnm2, m_Star2);
3043+
double m_DOmega1Dt_tidal = CalculateDOmegaTidalDt(m_ImKnm1_tidal, m_Star1);
3044+
double m_DOmega2Dt_tidal = CalculateDOmegaTidalDt(m_ImKnm2_tidal, m_Star2);
30133045

30143046
// Ensure that the change in orbital and spin properties due to tides in a single timestep is constrained (to 1 percent by default)
30153047
// Limit the spin evolution of each star based on the orbital frequency rather than its spin frequency, since tides should not cause major problems until synchronization.
3016-
double Dt_SemiMajorAxis1Tidal = utils::Compare(DSemiMajorAxis1DtTidal, 0.0) == 0 ? dt : std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * m_SemiMajorAxis / DSemiMajorAxis1DtTidal) * YEAR_TO_MYR;
3017-
double Dt_SemiMajorAxis2Tidal = utils::Compare(DSemiMajorAxis2DtTidal, 0.0) == 0 ? dt : std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * m_SemiMajorAxis / DSemiMajorAxis2DtTidal) * YEAR_TO_MYR;
3048+
double Dt_SemiMajorAxis1Tidal = utils::Compare(m_DSemiMajorAxis1Dt_tidal, 0.0) == 0 ? dt : std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * m_SemiMajorAxis / m_DSemiMajorAxis1Dt_tidal) * YEAR_TO_MYR;
3049+
double Dt_SemiMajorAxis2Tidal = utils::Compare(m_DSemiMajorAxis2Dt_tidal, 0.0) == 0 ? dt : std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * m_SemiMajorAxis / m_DSemiMajorAxis2Dt_tidal) * YEAR_TO_MYR;
30183050
double Dt_SemiMajorAxisTidal = std::min(Dt_SemiMajorAxis1Tidal, Dt_SemiMajorAxis2Tidal);
30193051

3020-
double Dt_Eccentricity1Tidal = utils::Compare(DEccentricity1DtTidal, 0.0) == 0 ? dt : std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * m_Eccentricity / DEccentricity1DtTidal) * YEAR_TO_MYR;
3021-
double Dt_Eccentricity2Tidal = utils::Compare(DEccentricity2DtTidal, 0.0) == 0 ? dt : std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * m_Eccentricity / DEccentricity2DtTidal) * YEAR_TO_MYR;
3052+
double Dt_Eccentricity1Tidal = utils::Compare(m_DEccentricity1Dt_tidal, 0.0) == 0 ? dt : std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * m_Eccentricity / m_DEccentricity1Dt_tidal) * YEAR_TO_MYR;
3053+
double Dt_Eccentricity2Tidal = utils::Compare(m_DEccentricity2Dt_tidal, 0.0) == 0 ? dt : std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * m_Eccentricity / m_DEccentricity2Dt_tidal) * YEAR_TO_MYR;
30223054
double Dt_EccentricityTidal = std::min(Dt_Eccentricity1Tidal, Dt_Eccentricity2Tidal);
30233055

3024-
double Dt_Omega1Tidal = utils::Compare(DOmega1Dt_tidal, 0.0) == 0 ? dt : std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * omega / DOmega1Dt_tidal) * YEAR_TO_MYR;
3025-
double Dt_Omega2Tidal = utils::Compare(DOmega2Dt_tidal, 0.0) == 0 ? dt : std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * omega / DOmega2Dt_tidal) * YEAR_TO_MYR;
3056+
double Dt_Omega1Tidal = utils::Compare(m_DOmega1Dt_tidal, 0.0) == 0 ? dt : std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * omega / m_DOmega1Dt_tidal) * YEAR_TO_MYR;
3057+
double Dt_Omega2Tidal = utils::Compare(m_DOmega2Dt_tidal, 0.0) == 0 ? dt : std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * omega / m_DOmega2Dt_tidal) * YEAR_TO_MYR;
30263058
double Dt_OmegaTidal = std::min(Dt_Omega1Tidal, Dt_Omega2Tidal);
30273059

30283060
dt = std::min(dt, std::min(Dt_SemiMajorAxisTidal, std::min(Dt_EccentricityTidal, Dt_OmegaTidal)));

0 commit comments

Comments
 (0)