Skip to content

Commit

Permalink
Merge pull request #1025 from jeffriley/issue#1022
Browse files Browse the repository at this point in the history
Fix for issue #1022
  • Loading branch information
ilyamandel authored Dec 11, 2023
2 parents 9fa0baf + 501ae24 commit 200b1a1
Show file tree
Hide file tree
Showing 3 changed files with 24 additions and 21 deletions.
2 changes: 1 addition & 1 deletion src/BaseBinaryStar.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1893,7 +1893,7 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) {
// NOTE: Critical mass ratio is defined as mAccretor/mDonor
double qCrit = m_Donor->CalculateCriticalMassRatio(m_Accretor->IsDegenerate());

isUnstable = (m_Accretor->Mass()/m_Donor->Mass()) < qCrit;
isUnstable = (m_Accretor->Mass() / m_Donor->Mass()) < qCrit;
m_FractionAccreted = 1.0; // Accretion is assumed fully conservative for qCrit calculations
}
else { // Determine stability based on zetas
Expand Down
38 changes: 19 additions & 19 deletions src/BaseStar.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1364,11 +1364,11 @@ double BaseStar::CalculateCriticalMassRatio(const bool p_AccretorIsDegenerate) {
* from Ge et al. (2020), GE20 or GE20_IC. The first is the full adiabatic response, the second assumes
* artificially isentropic envelopes.
*
* double BaseStar::InterpolateGe20QCrit( const QCRIT_PRESCRIPTION p_qCritPrescription)
* double BaseStar::InterpolateGe20QCrit(const QCRIT_PRESCRIPTION p_qCritPrescription)
*
* @return Interpolated value of either the critical mass ratio or zeta for given stellar mass / radius
*/
double BaseStar::InterpolateGe20QCrit( const QCRIT_PRESCRIPTION p_qCritPrescription) {
double BaseStar::InterpolateGe20QCrit(const QCRIT_PRESCRIPTION p_qCritPrescription) {

// Get vector of masses from GE20_QCRIT_AND_ZETA
std::vector<double> massesFromGe20 = std::get<0>(GE20_QCRIT_AND_ZETA);
Expand All @@ -1384,8 +1384,8 @@ double BaseStar::InterpolateGe20QCrit( const QCRIT_PRESCRIPTION p_qCritPrescript
upperMassInd = 0;
}
else if (upperMassInd == -1) {
lowerMassInd = massesFromGe20.size();
upperMassInd = massesFromGe20.size();
lowerMassInd = massesFromGe20.size() - 1;
upperMassInd = massesFromGe20.size() - 1;
}

// Get vector of radii from GE20_QCRIT_AND_ZETA for the lower and upper mass indices
Expand All @@ -1398,12 +1398,12 @@ double BaseStar::InterpolateGe20QCrit( const QCRIT_PRESCRIPTION p_qCritPrescript

// One of the following must be set
if (p_qCritPrescription == QCRIT_PRESCRIPTION::GE20) {
qCritVectorLowerMass = std::get<1>(radiiQCritsZetasFromGe20[lowerMassInd]);
qCritVectorUpperMass = std::get<1>(radiiQCritsZetasFromGe20[upperMassInd]);
qCritVectorLowerMass = std::get<1>(radiiQCritsZetasFromGe20[lowerMassInd]);
qCritVectorUpperMass = std::get<1>(radiiQCritsZetasFromGe20[upperMassInd]);
}
else if (p_qCritPrescription == QCRIT_PRESCRIPTION::GE20_IC) {
qCritVectorLowerMass = std::get<2>(radiiQCritsZetasFromGe20[lowerMassInd]);
qCritVectorUpperMass = std::get<2>(radiiQCritsZetasFromGe20[upperMassInd]);
qCritVectorLowerMass = std::get<2>(radiiQCritsZetasFromGe20[lowerMassInd]);
qCritVectorUpperMass = std::get<2>(radiiQCritsZetasFromGe20[upperMassInd]);
}

// Get vector of radii from GE20_QCRIT_AND_ZETA for both lower and upper masses
Expand All @@ -1416,8 +1416,8 @@ double BaseStar::InterpolateGe20QCrit( const QCRIT_PRESCRIPTION p_qCritPrescript
upperRadiusLowerMassInd = 0;
}
else if (upperRadiusLowerMassInd == -1) {
lowerRadiusLowerMassInd = logRadiusVectorLowerMass.size();
upperRadiusLowerMassInd = logRadiusVectorLowerMass.size();
lowerRadiusLowerMassInd = logRadiusVectorLowerMass.size() - 1;
upperRadiusLowerMassInd = logRadiusVectorLowerMass.size() - 1;
}

std::vector<int> indR1 = utils::binarySearch(logRadiusVectorUpperMass, log10(m_Radius));
Expand All @@ -1429,8 +1429,8 @@ double BaseStar::InterpolateGe20QCrit( const QCRIT_PRESCRIPTION p_qCritPrescript
upperRadiusUpperMassInd = 0;
}
else if (upperRadiusUpperMassInd == -1) {
lowerRadiusUpperMassInd = logRadiusVectorUpperMass.size();
upperRadiusUpperMassInd = logRadiusVectorUpperMass.size();
lowerRadiusUpperMassInd = logRadiusVectorUpperMass.size() - 1;
upperRadiusUpperMassInd = logRadiusVectorUpperMass.size() - 1;
}

// Set the 4 boundary points for the 2D interpolation
Expand All @@ -1441,18 +1441,18 @@ double BaseStar::InterpolateGe20QCrit( const QCRIT_PRESCRIPTION p_qCritPrescript

double lowerMass = massesFromGe20[lowerMassInd];
double upperMass = massesFromGe20[upperMassInd];
double lowerRadiusLowerMass = pow(10, logRadiusVectorLowerMass[lowerRadiusLowerMassInd]);
double upperRadiusLowerMass = pow(10, logRadiusVectorLowerMass[upperRadiusLowerMassInd]);
double lowerRadiusUpperMass = pow(10, logRadiusVectorUpperMass[lowerRadiusUpperMassInd]);
double upperRadiusUpperMass = pow(10, logRadiusVectorUpperMass[upperRadiusUpperMassInd]);

double lowerRadiusLowerMass = PPOW(10.0, logRadiusVectorLowerMass[lowerRadiusLowerMassInd]);
double upperRadiusLowerMass = PPOW(10.0, logRadiusVectorLowerMass[upperRadiusLowerMassInd]);
double lowerRadiusUpperMass = PPOW(10.0, logRadiusVectorUpperMass[lowerRadiusUpperMassInd]);
double upperRadiusUpperMass = PPOW(10.0, logRadiusVectorUpperMass[upperRadiusUpperMassInd]);

// Interpolate on the radii first, then the masses
double qCritLowerMass = qLowLow + (upperRadiusLowerMass - m_Radius)/(upperRadiusLowerMass - lowerRadiusLowerMass) * (qLowUpp - qLowLow);
double qCritUpperMass = qUppLow + (upperRadiusUpperMass - m_Radius)/(upperRadiusUpperMass - lowerRadiusUpperMass) * (qUppUpp - qUppLow);
double qCritLowerMass = qLowLow + (upperRadiusLowerMass - m_Radius) / (upperRadiusLowerMass - lowerRadiusLowerMass) * (qLowUpp - qLowLow);
double qCritUpperMass = qUppLow + (upperRadiusUpperMass - m_Radius) / (upperRadiusUpperMass - lowerRadiusUpperMass) * (qUppUpp - qUppLow);
double interpolatedQCrit = qCritLowerMass + (upperMass - m_Mass)/(upperMass - lowerMass) * (qCritUpperMass - qCritLowerMass);

return interpolatedQCrit;

}


Expand Down
5 changes: 4 additions & 1 deletion src/changelog.h
Original file line number Diff line number Diff line change
Expand Up @@ -1071,8 +1071,11 @@
// - Fixed CalculateOrbitalAngularMomentum() (now uses eccentricity)
// - Added links to online documentation to splash string
// - Constants 'G1' and 'G_SN' renamed to 'G_AU_Msol_yr' and 'G_km_Msol_s' respectively
// 02.41.01 JR - Dec 11, 2023 - Defect repair, a little code cleanup:
// - Fix for issue #1022 - incorrect index used for last array entry.
// - A little code cleanup

const std::string VERSION_STRING = "02.41.00";
const std::string VERSION_STRING = "02.41.01";


# endif // __changelog_h__

0 comments on commit 200b1a1

Please sign in to comment.