Skip to content

Commit ce60b5a

Browse files
veome22veome22
committed
added Zahn 1977 and Hurley 2002 prescription for testing
1 parent 7af07e6 commit ce60b5a

File tree

2 files changed

+138
-5
lines changed

2 files changed

+138
-5
lines changed

‎src/BaseBinaryStar.cpp‎

Lines changed: 135 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -2857,6 +2857,139 @@ void BaseBinaryStar::ProcessTides(const double p_Dt) {
28572857
}
28582858
} break;
28592859

2860+
case TIDES_PRESCRIPTION::ZAHN1977: {
2861+
// Evolve binary semi-major axis, eccentricity, and spin of each star based on Zahn (1977)
2862+
// Porb = 2*np.pi / omega_orb
2863+
// Pspin = 2*np.pi / omega_spin
2864+
// Ptid_inv = np.abs(1/Porb - 1/Pspin) # Eq. (33)
2865+
// Ptid = 1/Ptid_inv # Eq. (33)
2866+
2867+
// fconv = np.minimum(1, (Ptid/(2*tau_conv))**2) # Eq. (32)
2868+
// k_over_T_c = (2/21) * (fconv/tau_conv) * Menv/M # Eq. (30)
2869+
// double ImKnm1_Zahn = Imk22_eq_zahn = 2 * (k_over_T_c) * (R_AU**3 / (G_AU_Msol_yr * M)) * (2*omega_orb - 2*omega_spin)
2870+
double envMass1, envMassMax1;
2871+
std::tie(envMass1, envMassMax1) = m_Star1->CalculateConvectiveEnvelopeMass();
2872+
2873+
double envMass2, envMassMax2;
2874+
std::tie(envMass2, envMassMax2) = m_Star2->CalculateConvectiveEnvelopeMass();
2875+
2876+
double pOrbital = _2_PI / omega;
2877+
double pSpin1 = _2_PI / m_Star1->Omega();
2878+
double pSpin2 = _2_PI / m_Star2->Omega();
2879+
double pTid1 = 1 / std::abs(1.0 / pOrbital - 1.0 / pSpin1);
2880+
double pTid2 = 1 / std::abs(1.0 / pOrbital - 1.0 / pSpin2);
2881+
2882+
double tauConv1 = m_Star1->CalculateEddyTurnoverTimescale();
2883+
double tauConv2 = m_Star2->CalculateEddyTurnoverTimescale();
2884+
2885+
double fConv1 = std::min(1.0, std::pow(pTid1 / (2.0 * tauConv1), 2));
2886+
double fConv2 = std::min(1.0, std::pow(pTid2 / (2.0 * tauConv2), 2));
2887+
double kOverT1 = (2.0 / 21.0) * (fConv1 / tauConv1) * (envMass1 / m_Star1->Mass());
2888+
double kOverT2 = (2.0 / 21.0) * (fConv2 / tauConv2) * (envMass2 / m_Star2->Mass());
2889+
2890+
double R1_AU = m_Star1->Radius() * RSOL_TO_AU;
2891+
double R2_AU = m_Star2->Radius() * RSOL_TO_AU;
2892+
2893+
double omega1 = m_Star1->Omega();
2894+
double omega2 = m_Star2->Omega();
2895+
2896+
double E2_Dynamical_Zahn1 = 1.592 * std::pow(10, -9) * std::pow(m_Star1->Mass(), 2.84); // Eq. 43
2897+
double E2_Dynamical_Zahn2 = 1.592 * std::pow(10, -9) * std::pow(m_Star2->Mass(), 2.84); // Eq. 43
2898+
2899+
// double DSemiMajorAxis1Dt_Zahn = - m_SemiMajorAxis * 12.0 * kOverT1 * (q1 * (1.0 + q1)) * std::pow(R1_AU / m_SemiMajorAxis, 8) * ( (1.0 - omega1_over_omega) + (m_Eccentricity * m_Eccentricity * (23.0 - 27.0/2.0 * omega1_over_omega)) ) ; // Eq. (4.3)
2900+
// double DEccentricity1Dt_Zahn = -3.0/2.0 * kOverT1 * (q1 * (1.0 + q1)) * std::pow(R1_AU / m_SemiMajorAxis, 8) * m_Eccentricity * ( 27.0 - 33.0/2.0 * omega1_over_omega ); // Eq. (4.4)
2901+
// double DOmega1Dt_Zahn = (1.0 / m_Star1->CalculateMomentOfInertiaAU()) * 6.0 * kOverT1 * (q1 * q1) * std::pow(R1_AU / m_SemiMajorAxis, 6) * ((omega - omega1) + m_Eccentricity * m_Eccentricity * (26.0 * omega - 15.0 * omega1) ); // Eq. (4.5)
2902+
2903+
double w10 = 1.0 * omega;
2904+
double w12_1 = (1.0 * omega - 2.0 * omega1);
2905+
double w12_2 = (1.0 * omega - 2.0 * omega2);
2906+
2907+
double w22_1 = (2.0 * omega - 2.0 * omega1);
2908+
double w22_2 = (2.0 * omega - 2.0 * omega2);
2909+
2910+
double w32_1 = (3.0 * omega - 2.0 * omega1);
2911+
double w32_2 = (3.0 * omega - 2.0 * omega2);
2912+
2913+
double ImK22_Zahn_Equilibrium1 = 2.0 * kOverT1 * (R1_AU * R1_AU * R1_AU / (G_AU_Msol_yr * m_Star1->Mass())) * w22_1; // Eq. (4.1), (4.2)
2914+
ImK22_Zahn_Equilibrium1 = std::isnan(ImK22_Zahn_Equilibrium1) ? 0.0 : ImK22_Zahn_Equilibrium1; // check for NaN
2915+
double ImK22_Zahn_Equilibrium2 = 2.0 * kOverT2 * (R2_AU * R2_AU * R2_AU / (G_AU_Msol_yr * m_Star2->Mass())) * w22_2; // Eq. (4.1), (4.2)
2916+
ImK22_Zahn_Equilibrium2 = std::isnan(ImK22_Zahn_Equilibrium2) ? 0.0 : ImK22_Zahn_Equilibrium2; // check for NaN
2917+
2918+
double ImK10_Zahn_Equilibrium1 = 2.0 * kOverT1 * (R1_AU * R1_AU * R1_AU / (G_AU_Msol_yr * m_Star1->Mass())) * w10;
2919+
ImK10_Zahn_Equilibrium1 = std::isnan(ImK10_Zahn_Equilibrium1) ? 0.0 : ImK10_Zahn_Equilibrium1; // check for NaN
2920+
double ImK10_Zahn_Equilibrium2 = 2.0 * kOverT2 * (R2_AU * R2_AU * R2_AU / (G_AU_Msol_yr * m_Star2->Mass())) * w10;
2921+
ImK10_Zahn_Equilibrium2 = std::isnan(ImK10_Zahn_Equilibrium2) ? 0.0 : ImK10_Zahn_Equilibrium2; // check for NaN
2922+
2923+
double ImK12_Zahn_Equilibrium1 = 2.0 * kOverT1 * (R1_AU * R1_AU * R1_AU / (G_AU_Msol_yr * m_Star1->Mass())) * w12_1;
2924+
ImK12_Zahn_Equilibrium1 = std::isnan(ImK12_Zahn_Equilibrium1) ? 0.0 : ImK12_Zahn_Equilibrium1; // check for NaN
2925+
double ImK12_Zahn_Equilibrium2 = 2.0 * kOverT2 * (R2_AU * R2_AU * R2_AU / (G_AU_Msol_yr * m_Star2->Mass())) * w12_2;
2926+
ImK12_Zahn_Equilibrium2 = std::isnan(ImK12_Zahn_Equilibrium2) ? 0.0 : ImK12_Zahn_Equilibrium2; // check for NaN
2927+
2928+
double ImK32_Zahn_Equilibrium1 = 2.0 * kOverT1 * (R1_AU * R1_AU * R1_AU / (G_AU_Msol_yr * m_Star1->Mass())) * w32_1;
2929+
ImK32_Zahn_Equilibrium1 = std::isnan(ImK32_Zahn_Equilibrium1) ? 0.0 : ImK32_Zahn_Equilibrium1; // check for NaN
2930+
double ImK32_Zahn_Equilibrium2 = 2.0 * kOverT2 * (R2_AU * R2_AU * R2_AU / (G_AU_Msol_yr * m_Star2->Mass())) * w32_2;
2931+
ImK32_Zahn_Equilibrium2 = std::isnan(ImK32_Zahn_Equilibrium2) ? 0.0 : ImK32_Zahn_Equilibrium2; // check for NaN
2932+
2933+
2934+
double ImK22_Zahn_Dynamical1 = std::copysign(1.0, w22_1) * E2_Dynamical_Zahn1 * std::pow(std::sqrt(R1_AU * R1_AU * R1_AU / (G_AU_Msol_yr * m_Star1->Mass())) * std::abs(w22_1), 8.0/3.0); // Eq. (5.5)
2935+
ImK22_Zahn_Dynamical1 = std::isnan(ImK22_Zahn_Dynamical1) ? 0.0 : ImK22_Zahn_Dynamical1; // check for NaN
2936+
double ImK22_Zahn_Dynamical2 = std::copysign(1.0, w22_2) * E2_Dynamical_Zahn2 * std::pow(std::sqrt(R2_AU * R2_AU * R2_AU / (G_AU_Msol_yr * m_Star2->Mass())) * std::abs(w22_2), 8.0/3.0); // Eq. (5.5)
2937+
ImK22_Zahn_Dynamical2 = std::isnan(ImK22_Zahn_Dynamical2) ? 0.0 : ImK22_Zahn_Dynamical2; // check for NaN
2938+
2939+
double ImK10_Zahn_Dynamical1 = std::copysign(1.0, w10) * E2_Dynamical_Zahn1 * std::pow(std::sqrt(R1_AU * R1_AU * R1_AU / (G_AU_Msol_yr * m_Star1->Mass())) * std::abs(w10), 8.0/3.0); // Eq. (5.5)
2940+
ImK10_Zahn_Dynamical1 = std::isnan(ImK10_Zahn_Dynamical1) ? 0.0 : ImK10_Zahn_Dynamical1; // check for NaN
2941+
double ImK10_Zahn_Dynamical2 =std::copysign(1.0, w22_1) * E2_Dynamical_Zahn2 * std::pow(std::sqrt(R2_AU * R2_AU * R2_AU / (G_AU_Msol_yr * m_Star2->Mass())) * std::abs(w10), 8.0/3.0); // Eq. (5.5)
2942+
ImK10_Zahn_Dynamical2 = std::isnan(ImK10_Zahn_Dynamical2) ? 0.0 : ImK10_Zahn_Dynamical2; // check for NaN
2943+
2944+
double ImK12_Zahn_Dynamical1 = std::copysign(1.0, w12_1) * E2_Dynamical_Zahn1 * std::pow(std::sqrt(R1_AU * R1_AU * R1_AU / (G_AU_Msol_yr * m_Star1->Mass())) * std::abs(w12_1), 8.0/3.0); // Eq. (5.5)
2945+
ImK12_Zahn_Dynamical1 = std::isnan(ImK12_Zahn_Dynamical1) ? 0.0 : ImK12_Zahn_Dynamical1; // check for NaN
2946+
double ImK12_Zahn_Dynamical2 = std::copysign(1.0, w12_2) * E2_Dynamical_Zahn2 * std::pow(std::sqrt(R2_AU * R2_AU * R2_AU / (G_AU_Msol_yr * m_Star2->Mass())) * std::abs(w12_2), 8.0/3.0); // Eq. (5.5)
2947+
ImK12_Zahn_Dynamical2 = std::isnan(ImK12_Zahn_Dynamical2) ? 0.0 : ImK12_Zahn_Dynamical2; // check for NaN
2948+
2949+
double ImK32_Zahn_Dynamical1 = std::copysign(1.0, w32_1) * E2_Dynamical_Zahn1 * std::pow(std::sqrt(R1_AU * R1_AU * R1_AU / (G_AU_Msol_yr * m_Star1->Mass())) * std::abs(w32_1), 8.0/3.0); // Eq. (5.5)
2950+
ImK32_Zahn_Dynamical1 = std::isnan(ImK32_Zahn_Dynamical1) ? 0.0 : ImK32_Zahn_Dynamical1; // check for NaN
2951+
double ImK32_Zahn_Dynamical2 = std::copysign(1.0, w32_2) * E2_Dynamical_Zahn2 * std::pow(std::sqrt(R2_AU * R2_AU * R2_AU / (G_AU_Msol_yr * m_Star2->Mass())) * std::abs(w32_2), 8.0/3.0); // Eq. (5.5)
2952+
ImK32_Zahn_Dynamical2 = std::isnan(ImK32_Zahn_Dynamical2) ? 0.0 : ImK32_Zahn_Dynamical2; // check for NaN
2953+
2954+
DBL_DBL_DBL_DBL ImKnm1_tidal = std::make_tuple(ImK10_Zahn_Equilibrium1+ImK10_Zahn_Dynamical1, ImK12_Zahn_Equilibrium1+ImK12_Zahn_Dynamical1, ImK22_Zahn_Equilibrium1 + ImK22_Zahn_Dynamical1, ImK32_Zahn_Equilibrium1+ImK32_Zahn_Dynamical1);
2955+
DBL_DBL_DBL_DBL ImKnm2_tidal = std::make_tuple(ImK10_Zahn_Equilibrium2+ImK10_Zahn_Dynamical2, ImK12_Zahn_Equilibrium2+ImK12_Zahn_Dynamical2, ImK22_Zahn_Equilibrium2 + ImK22_Zahn_Dynamical2, ImK32_Zahn_Equilibrium2 + ImK32_Zahn_Dynamical2);
2956+
2957+
2958+
double DSemiMajorAxis1Dt_tidal = CalculateDSemiMajorAxisTidalDt(ImKnm1_tidal, m_Star1); // change in semi-major axis from star1
2959+
double DSemiMajorAxis2Dt_tidal = CalculateDSemiMajorAxisTidalDt(ImKnm2_tidal, m_Star2); // change in semi-major axis from star2
2960+
2961+
double DEccentricity1Dt_tidal = CalculateDEccentricityTidalDt(ImKnm1_tidal, m_Star1); // change in eccentricity from star1
2962+
double DEccentricity2Dt_tidal = CalculateDEccentricityTidalDt(ImKnm2_tidal, m_Star2); // change in eccentricity from star2
2963+
2964+
double DOmega1Dt_tidal = CalculateDOmegaTidalDt(ImKnm1_tidal, m_Star1); // change in spin from star1
2965+
double DOmega2Dt_tidal = CalculateDOmegaTidalDt(ImKnm2_tidal, m_Star2); // change in spin from star2
2966+
2967+
// limit change in stellar and orbital properties from tides to a maximum fraction of the current value
2968+
double fraction_tidal_change = 1.0;
2969+
fraction_tidal_change = std::min(fraction_tidal_change, std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * OrbitalAngularVelocity() / (DOmega1Dt_tidal * p_Dt * MYR_TO_YEAR)));
2970+
fraction_tidal_change = std::min(fraction_tidal_change, std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * OrbitalAngularVelocity() / (DOmega2Dt_tidal * p_Dt * MYR_TO_YEAR)));
2971+
fraction_tidal_change = std::min(fraction_tidal_change, std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * m_SemiMajorAxis / ((DSemiMajorAxis1Dt_tidal + DSemiMajorAxis2Dt_tidal) * p_Dt * MYR_TO_YEAR)));
2972+
fraction_tidal_change = std::min(fraction_tidal_change, std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * m_Eccentricity / ((DEccentricity1Dt_tidal + DEccentricity2Dt_tidal) * p_Dt * MYR_TO_YEAR)));
2973+
2974+
m_Star1->SetOmega(m_Star1->Omega() + fraction_tidal_change * (DOmega1Dt_tidal * p_Dt * MYR_TO_YEAR)); // evolve star 1 spin
2975+
m_Star2->SetOmega(m_Star2->Omega() + fraction_tidal_change * (DOmega2Dt_tidal * p_Dt * MYR_TO_YEAR)); // evolve star 2 spin
2976+
m_SemiMajorAxis = m_SemiMajorAxis + fraction_tidal_change * ((DSemiMajorAxis1Dt_tidal + DSemiMajorAxis2Dt_tidal) * p_Dt * MYR_TO_YEAR); // evolve separation
2977+
m_Eccentricity = m_Eccentricity + fraction_tidal_change * ((DEccentricity1Dt_tidal + DEccentricity2Dt_tidal) * p_Dt * MYR_TO_YEAR); // evolve eccentricity
2978+
2979+
m_CircularizationTimescale = - m_Eccentricity / (DEccentricity1Dt_tidal + DEccentricity2Dt_tidal) * YEAR_TO_MYR; // Circularization timescale in Myr (for output files)
2980+
m_CircularizationTimescale = (std::isnan(m_CircularizationTimescale) || std::isinf(m_CircularizationTimescale))? 0.0 : m_CircularizationTimescale; // check for NaN or Inf for circular binaries
2981+
2982+
m_SynchronizationTimescale1 = - (m_Star1->Omega() - omega) / DOmega1Dt_tidal * YEAR_TO_MYR; // Synchronization timescale for Star1 in Myr (for output files)
2983+
m_SynchronizationTimescale1 = (std::isnan(m_SynchronizationTimescale1) || std::isinf(m_SynchronizationTimescale1))? 0.0 : m_SynchronizationTimescale1; // check for NaN or Inf for synchronized binaries
2984+
2985+
m_SynchronizationTimescale2 = - (m_Star2->Omega() - omega) / DOmega2Dt_tidal * YEAR_TO_MYR; // Synchronization timescale for Star2 in Myr (for output files)
2986+
m_SynchronizationTimescale2 = (std::isnan(m_SynchronizationTimescale2) || std::isinf(m_SynchronizationTimescale2))? 0.0 : m_SynchronizationTimescale2; // check for NaN or Inf for synchronized binaries
2987+
2988+
m_TotalAngularMomentum = CalculateAngularMomentum(); // re-calculate angular momenta
2989+
m_OrbitalAngularMomentum = CalculateOrbitalAngularMomentum(m_Star1->Mass(), m_Star2->Mass(), m_SemiMajorAxis, m_Eccentricity);
2990+
2991+
} break;
2992+
28602993
default:
28612994
// the only way this can happen is if someone added a TIDES_PRESCRIPTION
28622995
// and it isn't accounted for in this code. We should not default here, with or without a warning.
@@ -3179,9 +3312,8 @@ void BaseBinaryStar::EvaluateBinary(const double p_Dt) {
31793312
}
31803313

31813314
CalculateEnergyAndAngularMomentum(); // perform energy and angular momentum calculations
3182-
3183-
ProcessTides(p_Dt); // process tides if required
3184-
3315+
ProcessTides(p_Dt);
3316+
// process tides if required
31853317
// assign new values to "previous" values, for following timestep
31863318
m_EccentricityPrev = m_Eccentricity;
31873319
m_SemiMajorAxisPrev = m_SemiMajorAxis;

‎src/typedefs.h‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -937,11 +937,12 @@ const COMPASUnorderedMap<STELLAR_POPULATION, std::string> STELLAR_POPULATION_LAB
937937
};
938938

939939
// tides prescriptions
940-
enum class TIDES_PRESCRIPTION: int { NONE, PERFECT, KAPIL2025 };
940+
enum class TIDES_PRESCRIPTION: int { NONE, PERFECT, KAPIL2025, ZAHN1977 };
941941
const COMPASUnorderedMap<TIDES_PRESCRIPTION, std::string> TIDES_PRESCRIPTION_LABEL = {
942942
{ TIDES_PRESCRIPTION::NONE, "NONE" },
943943
{ TIDES_PRESCRIPTION::PERFECT, "PERFECT" },
944-
{ TIDES_PRESCRIPTION::KAPIL2025, "KAPIL2025" }
944+
{ TIDES_PRESCRIPTION::KAPIL2025, "KAPIL2025" },
945+
{ TIDES_PRESCRIPTION::ZAHN1977, "ZAHN1977" }
945946
};
946947

947948
// symbolic names for timescales

0 commit comments

Comments
 (0)