diff --git a/src/BH.h b/src/BH.h index b597c966a..ea99dd1d5 100755 --- a/src/BH.h +++ b/src/BH.h @@ -56,7 +56,7 @@ class BH: virtual public BaseStar, public Remnants { double CalculateMassLossRate() { return 0.0; } // Ensure that BHs don't lose mass in winds - double CalculateMomentOfInertia(const double p_RemnantRadius = 0.0) const { return (2.0 / 5.0) * m_Mass * m_Radius * m_Radius; } + double CalculateMomentOfInertia(const double p_RemnantRadius = 0.0) const { return -(2.0 / 5.0) * m_Mass * m_Radius * m_Radius; } double CalculateMomentOfInertiaAU(const double p_RemnantRadius = 0.0) const { return CalculateMomentOfInertia(p_RemnantRadius * RSOL_TO_AU) * RSOL_TO_AU * RSOL_TO_AU; } double CalculateRadiusOnPhase() const { return CalculateRadiusOnPhase_Static(m_Mass); } // Use class member variables - returns radius in Rsol diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 103f5ead3..525fed3ab 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -161,8 +161,8 @@ BaseBinaryStar::BaseBinaryStar(const unsigned long int p_Seed, const long int p_ double starToRocheLobeRadiusRatio1 = (m_Star1->Radius() * RSOL_TO_AU) / (m_SemiMajorAxis * (1.0 - m_Eccentricity) * CalculateRocheLobeRadius_Static(mass1, mass2)); double starToRocheLobeRadiusRatio2 = (m_Star2->Radius() * RSOL_TO_AU) / (m_SemiMajorAxis * (1.0 - m_Eccentricity) * CalculateRocheLobeRadius_Static(mass2, mass1)); - m_Flags.massesEquilibrated = false; // default - m_Flags.massesEquilibratedAtBirth = false; // default + m_Flags.massesEquilibrated = false; // default + m_Flags.massesEquilibratedAtBirth = false; // default rlof = utils::Compare(starToRocheLobeRadiusRatio1, 1.0) > 0 || utils::Compare(starToRocheLobeRadiusRatio2, 1.0) > 0; // either star overflowing Roche Lobe? @@ -272,9 +272,19 @@ void BaseBinaryStar::SetRemainingValues() { m_SemiMajorAxisAtDCOFormation = DEFAULT_INITIAL_DOUBLE_VALUE; m_EccentricityAtDCOFormation = DEFAULT_INITIAL_DOUBLE_VALUE; - // if CHE enabled, update rotational frequency for constituent stars - assume tidally locked + double gyrationRadius1 = m_Star1->CalculateGyrationRadius(); + double gyrationRadius2 = m_Star2->CalculateGyrationRadius(); + + m_TotalEnergy = CalculateTotalEnergy(m_SemiMajorAxis, m_Star1->Mass(), m_Star2->Mass(), m_Star1->RZAMS(), m_Star2->RZAMS(), m_Star1->Omega(), m_Star2->Omega(), gyrationRadius1, gyrationRadius2); + + m_TotalAngularMomentum = CalculateAngularMomentum(m_SemiMajorAxis, m_Eccentricity, m_Star1->Mass(), m_Star2->Mass(), m_Star1->RZAMS(), m_Star2->RZAMS(), m_Star1->Omega(), m_Star2->Omega(), gyrationRadius1, gyrationRadius2); + m_TotalAngularMomentumPrev = m_TotalAngularMomentum; + + m_Omega = 0.0; - if (OPTIONS->CHEMode() != CHE_MODE::NONE) { + if (OPTIONS->CHEMode() != CHE_MODE::NONE) { // CHE enabled? + + // CHE enabled, update rotational frequency for constituent stars - assume tidally locked m_Star1->SetOmega(OrbitalAngularVelocity()); m_Star2->SetOmega(OrbitalAngularVelocity()); @@ -312,209 +322,185 @@ void BaseBinaryStar::SetRemainingValues() { } } - double gyrationRadius1 = m_Star1->CalculateGyrationRadius(); - double gyrationRadius2 = m_Star2->CalculateGyrationRadius(); - - m_TotalEnergy = CalculateTotalEnergy(m_SemiMajorAxis, - m_Star1->Mass(), - m_Star2->Mass(), - m_Star1->RZAMS(), - m_Star2->RZAMS(), - m_Star1->Omega(), - m_Star2->Omega(), - gyrationRadius1, - gyrationRadius2); - - m_TotalAngularMomentum = CalculateAngularMomentum(m_SemiMajorAxis, - m_Eccentricity, - m_Star1->Mass(), - m_Star2->Mass(), - m_Star1->RZAMS(), - m_Star2->RZAMS(), - m_Star1->Omega(), - m_Star2->Omega(), - gyrationRadius1, - gyrationRadius2); + double totalMass = m_Star1->Mass() + m_Star2->Mass(); + double reducedMass = (m_Star1->Mass() * m_Star2->Mass()) / totalMass; + m_OrbitalEnergy = CalculateOrbitalEnergy(reducedMass, totalMass, m_SemiMajorAxis); + m_OrbitalEnergyPrev = m_OrbitalEnergy; - double totalMass = m_Star1->Mass() + m_Star2->Mass(); - double reducedMass = (m_Star1->Mass() * m_Star2->Mass()) / totalMass; - m_OrbitalEnergy = CalculateOrbitalEnergy(reducedMass, totalMass, m_SemiMajorAxis); - m_OrbitalEnergyPrev = m_OrbitalEnergy; + m_OrbitalAngularMomentum = CalculateOrbitalAngularMomentum(reducedMass, totalMass, m_SemiMajorAxis); + m_OrbitalAngularMomentumPrev = m_OrbitalAngularMomentum; - m_OrbitalAngularMomentum = CalculateOrbitalAngularMomentum(reducedMass, totalMass, m_SemiMajorAxis); - m_OrbitalAngularMomentumPrev = m_OrbitalAngularMomentum; + m_Time = DEFAULT_INITIAL_DOUBLE_VALUE; + m_Dt = DEFAULT_INITIAL_DOUBLE_VALUE; + m_TimePrev = DEFAULT_INITIAL_DOUBLE_VALUE; + m_DCOFormationTime = DEFAULT_INITIAL_DOUBLE_VALUE; - m_Time = DEFAULT_INITIAL_DOUBLE_VALUE; - m_Dt = DEFAULT_INITIAL_DOUBLE_VALUE; - m_TimePrev = DEFAULT_INITIAL_DOUBLE_VALUE; - m_DCOFormationTime = DEFAULT_INITIAL_DOUBLE_VALUE; + m_aMassLossDiff = DEFAULT_INITIAL_DOUBLE_VALUE; - m_aMassLossDiff = DEFAULT_INITIAL_DOUBLE_VALUE; + m_aMassTransferDiff = DEFAULT_INITIAL_DOUBLE_VALUE; - m_aMassTransferDiff = DEFAULT_INITIAL_DOUBLE_VALUE; + m_MassTransferTrackerHistory = MT_TRACKING::NO_MASS_TRANSFER; + m_MassTransfer = false; - m_MassTransferTrackerHistory = MT_TRACKING::NO_MASS_TRANSFER; - m_MassTransfer = false; + m_JLoss = OPTIONS->MassTransferJloss(); - m_JLoss = OPTIONS->MassTransferJloss(); - - m_FractionAccreted = OPTIONS->MassTransferFractionAccreted(); + m_FractionAccreted = OPTIONS->MassTransferFractionAccreted(); // Common Envelope - m_CEDetails.CEEcount = 0; - m_CEDetails.CEEnow = false; - m_CEDetails.doubleCoreCE = false; - m_CEDetails.optimisticCE = false; - m_CEDetails.postCEE.eccentricity = DEFAULT_INITIAL_DOUBLE_VALUE; - m_CEDetails.postCEE.rocheLobe1to2 = DEFAULT_INITIAL_DOUBLE_VALUE; - m_CEDetails.postCEE.rocheLobe2to1 = DEFAULT_INITIAL_DOUBLE_VALUE; - m_CEDetails.postCEE.semiMajorAxis = DEFAULT_INITIAL_DOUBLE_VALUE; - m_CEDetails.preCEE.eccentricity = DEFAULT_INITIAL_DOUBLE_VALUE; - m_CEDetails.preCEE.rocheLobe1to2 = DEFAULT_INITIAL_DOUBLE_VALUE; - m_CEDetails.preCEE.rocheLobe2to1 = DEFAULT_INITIAL_DOUBLE_VALUE; - m_CEDetails.preCEE.semiMajorAxis = DEFAULT_INITIAL_DOUBLE_VALUE; - - m_Flags.stellarMerger = false; - m_Flags.stellarMergerAtBirth = false; - - m_Mass1Final = DEFAULT_INITIAL_DOUBLE_VALUE; - m_Mass2Final = DEFAULT_INITIAL_DOUBLE_VALUE; - m_MassEnv1 = DEFAULT_INITIAL_DOUBLE_VALUE; - m_MassEnv2 = DEFAULT_INITIAL_DOUBLE_VALUE; - - m_ZetaLobe = DEFAULT_INITIAL_DOUBLE_VALUE; - m_ZetaStar = DEFAULT_INITIAL_DOUBLE_VALUE; + m_CEDetails.CEEcount = 0; + m_CEDetails.CEEnow = false; + m_CEDetails.doubleCoreCE = false; + m_CEDetails.optimisticCE = false; + m_CEDetails.postCEE.eccentricity = DEFAULT_INITIAL_DOUBLE_VALUE; + m_CEDetails.postCEE.rocheLobe1to2 = DEFAULT_INITIAL_DOUBLE_VALUE; + m_CEDetails.postCEE.rocheLobe2to1 = DEFAULT_INITIAL_DOUBLE_VALUE; + m_CEDetails.postCEE.semiMajorAxis = DEFAULT_INITIAL_DOUBLE_VALUE; + m_CEDetails.preCEE.eccentricity = DEFAULT_INITIAL_DOUBLE_VALUE; + m_CEDetails.preCEE.rocheLobe1to2 = DEFAULT_INITIAL_DOUBLE_VALUE; + m_CEDetails.preCEE.rocheLobe2to1 = DEFAULT_INITIAL_DOUBLE_VALUE; + m_CEDetails.preCEE.semiMajorAxis = DEFAULT_INITIAL_DOUBLE_VALUE; + + m_Flags.stellarMerger = false; + m_Flags.stellarMergerAtBirth = false; + + m_Mass1Final = DEFAULT_INITIAL_DOUBLE_VALUE; + m_Mass2Final = DEFAULT_INITIAL_DOUBLE_VALUE; + m_MassEnv1 = DEFAULT_INITIAL_DOUBLE_VALUE; + m_MassEnv2 = DEFAULT_INITIAL_DOUBLE_VALUE; + + m_ZetaLobe = DEFAULT_INITIAL_DOUBLE_VALUE; + m_ZetaStar = DEFAULT_INITIAL_DOUBLE_VALUE; // Initialise other parameters to 0 - m_uK = DEFAULT_INITIAL_DOUBLE_VALUE; - m_CosIPrime = DEFAULT_INITIAL_DOUBLE_VALUE; - m_IPrime = DEFAULT_INITIAL_DOUBLE_VALUE; - m_TimeToCoalescence = DEFAULT_INITIAL_DOUBLE_VALUE; + m_uK = DEFAULT_INITIAL_DOUBLE_VALUE; + m_CosIPrime = DEFAULT_INITIAL_DOUBLE_VALUE; + m_IPrime = DEFAULT_INITIAL_DOUBLE_VALUE; + m_TimeToCoalescence = DEFAULT_INITIAL_DOUBLE_VALUE; - m_SupernovaState = SN_STATE::NONE; + m_SupernovaState = SN_STATE::NONE; - m_Flags.mergesInHubbleTime = false; - m_Unbound = false; + m_Flags.mergesInHubbleTime = false; + m_Unbound = false; - m_SystemicVelocity = Vector3d(); - m_OrbitalAngularMomentumVector = Vector3d(); - m_ThetaE = DEFAULT_INITIAL_DOUBLE_VALUE; - m_PhiE = DEFAULT_INITIAL_DOUBLE_VALUE; - m_PsiE = DEFAULT_INITIAL_DOUBLE_VALUE; + m_SystemicVelocity = Vector3d(); + m_OrbitalAngularMomentumVector = Vector3d(); + m_ThetaE = DEFAULT_INITIAL_DOUBLE_VALUE; + m_PhiE = DEFAULT_INITIAL_DOUBLE_VALUE; + m_PsiE = DEFAULT_INITIAL_DOUBLE_VALUE; - m_SynchronizationTimescale = DEFAULT_INITIAL_DOUBLE_VALUE; - m_CircularizationTimescale = DEFAULT_INITIAL_DOUBLE_VALUE; + m_SynchronizationTimescale = DEFAULT_INITIAL_DOUBLE_VALUE; + m_CircularizationTimescale = DEFAULT_INITIAL_DOUBLE_VALUE; // RLOF details - m_RLOFDetails.experiencedRLOF = false; - m_RLOFDetails.immediateRLOFPostCEE = false; - m_RLOFDetails.isRLOF = false; - m_RLOFDetails.simultaneousRLOF = false; - m_RLOFDetails.stableRLOFPostCEE = false; + m_RLOFDetails.experiencedRLOF = false; + m_RLOFDetails.immediateRLOFPostCEE = false; + m_RLOFDetails.isRLOF = false; + m_RLOFDetails.simultaneousRLOF = false; + m_RLOFDetails.stableRLOFPostCEE = false; // RLOF details - properties 1 - m_RLOFDetails.props1.id = -1l; + m_RLOFDetails.props1.id = -1l; - m_RLOFDetails.props1.stellarType1 = STELLAR_TYPE::NONE; - m_RLOFDetails.props1.stellarType2 = STELLAR_TYPE::NONE; + m_RLOFDetails.props1.stellarType1 = STELLAR_TYPE::NONE; + m_RLOFDetails.props1.stellarType2 = STELLAR_TYPE::NONE; - m_RLOFDetails.props1.mass1 = DEFAULT_INITIAL_DOUBLE_VALUE; - m_RLOFDetails.props1.mass2 = DEFAULT_INITIAL_DOUBLE_VALUE; + m_RLOFDetails.props1.mass1 = DEFAULT_INITIAL_DOUBLE_VALUE; + m_RLOFDetails.props1.mass2 = DEFAULT_INITIAL_DOUBLE_VALUE; - m_RLOFDetails.props1.radius1 = DEFAULT_INITIAL_DOUBLE_VALUE; - m_RLOFDetails.props1.radius2 = DEFAULT_INITIAL_DOUBLE_VALUE; + m_RLOFDetails.props1.radius1 = DEFAULT_INITIAL_DOUBLE_VALUE; + m_RLOFDetails.props1.radius2 = DEFAULT_INITIAL_DOUBLE_VALUE; - m_RLOFDetails.props1.starToRocheLobeRadiusRatio1 = DEFAULT_INITIAL_DOUBLE_VALUE; - m_RLOFDetails.props1.starToRocheLobeRadiusRatio2 = DEFAULT_INITIAL_DOUBLE_VALUE; + m_RLOFDetails.props1.starToRocheLobeRadiusRatio1 = DEFAULT_INITIAL_DOUBLE_VALUE; + m_RLOFDetails.props1.starToRocheLobeRadiusRatio2 = DEFAULT_INITIAL_DOUBLE_VALUE; - m_RLOFDetails.props1.semiMajorAxis = DEFAULT_INITIAL_DOUBLE_VALUE; - m_RLOFDetails.props1.eccentricity = DEFAULT_INITIAL_DOUBLE_VALUE; + m_RLOFDetails.props1.semiMajorAxis = DEFAULT_INITIAL_DOUBLE_VALUE; + m_RLOFDetails.props1.eccentricity = DEFAULT_INITIAL_DOUBLE_VALUE; - m_RLOFDetails.props1.eventCounter = DEFAULT_INITIAL_ULONGINT_VALUE; + m_RLOFDetails.props1.eventCounter = DEFAULT_INITIAL_ULONGINT_VALUE; - m_RLOFDetails.props1.time = DEFAULT_INITIAL_DOUBLE_VALUE; + m_RLOFDetails.props1.time = DEFAULT_INITIAL_DOUBLE_VALUE; - m_RLOFDetails.props1.isRLOF1 = false; - m_RLOFDetails.props1.isRLOF2 = false; + m_RLOFDetails.props1.isRLOF1 = false; + m_RLOFDetails.props1.isRLOF2 = false; - m_RLOFDetails.props1.isCE = false; + m_RLOFDetails.props1.isCE = false; // RLOF details - properties 2 m_RLOFDetails.props2.id = -1l; - m_RLOFDetails.props2.stellarType1 = STELLAR_TYPE::NONE; - m_RLOFDetails.props2.stellarType2 = STELLAR_TYPE::NONE; + m_RLOFDetails.props2.stellarType1 = STELLAR_TYPE::NONE; + m_RLOFDetails.props2.stellarType2 = STELLAR_TYPE::NONE; - m_RLOFDetails.props2.mass1 = DEFAULT_INITIAL_DOUBLE_VALUE; - m_RLOFDetails.props2.mass2 = DEFAULT_INITIAL_DOUBLE_VALUE; + m_RLOFDetails.props2.mass1 = DEFAULT_INITIAL_DOUBLE_VALUE; + m_RLOFDetails.props2.mass2 = DEFAULT_INITIAL_DOUBLE_VALUE; - m_RLOFDetails.props2.radius1 = DEFAULT_INITIAL_DOUBLE_VALUE; - m_RLOFDetails.props2.radius2 = DEFAULT_INITIAL_DOUBLE_VALUE; + m_RLOFDetails.props2.radius1 = DEFAULT_INITIAL_DOUBLE_VALUE; + m_RLOFDetails.props2.radius2 = DEFAULT_INITIAL_DOUBLE_VALUE; - m_RLOFDetails.props2.starToRocheLobeRadiusRatio1 = DEFAULT_INITIAL_DOUBLE_VALUE; - m_RLOFDetails.props2.starToRocheLobeRadiusRatio2 = DEFAULT_INITIAL_DOUBLE_VALUE; + m_RLOFDetails.props2.starToRocheLobeRadiusRatio1 = DEFAULT_INITIAL_DOUBLE_VALUE; + m_RLOFDetails.props2.starToRocheLobeRadiusRatio2 = DEFAULT_INITIAL_DOUBLE_VALUE; - m_RLOFDetails.props2.semiMajorAxis = DEFAULT_INITIAL_DOUBLE_VALUE; - m_RLOFDetails.props2.eccentricity = DEFAULT_INITIAL_DOUBLE_VALUE; + m_RLOFDetails.props2.semiMajorAxis = DEFAULT_INITIAL_DOUBLE_VALUE; + m_RLOFDetails.props2.eccentricity = DEFAULT_INITIAL_DOUBLE_VALUE; - m_RLOFDetails.props2.eventCounter = DEFAULT_INITIAL_ULONGINT_VALUE; + m_RLOFDetails.props2.eventCounter = DEFAULT_INITIAL_ULONGINT_VALUE; - m_RLOFDetails.props2.time = DEFAULT_INITIAL_DOUBLE_VALUE; + m_RLOFDetails.props2.time = DEFAULT_INITIAL_DOUBLE_VALUE; - m_RLOFDetails.props2.isRLOF1 = false; - m_RLOFDetails.props2.isRLOF2 = false; + m_RLOFDetails.props2.isRLOF1 = false; + m_RLOFDetails.props2.isRLOF2 = false; - m_RLOFDetails.props2.isCE = false; + m_RLOFDetails.props2.isCE = false; // RLOF details - pre/post-MT props pointers - m_RLOFDetails.propsPostMT = &m_RLOFDetails.props1; - m_RLOFDetails.propsPreMT = &m_RLOFDetails.props2; + m_RLOFDetails.propsPostMT = &m_RLOFDetails.props1; + m_RLOFDetails.propsPreMT = &m_RLOFDetails.props2; // BeBinary details - properties 1 - m_BeBinaryDetails.props1.id = -1l; + m_BeBinaryDetails.props1.id = -1l; - m_BeBinaryDetails.props1.dt = DEFAULT_INITIAL_DOUBLE_VALUE; - m_BeBinaryDetails.props1.totalTime = DEFAULT_INITIAL_DOUBLE_VALUE; + m_BeBinaryDetails.props1.dt = DEFAULT_INITIAL_DOUBLE_VALUE; + m_BeBinaryDetails.props1.totalTime = DEFAULT_INITIAL_DOUBLE_VALUE; - m_BeBinaryDetails.props1.massNS = DEFAULT_INITIAL_DOUBLE_VALUE; + m_BeBinaryDetails.props1.massNS = DEFAULT_INITIAL_DOUBLE_VALUE; - m_BeBinaryDetails.props1.companionMass = DEFAULT_INITIAL_DOUBLE_VALUE; - m_BeBinaryDetails.props1.companionLuminosity = DEFAULT_INITIAL_DOUBLE_VALUE; - m_BeBinaryDetails.props1.companionTeff = DEFAULT_INITIAL_DOUBLE_VALUE; - m_BeBinaryDetails.props1.companionRadius = DEFAULT_INITIAL_DOUBLE_VALUE; + m_BeBinaryDetails.props1.companionMass = DEFAULT_INITIAL_DOUBLE_VALUE; + m_BeBinaryDetails.props1.companionLuminosity = DEFAULT_INITIAL_DOUBLE_VALUE; + m_BeBinaryDetails.props1.companionTeff = DEFAULT_INITIAL_DOUBLE_VALUE; + m_BeBinaryDetails.props1.companionRadius = DEFAULT_INITIAL_DOUBLE_VALUE; - m_BeBinaryDetails.props1.semiMajorAxis = DEFAULT_INITIAL_DOUBLE_VALUE; - m_BeBinaryDetails.props1.eccentricity = DEFAULT_INITIAL_DOUBLE_VALUE; + m_BeBinaryDetails.props1.semiMajorAxis = DEFAULT_INITIAL_DOUBLE_VALUE; + m_BeBinaryDetails.props1.eccentricity = DEFAULT_INITIAL_DOUBLE_VALUE; // BeBinary details - properties 2 - m_BeBinaryDetails.props2.id = -1l; + m_BeBinaryDetails.props2.id = -1l; - m_BeBinaryDetails.props2.dt = DEFAULT_INITIAL_DOUBLE_VALUE; - m_BeBinaryDetails.props2.totalTime = DEFAULT_INITIAL_DOUBLE_VALUE; + m_BeBinaryDetails.props2.dt = DEFAULT_INITIAL_DOUBLE_VALUE; + m_BeBinaryDetails.props2.totalTime = DEFAULT_INITIAL_DOUBLE_VALUE; - m_BeBinaryDetails.props2.massNS = DEFAULT_INITIAL_DOUBLE_VALUE; + m_BeBinaryDetails.props2.massNS = DEFAULT_INITIAL_DOUBLE_VALUE; - m_BeBinaryDetails.props2.companionMass = DEFAULT_INITIAL_DOUBLE_VALUE; - m_BeBinaryDetails.props2.companionLuminosity = DEFAULT_INITIAL_DOUBLE_VALUE; - m_BeBinaryDetails.props2.companionTeff = DEFAULT_INITIAL_DOUBLE_VALUE; - m_BeBinaryDetails.props2.companionRadius = DEFAULT_INITIAL_DOUBLE_VALUE; + m_BeBinaryDetails.props2.companionMass = DEFAULT_INITIAL_DOUBLE_VALUE; + m_BeBinaryDetails.props2.companionLuminosity = DEFAULT_INITIAL_DOUBLE_VALUE; + m_BeBinaryDetails.props2.companionTeff = DEFAULT_INITIAL_DOUBLE_VALUE; + m_BeBinaryDetails.props2.companionRadius = DEFAULT_INITIAL_DOUBLE_VALUE; - m_BeBinaryDetails.props2.semiMajorAxis = DEFAULT_INITIAL_DOUBLE_VALUE; - m_BeBinaryDetails.props2.eccentricity = DEFAULT_INITIAL_DOUBLE_VALUE; + m_BeBinaryDetails.props2.semiMajorAxis = DEFAULT_INITIAL_DOUBLE_VALUE; + m_BeBinaryDetails.props2.eccentricity = DEFAULT_INITIAL_DOUBLE_VALUE; // BeBinary details - current/prev props pointers - m_BeBinaryDetails.currentProps = &m_BeBinaryDetails.props1; - m_BeBinaryDetails.previousProps = &m_BeBinaryDetails.props2; + m_BeBinaryDetails.currentProps = &m_BeBinaryDetails.props1; + m_BeBinaryDetails.previousProps = &m_BeBinaryDetails.props2; // pointers - m_Donor = nullptr; - m_Accretor = nullptr; + m_Donor = nullptr; + m_Accretor = nullptr; - m_Supernova = nullptr; - m_Companion = nullptr; + m_Supernova = nullptr; + m_Companion = nullptr; } @@ -833,7 +819,7 @@ bool BaseBinaryStar::IsHMXRBinary() const { if (m_Star1->StellarType() < STELLAR_TYPE::NEUTRON_STAR && utils::Compare(StarToRocheLobeRadiusRatio1(), MIN_HMXRB_STAR_TO_ROCHE_LOBE_RADIUS_RATIO) > 0) return true; if (m_Star2->StellarType() < STELLAR_TYPE::NEUTRON_STAR && utils::Compare(StarToRocheLobeRadiusRatio2(), MIN_HMXRB_STAR_TO_ROCHE_LOBE_RADIUS_RATIO) > 0) return true; } - return false; + return false; } @@ -858,7 +844,7 @@ bool BaseBinaryStar::PrintRLOFParameters(const RLOF_RECORD_TYPE p_RecordType) { if (m_Star1->IsRLOF() || m_Star2->IsRLOF()) { // print if either star is in RLOF m_RLOFDetails.propsPostMT->eventCounter += 1; // every time we print a MT event happened, increment counter - ok = LOGGING->LogRLOFParameters(this, p_RecordType); // yes - write to log file + ok = LOGGING->LogRLOFParameters(this, p_RecordType); // yes - write to log file } if (OPTIONS->HMXRBinaries()) { @@ -905,7 +891,7 @@ bool BaseBinaryStar::PrintBeBinary(const BE_BINARY_RECORD_TYPE p_RecordType) { */ void BaseBinaryStar::StashRLOFProperties(const MASS_TRANSFER_TIMING p_Which) { - if (!OPTIONS->RLOFPrinting()) return; // nothing to do + if (!OPTIONS->RLOFPrinting()) return; // nothing to do // set whether to update pre-MT or post-MT parameters depending on input argument RLOFPropertiesT* rlofPropertiesToReset; @@ -924,7 +910,7 @@ void BaseBinaryStar::StashRLOFProperties(const MASS_TRANSFER_TIMING p_Which) { rlofPropertiesToReset->stellarType1 = m_Star1->StellarType(); rlofPropertiesToReset->stellarType2 = m_Star2->StellarType(); rlofPropertiesToReset->eccentricity = m_Eccentricity; - rlofPropertiesToReset->semiMajorAxis = m_SemiMajorAxis * AU_TO_RSOL; // semi-major axis - change units to Rsol + rlofPropertiesToReset->semiMajorAxis = m_SemiMajorAxis * AU_TO_RSOL; // semi-major axis - change units to Rsol rlofPropertiesToReset->time = m_Time; rlofPropertiesToReset->timePrev = m_TimePrev; rlofPropertiesToReset->isRLOF1 = m_Star1->IsRLOF(); @@ -948,8 +934,7 @@ void BaseBinaryStar::StashBeBinaryProperties() { if (!OPTIONS->BeBinaries() || !IsBeBinary()) return; // nothing to do; // switch previous<->current (preserves existing current as (new) previous) - BeBinaryPropertiesT* tmp; - tmp = m_BeBinaryDetails.previousProps; // save pointer to existing previous props + BeBinaryPropertiesT* tmp = m_BeBinaryDetails.previousProps; // save pointer to existing previous props m_BeBinaryDetails.previousProps = m_BeBinaryDetails.currentProps; // existing current props become new previous props (values will be preserved) m_BeBinaryDetails.currentProps = tmp; // new current props points at existing previous (values will be replaced) @@ -1137,36 +1122,33 @@ void BaseBinaryStar::ResolveCoalescence() { * * Note: the systemic speed is only valid for intact binaries, and component speeds are only valid for disrupted binaries. * - * ///////////////////////////////// - * // Logic - * // - * // If (Unbound before SN): - * // - * // Must be 2nd SN, only need to update starSN component velocity (rotated into previous reference frame). - * // - * // Else: (Intact before SN) - * // - * // Evolve binary according to vector algebra to determine centerofmass velocity, h', e', a', and whether bound or unbound. - * // - * // Update binary systemic velocity (even if disrupted, just for consistency) - rotate into previous reference frame if needed. - * // - * // If now unbound: - * // - * // Set m_Unbound to True - should be the only place in the code this is done. - * // - * // Continue vector algebra to find v1inf and v2inf. - * // Add these values to previous component velocities (rotated if need be) which will be the systemic velocity if this is the 2nd SN. - * // - * // For unbound binary, new Euler Angles should be randomized (see vector3d.cpp). - * // - * // If still intact: - * // - * // Binary systemic velocity has already been set, so just set the component velocities to the same vector. - * // (this is to make it easier to add just a component velocity later). - * // - * // For intact binary, Euler Angles must be calculated according to the vector algebra (see vector3d.h). - * // - * ///////////////////////////////////////////////////////////////////////////// + * Logic: + * + * if (Unbound before SN): + * + * Must be 2nd SN, only need to update starSN component velocity (rotated into previous reference frame). + * + * else: (Intact before SN) + * + * Evolve binary according to vector algebra to determine centerofmass velocity, h', e', a', and whether bound or unbound. + * Update binary systemic velocity (even if disrupted, just for consistency) - rotate into previous reference frame if needed. + * + * if now unbound: + * + * Set m_Unbound to True - should be the only place in the code this is done. + * + * Continue vector algebra to find v1inf and v2inf. + * Add these values to previous component velocities (rotated if need be) which will be the systemic velocity if this is the 2nd SN. + * + * For unbound binary, new Euler Angles should be randomized (see vector3d.cpp). + * + * if still intact: + * + * Binary systemic velocity has already been set, so just set the component velocities to the same vector. + * (this is to make it easier to add just a component velocity later). + * + * For intact binary, Euler Angles must be calculated according to the vector algebra (see vector3d.h). + * * * * bool ResolveSupernova() @@ -1177,57 +1159,51 @@ bool BaseBinaryStar::ResolveSupernova() { if (!m_Supernova->IsSNevent()) { SHOW_WARN(ERROR::RESOLVE_SUPERNOVA_IMPROPERLY_CALLED); - return false; // not a supernova event - bail out + return false; // not a supernova event - bail out } // Set relevant preSN parameters - m_EccentricityPreSN = m_Eccentricity; + m_EccentricityPreSN = m_Eccentricity; m_SemiMajorAxisPreSN = m_SemiMajorAxis; - double totalMassPreSN = m_Supernova->SN_TotalMassAtCOFormation() + m_Companion->Mass(); // Total Mass preSN - double reducedMassPreSN = m_Supernova->SN_TotalMassAtCOFormation() * m_Companion->Mass() / totalMassPreSN; // Reduced Mass preSN - m_Supernova->SetOrbitalEnergyPreSN(CalculateOrbitalEnergy(reducedMassPreSN, totalMassPreSN, m_SemiMajorAxisPreSN)); // Orbital energy preSN + double totalMassPreSN = m_Supernova->SN_TotalMassAtCOFormation() + m_Companion->Mass(); // total Mass preSN + double reducedMassPreSN = m_Supernova->SN_TotalMassAtCOFormation() * m_Companion->Mass() / totalMassPreSN; // reduced Mass preSN + m_Supernova->SetOrbitalEnergyPreSN(CalculateOrbitalEnergy(reducedMassPreSN, totalMassPreSN, m_SemiMajorAxisPreSN)); // orbital energy preSN // Define the natal kick vector (see above for precise definitions of the angles) - double theta = m_Supernova->SN_Theta(); // Angle out of the binary plane - double phi = m_Supernova->SN_Phi(); // Angle in the binary plane - Vector3d natalKickVector = m_Supernova->SN_KickMagnitude() *Vector3d(cos(theta)*cos(phi), - cos(theta)*sin(phi), + double theta = m_Supernova->SN_Theta(); // angle out of the binary plane + double phi = m_Supernova->SN_Phi(); // angle in the binary plane + Vector3d natalKickVector = m_Supernova->SN_KickMagnitude() *Vector3d(cos(theta) * cos(phi), + cos(theta) * sin(phi), sin(theta)); // Check if the system is already unbound - if (IsUnbound()) { // Is system already unbound? + if (IsUnbound()) { // is system already unbound? - m_Supernova->UpdateComponentVelocity( natalKickVector.RotateVector(m_ThetaE, m_PhiE, m_PsiE)); // yes - only need to update the velocity of the star undergoing SN + m_Supernova->UpdateComponentVelocity( natalKickVector.RotateVector(m_ThetaE, m_PhiE, m_PsiE)); // yes - only need to update the velocity of the star undergoing SN // The quantities below are meaningless in this context, so they are set to nan to avoid misuse m_OrbitalVelocityPreSN = -nan(""); - m_uK = nan(""); // -- - Dimensionless kick magnitude - + m_uK = nan(""); // -- - Dimensionless kick magnitude } - else { // no - evaluate orbital changes and calculate velocities - - ////////////////////////////////////////////////////////////////////////////////////////////////// - // - // Evolve SN out of binary - // - ////////////////////////////////////////////////////////////////////////////////////////////////// - + else { // no - evaluate orbital changes and calculate velocities + // Evolve SN out of binary // Functions defined in vector3d.h - #define cross(x,y) linalg::cross(x,y) - #define dot(x,y) linalg::dot(x,y) - #define angleBetween(x,y) linalg::angleBetween(x,y) - #define mag Magnitude() - #define hat UnitVector() + // Defined here for convenience - undefined later + #define cross(x,y) linalg::cross(x,y) + #define dot(x,y) linalg::dot(x,y) + #define angleBetween(x,y) linalg::angleBetween(x,y) + #define mag Magnitude() + #define hat UnitVector() // Pre-SN parameters - double semiMajorAxisPrev_km = m_SemiMajorAxis * AU_TO_KM; // km - Semi-Major axis - double eccentricityPrev = m_Eccentricity; // -- - Eccentricity, written with a prev to distinguish from later use - double sqrt1MinusEccPrevSquared = std::sqrt(1 - eccentricityPrev * eccentricityPrev); // useful function of eccentricity + double semiMajorAxisPrev_km = m_SemiMajorAxis * AU_TO_KM; // km - Semi-Major axis + double eccentricityPrev = m_Eccentricity; // -- - Eccentricity, written with a prev to distinguish from later use + double sqrt1MinusEccPrevSquared = std::sqrt(1.0 - eccentricityPrev * eccentricityPrev); // useful function of eccentricity - double m1Prev = m_Supernova->SN_TotalMassAtCOFormation(); // Mo - SN star pre-SN mass - double m2Prev = m_Companion->Mass(); // Mo - CP star pre-SN mass - double totalMassPrev = m1Prev + m2Prev; // Mo - Total binary pre-SN mass + double m1Prev = m_Supernova->SN_TotalMassAtCOFormation(); // Mo - SN star pre-SN mass + double m2Prev = m_Companion->Mass(); // Mo - CP star pre-SN mass + double totalMassPrev = m1Prev + m2Prev; // Mo - Total binary pre-SN mass // Functions of eccentric anomaly m_Supernova->CalculateSNAnomalies(eccentricityPrev); @@ -1235,93 +1211,82 @@ bool BaseBinaryStar::ResolveSupernova() { double sinEccAnomaly = sin(m_Supernova->SN_EccentricAnomaly()); // Derived quantities - double omega = std::sqrt(G_SN*totalMassPrev / (semiMajorAxisPrev_km * semiMajorAxisPrev_km*semiMajorAxisPrev_km)); // rad/s - Keplerian orbital frequency + double aPrev = semiMajorAxisPrev_km; + double aPrev_2 = aPrev * aPrev; + double aPrev_3 = aPrev_2 * aPrev; + + double omega = std::sqrt(G_SN * totalMassPrev / aPrev_3); // rad/s - Keplerian orbital frequency - Vector3d separationVectorPrev = Vector3d( semiMajorAxisPrev_km * (cosEccAnomaly - eccentricityPrev), - semiMajorAxisPrev_km * (sinEccAnomaly) * sqrt1MinusEccPrevSquared, - 0.0 ); // km - Relative position vector, from m1Prev to m2Prev - double separationPrev = separationVectorPrev.mag; // km - Instantaneous Separation + Vector3d separationVectorPrev = Vector3d(aPrev * (cosEccAnomaly - eccentricityPrev), + aPrev * (sinEccAnomaly) * sqrt1MinusEccPrevSquared, + 0.0); // km - Relative position vector, from m1Prev to m2Prev + double separationPrev = separationVectorPrev.mag; // km - Instantaneous Separation + double fact1 = aPrev_2 * omega / separationPrev; - Vector3d relativeVelocityVectorPrev = Vector3d(-((semiMajorAxisPrev_km * semiMajorAxisPrev_km) * omega / separationPrev) * sinEccAnomaly, - ((semiMajorAxisPrev_km * semiMajorAxisPrev_km) * omega / separationPrev) * cosEccAnomaly * sqrt1MinusEccPrevSquared, - 0.0 ); // km/s - Relative velocity vector, in the m1Prev rest frame + Vector3d relativeVelocityVectorPrev = Vector3d(-fact1 * sinEccAnomaly, + fact1 * cosEccAnomaly * sqrt1MinusEccPrevSquared, + 0.0); // km/s - Relative velocity vector, in the m1Prev rest frame - Vector3d orbitalAngularMomentumVectorPrev = cross(separationVectorPrev, relativeVelocityVectorPrev); // km^2 s^-1 - Specific orbital angular momentum vector + Vector3d orbitalAngularMomentumVectorPrev = cross(separationVectorPrev, relativeVelocityVectorPrev); // km^2 s^-1 - Specific orbital angular momentum vector - Vector3d eccentricityVectorPrev = cross(relativeVelocityVectorPrev, orbitalAngularMomentumVectorPrev) / (G_SN * totalMassPrev) - separationVectorPrev.hat; // -- - Laplace-Runge-Lenz vector (magnitude = eccentricity) + Vector3d eccentricityVectorPrev = cross(relativeVelocityVectorPrev, orbitalAngularMomentumVectorPrev) / (G_SN * totalMassPrev) - separationVectorPrev.hat; // -- - Laplace-Runge-Lenz vector (magnitude = eccentricity) - m_OrbitalVelocityPreSN = relativeVelocityVectorPrev.mag; // km/s - Set the Pre-SN orbital velocity and - m_uK = m_Supernova->SN_KickMagnitude() / m_OrbitalVelocityPreSN; // -- - Dimensionless kick magnitude + m_OrbitalVelocityPreSN = relativeVelocityVectorPrev.mag; // km/s - Set the Pre-SN orbital velocity and + m_uK = m_Supernova->SN_KickMagnitude() / m_OrbitalVelocityPreSN; // -- - Dimensionless kick magnitude - ///////////////////////////////////////////////////////////////////////////////////////// // Note: In the following, // orbitalAngularMomentumVectorPrev defines the Z-axis, // eccentricityVectorPrev defines the X-axis, and // (orbitalAngularMomentumVectorPrev x eccentricityVectorPrev) defines the Y-axis - ///////////////////////////////////////////////////////////////////////////////////////// - - ///////////////////////////////////////////////////////////////////////////////////////// // Apply supernova natal kick and mass loss // // Note: the code allows for mass loss and kick in the companion // (due to ablation), though we currently do not apply these. - ///////////////////////////////////////////////////////////////////////////////////////// - Vector3d companionRecoilVector = Vector3d(0.0, 0.0, 0.0); // km/s - The recoil of the companion due to ablation - double m1 = m_Supernova->Mass(); // Mo - supernova star postSN mass - double m2 = m_Companion->Mass(); // Mo - companion star postSN mass - double totalMass = m1 + m2; // Mo - Total binary postSN mass + Vector3d companionRecoilVector = Vector3d(0.0, 0.0, 0.0); // km/s - The recoil of the companion due to ablation - double dm1 = (m1Prev - m1); // Mo - Mass difference of supernova star - double dm2 = (m2Prev - m2); // Mo - Mass difference of companion star + double m1 = m_Supernova->Mass(); // Mo - supernova star postSN mass + double m2 = m_Companion->Mass(); // Mo - companion star postSN mass + double totalMass = m1 + m2; // Mo - Total binary postSN mass + double fact2 = totalMassPrev * totalMass; + double dm1 = (m1Prev - m1); // Mo - Mass difference of supernova star + double dm2 = (m2Prev - m2); // Mo - Mass difference of companion star - Vector3d centerOfMassVelocity = (-m2Prev * dm1 / (totalMassPrev*totalMass) + m1Prev * dm2 / (totalMassPrev * totalMass)) * relativeVelocityVectorPrev - + (m1 / totalMass) * natalKickVector - + (m2 / totalMass) * companionRecoilVector; // km/s - PostSN center of mass velocity vector + Vector3d centerOfMassVelocity = (-m2Prev * dm1 / fact2 + m1Prev * dm2 / fact2) * relativeVelocityVectorPrev + + (m1 / totalMass) * natalKickVector + (m2 / totalMass) * companionRecoilVector; // km/s - PostSN center of mass velocity vector - Vector3d relativeVelocityVector = relativeVelocityVectorPrev + (natalKickVector - companionRecoilVector); // km/s - PostSN relative velocity vector + Vector3d relativeVelocityVector = relativeVelocityVectorPrev + (natalKickVector - companionRecoilVector); // km/s - PostSN relative velocity vector - Vector3d orbitalAngularMomentumVector = cross(separationVectorPrev, relativeVelocityVector); // km^2 s^-1 - PostSN specific orbital angular momentum vector - double orbitalAngularMomentum = orbitalAngularMomentumVector.mag; // km^2 s^-1 - PostSN specific orbital angular momentum - m_OrbitalAngularMomentumVector = orbitalAngularMomentumVector/orbitalAngularMomentum; // set unit vector here to make printing out the inclination vector easier + Vector3d orbitalAngularMomentumVector = cross(separationVectorPrev, relativeVelocityVector); // km^2 s^-1 - PostSN specific orbital angular momentum vector + double orbitalAngularMomentum = orbitalAngularMomentumVector.mag; // km^2 s^-1 - PostSN specific orbital angular momentum + m_OrbitalAngularMomentumVector = orbitalAngularMomentumVector / orbitalAngularMomentum; // set unit vector here to make printing out the inclination vector easier - Vector3d eccentricityVector = cross(relativeVelocityVector, orbitalAngularMomentumVector) / (G_SN * totalMass) - - separationVectorPrev / separationPrev; // PostSN Laplace-Runge-Lenz vector - m_Eccentricity = eccentricityVector.mag; // PostSN eccentricity - double eccSquared = m_Eccentricity * m_Eccentricity; // useful function of eccentricity + Vector3d eccentricityVector = cross(relativeVelocityVector, orbitalAngularMomentumVector) / + (G_SN * totalMass) - separationVectorPrev / separationPrev; // PostSN Laplace-Runge-Lenz vector + m_Eccentricity = eccentricityVector.mag; // PostSN eccentricity + double eccSquared = m_Eccentricity * m_Eccentricity; // useful function of eccentricity - double semiMajorAxis_km = (orbitalAngularMomentum*orbitalAngularMomentum) / (G_SN * totalMass * (1 - eccSquared)); // km - PostSN semi-major axis - m_SemiMajorAxis = semiMajorAxis_km * KM_TO_AU; // AU - PostSN semi-major axis + double semiMajorAxis_km = (orbitalAngularMomentum * orbitalAngularMomentum) / (G_SN * totalMass * (1.0 - eccSquared)); // km - PostSN semi-major axis + m_SemiMajorAxis = semiMajorAxis_km * KM_TO_AU; // AU - PostSN semi-major axis - - ///////////////////////////////////////////////////////////////////////////////////////// // Note: similar to above, // orbitalAngularMomentumVector defines the Z'-axis, // eccentricityVector defines the X'-axis, and // (orbitalAngularMomentumVector x eccentricityVector) defines the Y'-axis - ///////////////////////////////////////////////////////////////////////////////////////// - UpdateSystemicVelocity(centerOfMassVelocity.RotateVector(m_ThetaE, m_PhiE, m_PsiE)); // Update the system velocity with the new center of mass velocity - + UpdateSystemicVelocity(centerOfMassVelocity.RotateVector(m_ThetaE, m_PhiE, m_PsiE)); // update the system velocity with the new center of mass velocity - ///////////////////////////////////////////////////////////////////////////////////////// // Split off and evaluate depending on whether the binary is now bound or unbound - if (utils::Compare(m_Eccentricity, 1.0) >= 0) { - - //////////////////////////////////////// - // - // Binary has become unbound - // - //////////////////////////////////////// - + if (utils::Compare(m_Eccentricity, 1.0) >= 0) { // unbound? + // yes, unbound m_Unbound = true; // Calculate the asymptotic Center of Mass velocity - double relativeVelocityAtInfinity = (G_SN*totalMass/orbitalAngularMomentum) * std::sqrt(eccSquared - 1); + double relativeVelocityAtInfinity = (G_SN*totalMass/orbitalAngularMomentum) * std::sqrt(eccSquared - 1.0); Vector3d relativeVelocityVectorAtInfinity = relativeVelocityAtInfinity - * (-1 * (eccentricityVector.hat / m_Eccentricity) - + std::sqrt(1 - 1.0 / eccSquared) * cross(orbitalAngularMomentumVector.hat, eccentricityVector.hat)); + * (-1.0 * (eccentricityVector.hat / m_Eccentricity) + + std::sqrt(1.0 - 1.0 / eccSquared) * cross(orbitalAngularMomentumVector.hat, eccentricityVector.hat)); // Calculate the asymptotic velocities of Star1 (SN) and Star2 (CP) Vector3d component1VelocityVectorAtInfinity = (m2 / totalMass) * relativeVelocityVectorAtInfinity + centerOfMassVelocity; @@ -1332,101 +1297,84 @@ bool BaseBinaryStar::ResolveSupernova() { m_Companion->UpdateComponentVelocity(component2VelocityVectorAtInfinity.RotateVector(m_ThetaE, m_PhiE, m_PsiE)); // Set Euler Angles - m_ThetaE = angleBetween(orbitalAngularMomentumVectorPrev, orbitalAngularMomentumVector); // Angle between the angular momentum unit vectors, always well defined + m_ThetaE = angleBetween(orbitalAngularMomentumVectorPrev, orbitalAngularMomentumVector); // angle between the angular momentum unit vectors, always well defined m_PhiE = _2_PI * RAND->Random(); m_PsiE = _2_PI * RAND->Random(); } - else { - - //////////////////////////////////////// - // - // Binary is still bound - // - //////////////////////////////////////// + else { // no - binary still bound // Set the component velocites to the system velocity. System velocity was already correctly set above. m_Supernova->UpdateComponentVelocity(centerOfMassVelocity.RotateVector(m_ThetaE, m_PhiE, m_PsiE)); m_Companion->UpdateComponentVelocity(centerOfMassVelocity.RotateVector(m_ThetaE, m_PhiE, m_PsiE)); - //////////////////////////////////////////////////////////////////////////////////// // Calculate Euler angles - see RotateVector() in vector.cpp for details - - m_ThetaE = angleBetween(orbitalAngularMomentumVector, orbitalAngularMomentumVectorPrev); // Angle between the angular momentum unit vectors, always well defined + m_ThetaE = angleBetween(orbitalAngularMomentumVector, orbitalAngularMomentumVectorPrev); // angle between the angular momentum unit vectors, always well defined // If the new orbital A.M. is parallel or anti-parallel to the previous orbital A.M., // then the cross product is not well-defined, and we need to account for degeneracy between eccentricity vectors. // Also, if either eccentricity is 0.0, then the eccentricity vector is not well defined. - if ((utils::Compare(m_ThetaE, 0.0) == 0) && // Is orbitalAngularMomentumVectorPrev parallel to orbitalAngularMomentumVector ... - ((utils::Compare(eccentricityPrev, 0.0) > 0) && // ... - (utils::Compare(m_Eccentricity, 0.0) > 0))) { // ...and both eccentricityVectorPrev and eccentricityVector are well defined? + if ((utils::Compare(m_ThetaE, 0.0) == 0) && // is orbitalAngularMomentumVectorPrev parallel to orbitalAngularMomentumVector ... + ((utils::Compare(eccentricityPrev, 0.0) > 0) && (utils::Compare(m_Eccentricity, 0.0) > 0))) { // ...and both eccentricityVectorPrev and eccentricityVector are well defined? - double psiPlusPhi = angleBetween(eccentricityVector, eccentricityVectorPrev); // yes - then psi + phi is constant - m_PhiE = _2_PI * RAND->Random(); - m_PsiE = psiPlusPhi - m_PhiE; + double psiPlusPhi = angleBetween(eccentricityVector, eccentricityVectorPrev); // yes - then psi + phi is constant + m_PhiE = _2_PI * RAND->Random(); + m_PsiE = psiPlusPhi - m_PhiE; } - else if ((utils::Compare(m_ThetaE, M_PI) == 0) && // Is orbitalAngularMomentumVectorPrev anti-parallel to orbitalAngularMomentumVector ... - ((utils::Compare(eccentricityPrev, 0.0) > 0) && // ... - (utils::Compare(m_Eccentricity, 0.0) > 0))) { // ...and both eccentricityVectorPrev and eccentricityVector are well defined? - - // yes - then psi - phi is constant - double psiMinusPhi = angleBetween(eccentricityVector, eccentricityVectorPrev); - m_PhiE = _2_PI * RAND->Random(); - m_PsiE = psiMinusPhi + m_PhiE; + else if ((utils::Compare(m_ThetaE, M_PI) == 0) && // is orbitalAngularMomentumVectorPrev anti-parallel to orbitalAngularMomentumVector ... + ((utils::Compare(eccentricityPrev, 0.0) > 0) && (utils::Compare(m_Eccentricity, 0.0) > 0))) { // ...and both eccentricityVectorPrev and eccentricityVector are well defined? + + double psiMinusPhi = angleBetween(eccentricityVector, eccentricityVectorPrev); // yes - then psi - phi is constant + m_PhiE = _2_PI * RAND->Random(); + m_PsiE = psiMinusPhi + m_PhiE; } - else { // Neither - the cross product of the orbit normals is well-defined + else { // neither - the cross product of the orbit normals is well-defined - Vector3d orbitalPivotAxis = cross(orbitalAngularMomentumVectorPrev, orbitalAngularMomentumVector); // Cross product of the orbit normals + Vector3d orbitalPivotAxis = cross(orbitalAngularMomentumVectorPrev, orbitalAngularMomentumVector); // cross product of the orbit normals - if ( utils::Compare(eccentricityPrev, 0.0) == 0 ) { // Is eccentricityVectorPrev well-defined? - m_PhiE = _2_PI * RAND->Random(); // no - set phi random + if (utils::Compare(eccentricityPrev, 0.0) == 0 ) { // is eccentricityVectorPrev well-defined? + m_PhiE = _2_PI * RAND->Random(); // no - set phi random } - else { // yes - phi is +/- angle between eccentricityVectorPrev and orbitalPivotAxis + else { // yes - phi is +/- angle between eccentricityVectorPrev and orbitalPivotAxis - m_PhiE = utils::Compare( dot(eccentricityVectorPrev, orbitalAngularMomentumVector), 0.0) >= 0 ? // Are eccentricityVectorPrev and orbitalAngularMomentumVector in the same hemisphere? - angleBetween(eccentricityVectorPrev, orbitalPivotAxis): // yes - phi in [0,pi) - -angleBetween(eccentricityVectorPrev, orbitalPivotAxis); // no - phi in [-pi,0) + m_PhiE = utils::Compare( dot(eccentricityVectorPrev, orbitalAngularMomentumVector), 0.0) >= 0 // are eccentricityVectorPrev and orbitalAngularMomentumVector in the same hemisphere? + ? angleBetween(eccentricityVectorPrev, orbitalPivotAxis) // yes - phi in [0,pi) + : -angleBetween(eccentricityVectorPrev, orbitalPivotAxis); // no - phi in [-pi,0) } - if ( utils::Compare(m_Eccentricity, 0.0) == 0 ) { // Is eccentricityVector well-defined? - m_PsiE = _2_PI * RAND->Random(); // no - set psi random + if ( utils::Compare(m_Eccentricity, 0.0) == 0 ) { // is eccentricityVector well-defined? + m_PsiE = _2_PI * RAND->Random(); // no - set psi random } - else { // yes - psi is +/- angle between eccentricityVector and orbitalPivotAxis - - m_PsiE = utils::Compare( dot(eccentricityVector, orbitalAngularMomentumVectorPrev), 0.0) >= 0 ? // Are eccentricityVector and orbitalAngularMomentumVectorPrev in the same hemisphere? - angleBetween(eccentricityVector, orbitalPivotAxis): // yes - psi in [0,pi) - -angleBetween(eccentricityVector, orbitalPivotAxis); // no - psi in [-pi,0) + else { // yes - psi is +/- angle between eccentricityVector and orbitalPivotAxis + m_PsiE = utils::Compare( dot(eccentricityVector, orbitalAngularMomentumVectorPrev), 0.0) >= 0 // are eccentricityVector and orbitalAngularMomentumVectorPrev in the same hemisphere? + ? angleBetween(eccentricityVector, orbitalPivotAxis) // yes - psi in [0,pi) + : -angleBetween(eccentricityVector, orbitalPivotAxis); // no - psi in [-pi,0) } } // Note: There is some evidence for evolution of periapsis in mass transferring binaries (see e.g Dosopoulou & Kalogera 2016, 2018). - // This should be investigated in more depth, but until then, we assume that the periapsis *may* evolve, - // and accordingly randomize the angle of periapsis around the new orbital angular momentum, (i.e, Psi) - // - RTW 15/05/20 + // This should be investigated in more depth, but until then, we assume that the periapsis *may* evolve, and accordingly randomize + // the angle of periapsis around the new orbital angular momentum, (i.e, Psi) - RTW 15/05/20 m_PsiE = _2_PI * RAND->Random(); } - // Undefine the pre-processor commands - #undef cross - #undef dot - #undef angleBetween - #undef mag #undef hat + #undef mag + #undef angleBetween + #undef dot + #undef cross } - ////////////////////////// - // Do for all systems - // Set remaining post-SN values - double totalMass = m_Supernova->Mass() + m_Companion->Mass(); // Total Mass - double reducedMass = m_Supernova->Mass() * m_Companion->Mass() / totalMass; // Reduced Mass - m_Supernova->SetOrbitalEnergyPostSN(CalculateOrbitalEnergy(reducedMass, totalMass, m_SemiMajorAxis)); // Orbital energy + double totalMass = m_Supernova->Mass() + m_Companion->Mass(); // total Mass + double reducedMass = m_Supernova->Mass() * m_Companion->Mass() / totalMass; // reduced Mass + m_Supernova->SetOrbitalEnergyPostSN(CalculateOrbitalEnergy(reducedMass, totalMass, m_SemiMajorAxis)); // orbital energy - m_IPrime = m_ThetaE; // Inclination angle between preSN and postSN orbital planes + m_IPrime = m_ThetaE; // inclination angle between preSN and postSN orbital planes m_CosIPrime = cos(m_IPrime); - (void)PrintSupernovaDetails(); // Log record to supernovae logfile + (void)PrintSupernovaDetails(); // log record to supernovae logfile m_Supernova->ClearCurrentSNEvent(); return true; @@ -1501,16 +1449,16 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { double alphaCE = OPTIONS->CommonEnvelopeAlpha(); // CE efficiency parameter - double eccentricity = Eccentricity(); // current eccentricity (before CEE) - double semiMajorAxisRsol= SemiMajorAxisRsol(); // current semi-major axis in default units, Rsol (before CEE) - double periastronRsol = PeriastronRsol(); // periastron, Rsol (before CEE) - double rRLd1Rsol = periastronRsol * CalculateRocheLobeRadius_Static(m_Star1->Mass(), m_Star2->Mass()); // Roche-lobe radius at periastron in Rsol at the moment where CEE begins, seen by star1 - double rRLd2Rsol = periastronRsol * CalculateRocheLobeRadius_Static(m_Star2->Mass(), m_Star1->Mass()); // Roche-lobe radius at periastron in Rsol at the moment where CEE begins, seen by star2 + double eccentricity = Eccentricity(); // current eccentricity (before CEE) + double semiMajorAxisRsol = SemiMajorAxisRsol(); // current semi-major axis in default units, Rsol (before CEE) + double periastronRsol = PeriastronRsol(); // periastron, Rsol (before CEE) + double rRLd1Rsol = periastronRsol * CalculateRocheLobeRadius_Static(m_Star1->Mass(), m_Star2->Mass()); // Roche-lobe radius at periastron in Rsol at the moment where CEE begins, seen by star1 + double rRLd2Rsol = periastronRsol * CalculateRocheLobeRadius_Static(m_Star2->Mass(), m_Star1->Mass()); // Roche-lobe radius at periastron in Rsol at the moment where CEE begins, seen by star2 bool isDonorMS = false; // check for main sequence donor if (OPTIONS->AllowMainSequenceStarToSurviveCommonEnvelope()) { // allow main sequence stars to survive CEE? if (m_Star1->IsOneOf(ALL_MAIN_SEQUENCE)) { // yes - star1 MS_LTE_07, MS_GT_07 or NAKED_HELIUM_STAR_MS? - isDonorMS = isDonorMS || m_Star1->IsRLOF(); // yes - donor MS? + isDonorMS = isDonorMS || m_Star1->IsRLOF(); // yes - donor MS? m_Mass1Final = m_Star1->Mass(); // set mass m_MassEnv1 = 0.0; // no envelope } @@ -1520,7 +1468,7 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { } if (m_Star2->IsOneOf(ALL_MAIN_SEQUENCE)) { // star2 MS_LTE_07, MS_GT_07 or NAKED_HELIUM_STAR_MS? - isDonorMS = isDonorMS || m_Star2->IsRLOF(); // yes - donor MS? + isDonorMS = isDonorMS || m_Star2->IsRLOF(); // yes - donor MS? m_Mass2Final = m_Star2->Mass(); // yes - set mass m_MassEnv2 = 0.0; // no envelope } @@ -1563,16 +1511,16 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { // due to the CEE as described in Belczynsky et al. 2002, eq. (12) if( OPTIONS->CommonEnvelopeFormalism() == CE_FORMALISM::ENERGY ) { - double k1 = m_Star1->IsOneOf(COMPACT_OBJECTS) ? 0.0 : (2.0 / (lambda1 * alphaCE)) * m_Star1->Mass() * m_MassEnv1 / m_Star1->Radius(); - double k2 = m_Star2->IsOneOf(COMPACT_OBJECTS) ? 0.0 : (2.0 / (lambda2 * alphaCE)) * m_Star2->Mass() * m_MassEnv2 / m_Star2->Radius(); - double k3 = m_Star1->Mass() * m_Star2->Mass() / periastronRsol; //assumes immediate circularisation at periastron at start of CE - double k4 = (m_Mass1Final * m_Mass2Final); - double aFinalRsol = k4 / (k1 + k2 + k3); - m_SemiMajorAxis = aFinalRsol*RSOL_TO_AU; + double k1 = m_Star1->IsOneOf(COMPACT_OBJECTS) ? 0.0 : (2.0 / (lambda1 * alphaCE)) * m_Star1->Mass() * m_MassEnv1 / m_Star1->Radius(); + double k2 = m_Star2->IsOneOf(COMPACT_OBJECTS) ? 0.0 : (2.0 / (lambda2 * alphaCE)) * m_Star2->Mass() * m_MassEnv2 / m_Star2->Radius(); + double k3 = m_Star1->Mass() * m_Star2->Mass() / periastronRsol; //assumes immediate circularisation at periastron at start of CE + double k4 = (m_Mass1Final * m_Mass2Final); + double aFinalRsol = k4 / (k1 + k2 + k3); + m_SemiMajorAxis = aFinalRsol * RSOL_TO_AU; } - + // Two-stage common envelope, Hirai & Mandel (2022) - else if( OPTIONS->CommonEnvelopeFormalism() == CE_FORMALISM::TWO_STAGE ) { + else if ( OPTIONS->CommonEnvelopeFormalism() == CE_FORMALISM::TWO_STAGE ) { double convectiveEnvelopeMass1 = m_Star1->CalculateConvectiveEnvelopeMass(); double radiativeIntershellMass1 = m_MassEnv1 - convectiveEnvelopeMass1; double endOfFirstStageMass1 = m_Mass1Final + radiativeIntershellMass1; @@ -1581,35 +1529,34 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { double endOfFirstStageMass2 = m_Mass2Final + radiativeIntershellMass2; // Stage 1: convective envelope removal on a dynamical timescale; assumes lambda = 1.5, motivated by bottom panel of Figure 3 of Hirai & Mandel 2022, including internal energy - double k1 = m_Star1->IsOneOf(COMPACT_OBJECTS) ? 0.0 : (2.0 / (1.5 * alphaCE)) * m_Star1->Mass() * convectiveEnvelopeMass1 / m_Star1->Radius(); - double k2 = m_Star2->IsOneOf(COMPACT_OBJECTS) ? 0.0 : (2.0 / (1.5 * alphaCE)) * m_Star2->Mass() * convectiveEnvelopeMass2 / m_Star2->Radius(); - double k3 = m_Star1->Mass() * m_Star2->Mass() / periastronRsol; //assumes immediate circularisation at periastron at start of CE - double k4 = (endOfFirstStageMass1 * endOfFirstStageMass2); - double aFinalRsol = k4 / (k1 + k2 + k3); - m_SemiMajorAxis = aFinalRsol*RSOL_TO_AU; - + double k1 = m_Star1->IsOneOf(COMPACT_OBJECTS) ? 0.0 : (2.0 / (1.5 * alphaCE)) * m_Star1->Mass() * convectiveEnvelopeMass1 / m_Star1->Radius(); + double k2 = m_Star2->IsOneOf(COMPACT_OBJECTS) ? 0.0 : (2.0 / (1.5 * alphaCE)) * m_Star2->Mass() * convectiveEnvelopeMass2 / m_Star2->Radius(); + double k3 = m_Star1->Mass() * m_Star2->Mass() / periastronRsol; //assumes immediate circularisation at periastron at start of CE + double k4 = (endOfFirstStageMass1 * endOfFirstStageMass2); + double aFinalRsol = k4 / (k1 + k2 + k3); + m_SemiMajorAxis = aFinalRsol*RSOL_TO_AU; + // Stage 2: radiative envelope removal on a thermal timescale; assumed to be fully non-conservative if( utils::Compare(radiativeIntershellMass1, 0.0) > 0 ) { m_SemiMajorAxis = CalculateMassTransferOrbit(endOfFirstStageMass1, -radiativeIntershellMass1, *m_Star2, 0.0); } + if( utils::Compare(radiativeIntershellMass2, 0.0) > 0 ) { m_SemiMajorAxis = CalculateMassTransferOrbit(endOfFirstStageMass2, -radiativeIntershellMass2, *m_Star1, 0.0); } - } - - else { // Invalid CE formalism + } + else { // invalid CE formalism SHOW_WARN_STATIC(ERROR::UNKNOWN_CE_FORMALISM, // show warning "Orbital properties unchanged by CE", OBJECT_TYPE::BASE_BINARY_STAR, STELLAR_TYPE::BINARY_STAR); } - - double rRLdfin1 = m_SemiMajorAxis * CalculateRocheLobeRadius_Static(m_Mass1Final, m_Mass2Final); // Roche-lobe radius in AU after CEE, seen by star1 - double rRLdfin2 = m_SemiMajorAxis * CalculateRocheLobeRadius_Static(m_Mass2Final, m_Mass1Final); // Roche-lobe radius in AU after CEE, seen by star2 - double rRLdfin1Rsol = rRLdfin1 * AU_TO_RSOL; // Roche-lobe radius in Rsol after CEE, seen by star1 - double rRLdfin2Rsol = rRLdfin2 * AU_TO_RSOL; // Roche-lobe radius in Rsol after CEE, seen by star2 - m_Eccentricity = 0.0; // We assume that a common envelope event (CEE) circularises the binary + double rRLdfin1 = m_SemiMajorAxis * CalculateRocheLobeRadius_Static(m_Mass1Final, m_Mass2Final); // Roche-lobe radius in AU after CEE, seen by star1 + double rRLdfin2 = m_SemiMajorAxis * CalculateRocheLobeRadius_Static(m_Mass2Final, m_Mass1Final); // Roche-lobe radius in AU after CEE, seen by star2 + double rRLdfin1Rsol = rRLdfin1 * AU_TO_RSOL; // Roche-lobe radius in Rsol after CEE, seen by star1 + double rRLdfin2Rsol = rRLdfin2 * AU_TO_RSOL; // Roche-lobe radius in Rsol after CEE, seen by star2 + m_Eccentricity = 0.0; // we assume that a common envelope event (CEE) circularises the binary m_Star1->ResolveCommonEnvelopeAccretion(m_Mass1Final); // update star1's mass after CE accretion m_Star2->ResolveCommonEnvelopeAccretion(m_Mass2Final); // update star2's mass after CE accretion @@ -1664,9 +1611,9 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { m_Star2->SetPostCEEValues(); // squirrel away post CEE stellar values for star 2 SetPostCEEValues(m_SemiMajorAxis * AU_TO_RSOL, m_Eccentricity, rRLdfin1Rsol, rRLdfin2Rsol); // squirrel away post CEE binary values (checks for post-CE RLOF, so should be done at end) - if (m_RLOFDetails.immediateRLOFPostCEE == true && !OPTIONS->AllowImmediateRLOFpostCEToSurviveCommonEnvelope()) { // Is there immediate post-CE RLOF which is not allowed? + if (m_RLOFDetails.immediateRLOFPostCEE == true && !OPTIONS->AllowImmediateRLOFpostCEToSurviveCommonEnvelope()) { // is there immediate post-CE RLOF which is not allowed? m_MassTransferTrackerHistory = MT_TRACKING::MERGER; - m_Flags.stellarMerger = true; + m_Flags.stellarMerger = true; } (void)PrintCommonEnvelope(); // print (log) common envelope details @@ -1686,8 +1633,8 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { * @return Radius of Roche Lobe in units of the semi-major axis a */ double BaseBinaryStar::CalculateRocheLobeRadius_Static(const double p_MassPrimary, const double p_MassSecondary) { - double q = p_MassPrimary / p_MassSecondary; - double qCubeRoot = PPOW(q, 1.0 / 3.0); // cube roots are expensive, only compute once + double q = p_MassPrimary / p_MassSecondary; + double qCubeRoot = std::cbrt(q); // cube roots are expensive, only compute once return 0.49 / (0.6 + log(1.0 + qCubeRoot) / qCubeRoot / qCubeRoot); } @@ -1719,18 +1666,18 @@ double BaseBinaryStar::CalculateGammaAngularMomentumLoss(const double p_DonorMas case MT_ANGULAR_MOMENTUM_LOSS_PRESCRIPTION::JEANS : gamma = p_AccretorMass / p_DonorMass; break; // vicinity of the donor case MT_ANGULAR_MOMENTUM_LOSS_PRESCRIPTION::ISOTROPIC_RE_EMISSION: gamma = p_DonorMass / p_AccretorMass; break; // vicinity of the accretor case MT_ANGULAR_MOMENTUM_LOSS_PRESCRIPTION::CIRCUMBINARY_RING : gamma = (M_SQRT2 * (p_DonorMass + p_AccretorMass) * (p_DonorMass + p_AccretorMass)) / (p_DonorMass * p_AccretorMass); break; // Based on the assumption that a_ring ~= 2*a*, Vinciguerra+, 2020 - case MT_ANGULAR_MOMENTUM_LOSS_PRESCRIPTION::MACLEOD_LINEAR : { // Linear interpolation on separation between accretor and L2 point + case MT_ANGULAR_MOMENTUM_LOSS_PRESCRIPTION::MACLEOD_LINEAR : { // linear interpolation on separation between accretor and L2 point double q = p_AccretorMass / p_DonorMass; // interpolate in separation between a_acc and a_L2, both normalized to units of separation a - double aL2 = std::sqrt(M_SQRT2); // roughly, coincides with CIRCUMBINARY_RING def above - double aAcc = 1/(1+q); + double aL2 = std::sqrt(M_SQRT2); // roughly, coincides with CIRCUMBINARY_RING def above + double aAcc = 1.0 / (1.0 + q); double aGamma = aAcc + (aL2 - aAcc)*OPTIONS->MassTransferJlossMacLeodLinearFraction(); - gamma = aGamma*aGamma*(1+q)*(1+q)/q; + gamma = aGamma * aGamma * (1.0 + q) * (1.0 + q) / q; break; } case MT_ANGULAR_MOMENTUM_LOSS_PRESCRIPTION::ARBITRARY : gamma = OPTIONS->MassTransferJloss(); break; default: // unknown mass transfer angular momentum loss prescription - shouldn't happen - gamma = 1.0; // default value + gamma = 1.0; // default value m_Error = ERROR::UNKNOWN_MT_ANGULAR_MOMENTUM_LOSS_PRESCRIPTION; // set error value SHOW_WARN(m_Error); // warn that an error occurred } @@ -1761,30 +1708,28 @@ double BaseBinaryStar::CalculateMassTransferOrbit(const double p BinaryConstituentStar& p_Accretor, const double p_FractionAccreted) { - double semiMajorAxis = m_SemiMajorAxis; // new semi-major axis value - default is no change - double massA = p_Accretor.Mass(); // accretor mass - double massD = p_DonorMass; // donor mass - double massAtimesMassD = massA * massD; // accretor mass * donor mass - double massAplusMassD = massA + massD; // accretor mass + donor mass - double jOrb = (massAtimesMassD / massAplusMassD) * std::sqrt(semiMajorAxis * G1 * massAplusMassD); // orbital angular momentum - double jLoss; // specific angular momentum carried away by non-conservative mass transfer + double semiMajorAxis = m_SemiMajorAxis; // new semi-major axis value - default is no change + double massA = p_Accretor.Mass(); // accretor mass + double massD = p_DonorMass; // donor mass + double massAtimesMassD = massA * massD; // accretor mass * donor mass + double massAplusMassD = massA + massD; // accretor mass + donor mass + double jOrb = (massAtimesMassD / massAplusMassD) * std::sqrt(semiMajorAxis * G1 * massAplusMassD); // orbital angular momentum + double jLoss; // specific angular momentum carried away by non-conservative mass transfer - if (utils::Compare(p_DeltaMassDonor, 0.0) >= 0) { // no mass loss from donor, nothing to do here - return semiMajorAxis; - } - int numberIterations = fmax( floor (fabs(p_DeltaMassDonor/(MAXIMUM_MASS_TRANSFER_FRACTION_PER_STEP*massD))), 1); // number of iterations - - double dM = p_DeltaMassDonor / numberIterations; // mass change per time step + if (utils::Compare(p_DeltaMassDonor, 0.0) < 0) { // mass loss from donor? + // yes + int numberIterations = fmax( floor (fabs(p_DeltaMassDonor/(MAXIMUM_MASS_TRANSFER_FRACTION_PER_STEP * massD))), 1.0); // number of iterations + double dM = p_DeltaMassDonor / numberIterations; // mass change per time step - for(int i = 0; i < numberIterations ; i++) { - - jLoss = CalculateGammaAngularMomentumLoss(massD, massA); - jOrb = jOrb + ((jLoss * jOrb * (1.0 - p_FractionAccreted) / massAplusMassD) * dM); - semiMajorAxis = semiMajorAxis + (((-2.0 * dM / massD) * (1.0 - (p_FractionAccreted * (massD / massA)) - ((1.0 - p_FractionAccreted) * (jLoss + 0.5) * (massD / massAplusMassD)))) * semiMajorAxis); + for(int i = 0; i < numberIterations ; i++) { + jLoss = CalculateGammaAngularMomentumLoss(massD, massA); + jOrb = jOrb + ((jLoss * jOrb * (1.0 - p_FractionAccreted) / massAplusMassD) * dM); + semiMajorAxis = semiMajorAxis + (((-2.0 * dM / massD) * (1.0 - (p_FractionAccreted * (massD / massA)) - ((1.0 - p_FractionAccreted) * (jLoss + 0.5) * (massD / massAplusMassD)))) * semiMajorAxis); - massD = massD + dM; - massA = massA - (dM * p_FractionAccreted); - massAplusMassD = massA + massD; + massD = massD + dM; + massA = massA - (dM * p_FractionAccreted); + massAplusMassD = massA + massD; + } } return semiMajorAxis; @@ -1811,13 +1756,11 @@ double BaseBinaryStar::CalculateZetaRocheLobe(const double p_jLoss) const { double accretorMass = m_Accretor->Mass(); // accretor mass double beta = m_FractionAccreted; // fraction of mass accreted by accretor double gamma = p_jLoss; - - double q = donorMass / accretorMass; - - double q_1_3 = PPOW(q, 1.0 / 3.0); + double q = donorMass / accretorMass; + double cbrt_q = std::cbrt(q); double k1 = -2.0 * (1.0 - (beta * q) - (1.0 - beta) * (gamma + 0.5) * (q / (1.0 + q))); - double k2 = (2.0 / 3.0) - q_1_3 * (1.2 * q_1_3 + 1.0 / (1.0 + q_1_3)) / (3.0 * (0.6 * q_1_3 * q_1_3 + log(1.0 + q_1_3))); + double k2 = (2.0 / 3.0) - cbrt_q * (1.2 * cbrt_q + 1.0 / (1.0 + cbrt_q)) / (3.0 * (0.6 * cbrt_q * cbrt_q + log(1.0 + cbrt_q))); double k3 = 1.0 + (beta * q); return k1 + (k2 * k3); @@ -1903,8 +1846,8 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) { // Calculate accretion fraction if stable // This passes the accretor's Roche lobe radius to m_Accretor->CalculateThermalMassAcceptanceRate() // just in case MT_THERMALLY_LIMITED_VARIATION::RADIUS_TO_ROCHELOBE is used; otherwise, the radius input is ignored - double accretorRLradius = CalculateRocheLobeRadius_Static(m_Accretor->Mass(), m_Donor->Mass()) * AU_TO_RSOL * m_SemiMajorAxis * (1.0 - m_Eccentricity); - bool donorIsHeRich = m_Donor->IsOneOf(He_RICH_TYPES); + double accretorRLradius = CalculateRocheLobeRadius_Static(m_Accretor->Mass(), m_Donor->Mass()) * AU_TO_RSOL * m_SemiMajorAxis * (1.0 - m_Eccentricity); + bool donorIsHeRich = m_Donor->IsOneOf(He_RICH_TYPES); std::tie(std::ignore, m_FractionAccreted) = m_Accretor->CalculateMassAcceptanceRate(m_Donor->CalculateThermalMassLossRate(), m_Accretor->CalculateThermalMassAcceptanceRate(accretorRLradius), donorIsHeRich); @@ -1934,7 +1877,7 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) { isUnstable = true; if (!m_Donor->IsOneOf(GIANTS)) m_Flags.stellarMerger = true; } - else if (OPTIONS->QCritPrescription() != QCRIT_PRESCRIPTION::NONE) { // Determine stability based on critical mass ratios + else if (OPTIONS->QCritPrescription() != QCRIT_PRESCRIPTION::NONE) { // Determine stability based on critical mass ratios // NOTE: Critical mass ratio is defined as mAccretor/mDonor double qCrit = m_Donor->CalculateCriticalMassRatio(m_Accretor->IsDegenerate()); @@ -1961,11 +1904,11 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) { : MT_TRACKING::STABLE_2_TO_1_SURV; double massDiffDonor; - double envMassDonor = m_Donor->Mass() - m_Donor->CoreMass(); + double envMassDonor = m_Donor->Mass() - m_Donor->CoreMass(); bool isEnvelopeRemoved = false; if (utils::Compare(m_Donor->CoreMass(), 0) > 0 && utils::Compare(envMassDonor, 0) > 0) { // donor has a core and an envelope - massDiffDonor = -envMassDonor; // set donor mass loss to (negative of) the envelope mass + massDiffDonor = -envMassDonor; // set donor mass loss to (negative of) the envelope mass isEnvelopeRemoved = true; } else{ // donor has no envelope @@ -1974,10 +1917,10 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) { } double massGainAccretor = -massDiffDonor * m_FractionAccreted; // set accretor mass gain to mass loss * conservativeness - m_Donor->SetMassTransferDiffAndResolveWDShellChange(massDiffDonor); // set new mass of donor - m_Accretor->SetMassTransferDiffAndResolveWDShellChange(massGainAccretor); // set new mass of accretor + m_Donor->SetMassTransferDiffAndResolveWDShellChange(massDiffDonor); // set new mass of donor + m_Accretor->SetMassTransferDiffAndResolveWDShellChange(massGainAccretor); // set new mass of accretor - aFinal = CalculateMassTransferOrbit(m_Donor->Mass(), massDiffDonor, *m_Accretor, m_FractionAccreted); // calculate new orbit + aFinal = CalculateMassTransferOrbit(m_Donor->Mass(), massDiffDonor, *m_Accretor, m_FractionAccreted); // calculate new orbit m_aMassTransferDiff = aFinal - aInitial; // set change in orbit (semi-major axis) STELLAR_TYPE stellarTypeDonor = m_Donor->StellarType(); // donor stellar type before resolving envelope loss @@ -2070,7 +2013,7 @@ void BaseBinaryStar::InitialiseMassTransfer() { // If you don't do this, you end up modifying pre-MT pre-circularisation orbit // JR: todo: check that this is proper functionality, or just a kludge - if kludge, resolve it m_SemiMajorAxisPrev = m_SemiMajorAxis; - m_EccentricityPrev = m_Eccentricity; + m_EccentricityPrev = m_Eccentricity; } } } @@ -2179,6 +2122,7 @@ double BaseBinaryStar::CalculateAngularMomentum(const double p_SemiMajorAxis, const double p_Star2_SpinAngularVelocity, const double p_Star1_GyrationRadius, const double p_Star2_GyrationRadius) const { + double m1 = p_Star1Mass; double m2 = p_Star2Mass; @@ -2213,17 +2157,18 @@ void BaseBinaryStar::CalculateEnergyAndAngularMomentum() { if (m_Star1->IsOneOf({ STELLAR_TYPE::MASSLESS_REMNANT }) || m_Star2->IsOneOf({ STELLAR_TYPE::MASSLESS_REMNANT })) return; // Calculate orbital energy and angular momentum - m_OrbitalEnergyPrev = m_OrbitalEnergy; - m_OrbitalAngularMomentumPrev = m_OrbitalAngularMomentum; + m_OrbitalEnergyPrev = m_OrbitalEnergy; + m_OrbitalAngularMomentumPrev = m_OrbitalAngularMomentum; + m_TotalAngularMomentumPrev = m_TotalAngularMomentum; - double totalMass = m_Star1->Mass() + m_Star2->Mass(); - double reducedMass = (m_Star1->Mass() * m_Star2->Mass()) / totalMass; - m_OrbitalEnergy = CalculateOrbitalEnergy(reducedMass, totalMass, m_SemiMajorAxis); - m_OrbitalAngularMomentum = CalculateOrbitalAngularMomentum(reducedMass, totalMass, m_SemiMajorAxis); + double totalMass = m_Star1->Mass() + m_Star2->Mass(); + double reducedMass = (m_Star1->Mass() * m_Star2->Mass()) / totalMass; + m_OrbitalEnergy = CalculateOrbitalEnergy(reducedMass, totalMass, m_SemiMajorAxis); + m_OrbitalAngularMomentum = CalculateOrbitalAngularMomentum(reducedMass, totalMass, m_SemiMajorAxis); // Calculate total energy and angular momentum using regular conservation of energy, especially useful for checking tides and rotational effects - m_TotalEnergy = CalculateTotalEnergy(); - m_TotalAngularMomentum = CalculateAngularMomentum(); + m_TotalEnergy = CalculateTotalEnergy(); + m_TotalAngularMomentum = CalculateAngularMomentum(); } @@ -2245,14 +2190,14 @@ void BaseBinaryStar::ResolveMassChanges() { STELLAR_TYPE stellarType2 = m_Star2->StellarTypePrev(); // star 2 stellar type before updating attributes // update mass of star1 according to mass loss and mass transfer, then update age accordingly - (void)m_Star1->UpdateAttributes(m_Star1->MassPrev() - m_Star1->Mass() + m_Star1->MassLossDiff() + m_Star1->MassTransferDiff(), 0.0); // update mass for star1 + (void)m_Star1->UpdateAttributes(m_Star1->MassPrev() - m_Star1->Mass() + m_Star1->MassLossDiff() + m_Star1->MassTransferDiff(), 0.0); // update mass for star1 m_Star1->UpdateInitialMass(); // update effective initial mass of star1 (MS, HG & HeMS) m_Star1->UpdateAgeAfterMassLoss(); // update age of star1 m_Star1->ApplyMassTransferRejuvenationFactor(); // apply age rejuvenation factor for star1 m_Star1->UpdateAttributes(0.0, 0.0, true); // rinse and repeat for star2 - (void)m_Star2->UpdateAttributes(m_Star2->MassPrev() - m_Star2->Mass() + m_Star2->MassLossDiff() + m_Star2->MassTransferDiff(), 0.0); // update mass for star2 + (void)m_Star2->UpdateAttributes(m_Star2->MassPrev() - m_Star2->Mass() + m_Star2->MassLossDiff() + m_Star2->MassTransferDiff(), 0.0); // update mass for star2 m_Star2->UpdateInitialMass(); // update effective initial mass of star 2 (MS, HG & HeMS) m_Star2->UpdateAgeAfterMassLoss(); // update age of star2 m_Star2->ApplyMassTransferRejuvenationFactor(); // apply age rejuvenation factor for star2 @@ -2263,7 +2208,7 @@ void BaseBinaryStar::ResolveMassChanges() { //Envelope ejection for convective envelope stars exceeding threshold luminosity to mass ratio: assume the entire envelope was lost on timescales long relative to the orbit if(m_Star1->EnvelopeJustExpelledByPulsations() || m_Star2->EnvelopeJustExpelledByPulsations()) { - m_SemiMajorAxis /= (2.0 - ((m_Star1->MassPrev() + m_Star2->MassPrev()) / (m_Star1->Mass() + m_Star2->Mass()))); // update separation in response to pulsational mass loss + m_SemiMajorAxis /= (2.0 - ((m_Star1->MassPrev() + m_Star2->MassPrev()) / (m_Star1->Mass() + m_Star2->Mass()))); // update separation in response to pulsational mass loss m_Star1->ResetEnvelopeExpulsationByPulsations(); m_Star2->ResetEnvelopeExpulsationByPulsations(); } @@ -2302,54 +2247,111 @@ void BaseBinaryStar::ResolveMassChanges() { */ void BaseBinaryStar::EvaluateBinary(const double p_Dt) { - CalculateMassTransfer(p_Dt); // calculate mass transfer if necessary + CalculateMassTransfer(p_Dt); // calculate mass transfer if necessary - (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::POST_MT); // print (log) detailed output + (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::POST_MT); // print (log) detailed output - CalculateWindsMassLoss(); // calculate mass loss dues to winds + CalculateWindsMassLoss(); // calculate mass loss dues to winds - (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::POST_WINDS); // print (log) detailed output + (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::POST_WINDS); // print (log) detailed output - if ((m_CEDetails.CEEnow || StellarMerger()) && // CEE or merger? - !(OPTIONS->CHEMode() != CHE_MODE::NONE && HasTwoOf({STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS}))) { // yes - avoid CEE if CH+CH + if ((m_CEDetails.CEEnow || StellarMerger()) && // CEE or merger? + !(OPTIONS->CHEMode() != CHE_MODE::NONE && HasTwoOf({STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS}))) { // yes - avoid CEE if CH+CH - ResolveCommonEnvelopeEvent(); // resolve CEE - immediate event - (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::POST_CEE); // print (log) detailed output + ResolveCommonEnvelopeEvent(); // resolve CEE - immediate event + (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::POST_CEE); // print (log) detailed output } else if (m_Star1->IsSNevent() || m_Star2->IsSNevent()) { - EvaluateSupernovae(); // evaluate supernovae (both stars) - immediate event - (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::POST_SN); // print (log) detailed output + EvaluateSupernovae(); // evaluate supernovae (both stars) - immediate event + (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::POST_SN); // print (log) detailed output if (HasOneOf({ STELLAR_TYPE::NEUTRON_STAR })) { - (void)PrintPulsarEvolutionParameters(PULSAR_RECORD_TYPE::DEFAULT); // print (log) pulsar evolution parameters + (void)PrintPulsarEvolutionParameters(PULSAR_RECORD_TYPE::DEFAULT); // print (log) pulsar evolution parameters } } else { - ResolveMassChanges(); // apply mass loss and mass transfer as necessary - (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::POST_MASS_RESOLUTION); // print (log) detailed output + ResolveMassChanges(); // apply mass loss and mass transfer as necessary + (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::POST_MASS_RESOLUTION); // print (log) detailed output - if (HasStarsTouching()) { // if stars emerged from mass transfer as touching, it's a merger + if (HasStarsTouching()) { // if stars emerged from mass transfer as touching, it's a merger m_Flags.stellarMerger = true; // Set Roche lobe flags for both stars so that they show correct RLOF status - m_Star1->SetRocheLobeFlags(m_CEDetails.CEEnow, m_SemiMajorAxis, m_Eccentricity); // set Roche lobe flags for star1 - m_Star2->SetRocheLobeFlags(m_CEDetails.CEEnow, m_SemiMajorAxis, m_Eccentricity); // set Roche lobe flags for star2 - (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::POST_MASS_RESOLUTION_MERGER); // print (log) detailed output + m_Star1->SetRocheLobeFlags(m_CEDetails.CEEnow, m_SemiMajorAxis, m_Eccentricity); // set Roche lobe flags for star1 + m_Star2->SetRocheLobeFlags(m_CEDetails.CEEnow, m_SemiMajorAxis, m_Eccentricity); // set Roche lobe flags for star2 + (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::POST_MASS_RESOLUTION_MERGER); // print (log) detailed output } } if ((m_Star1->IsSNevent() || m_Star2->IsSNevent())) { - EvaluateSupernovae(); // evaluate supernovae (both stars) if mass changes are responsible for a supernova - (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::POST_SN); // print (log) detailed output + EvaluateSupernovae(); // evaluate supernovae (both stars) if mass changes are responsible for a supernova + (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::POST_SN); // print (log) detailed output if (HasOneOf({ STELLAR_TYPE::NEUTRON_STAR })) { - (void)PrintPulsarEvolutionParameters(PULSAR_RECORD_TYPE::DEFAULT); // print (log) pulsar evolution parameters + (void)PrintPulsarEvolutionParameters(PULSAR_RECORD_TYPE::DEFAULT); // print (log) pulsar evolution parameters } } - // assign new values to "previous" values, for following timestep - m_EccentricityPrev = m_Eccentricity; - m_SemiMajorAxisPrev = m_SemiMajorAxis; + CalculateEnergyAndAngularMomentum(); // perform energy and angular momentum calculations + + if (OPTIONS->EnableTides() && !m_Unbound) { + + /* + std::cout << "\nTime = " << m_Time << "\n"; + std::cout << "Total angular momentum before = " << m_TotalAngularMomentum << "\n"; + std::cout << "Semi-major axis before = " << m_SemiMajorAxis << "\n"; + std::cout << "Eccentricity before = " << m_Eccentricity << "\n"; + std::cout << "Omega (star1) before = " << m_Star1->Omega() << "\n"; + std::cout << "Omega (star2) before = " << m_Star2->Omega() << "\n"; + std::cout << "Omega (binary) before = " << m_Omega << "\n"; + std::cout << "MoI (star1) before = " << m_Star1->CalculateMomentOfInertia() << "\n"; + std::cout << "MoI (star2) before = " << m_Star2->CalculateMomentOfInertia() << "\n"; + std::cout << "Mass (star1) before = " << m_Star1->Mass() << "\n"; + std::cout << "Mass (star2) before = " << m_Star2->Mass() << "\n"; + std::cout << "Radius (star1) before = " << m_Star1->Radius() << "\n"; + std::cout << "Radius (star2) before = " << m_Star2->Radius() << "\n"; + std::cout << "Gyration radius (star1) before = " << m_Star1->CalculateGyrationRadius() << "\n"; + std::cout << "Gyration radius (star2) before = " << m_Star2->CalculateGyrationRadius() << "\n"; + */ + + // find omega assuming synchronisation + // rough guess is (m_Star1->Omega() + m_Star2->Omega()) / 2.0 + m_Omega = OmegaAfterSynchronisation(m_Star1->Mass(), + m_Star2->Mass(), + m_Star1->CalculateMomentOfInertiaAU(), + m_Star2->CalculateMomentOfInertiaAU(), + m_TotalAngularMomentum, + (m_Star1->Omega() + m_Star2->Omega()) / 2.0); + + if (m_Omega > 0.0) { // root found? + // yes + m_Star1->SetOmega(m_Omega); // synchronise star 1 + m_Star2->SetOmega(m_Omega); // synchronise star 2 + + m_SemiMajorAxis = std::cbrt(G1) * std::cbrt((m_Star1->Mass() + m_Star2->Mass())) / PPOW(m_Omega, 2.0 / 3.0); // re-calculate semi-major axis + m_Eccentricity = 0.0; // circularise + m_TotalAngularMomentum = CalculateAngularMomentum(); // re-calculate total angular momentum + + // assign new values to "previous" values, for following timestep + m_EccentricityPrev = m_Eccentricity; + m_SemiMajorAxisPrev = m_SemiMajorAxis; + } - CalculateEnergyAndAngularMomentum(); // perform energy and angular momentum calculations + /* + std::cout << "Total angular momentum after = " << m_TotalAngularMomentum << "\n"; + std::cout << "Semi-major axis after = " << m_SemiMajorAxis << "\n"; + std::cout << "Eccentricity after = " << m_Eccentricity << "\n"; + std::cout << "Omega (star1) after = " << m_Star1->Omega() << "\n"; + std::cout << "Omega (star2) after = " << m_Star2->Omega() << "\n"; + std::cout << "Omega (binary) after = " << m_Omega << "\n"; + std::cout << "MoI (star1) after = " << m_Star1->CalculateMomentOfInertia() << "\n"; + std::cout << "MoI (star2) after = " << m_Star2->CalculateMomentOfInertia() << "\n"; + std::cout << "Mass (star1) after = " << m_Star1->Mass() << "\n"; + std::cout << "Mass (star2) after = " << m_Star2->Mass() << "\n"; + std::cout << "Radius (star1) after = " << m_Star1->Radius() << "\n"; + std::cout << "Radius (star2) after = " << m_Star2->Radius() << "\n"; + std::cout << "Gyration radius (star1) after = " << m_Star1->CalculateGyrationRadius() << "\n"; + std::cout << "Gyration radius (star2) after = " << m_Star2->CalculateGyrationRadius() << "\n"; + */ + } m_Star1->UpdateMagneticFieldAndSpin(m_CEDetails.CEEnow, m_Dt * MYR_TO_YEAR * SECONDS_IN_YEAR, EPSILON_PULSAR); // update pulsar parameters for star1 m_Star2->UpdateMagneticFieldAndSpin(m_CEDetails.CEEnow, m_Dt * MYR_TO_YEAR * SECONDS_IN_YEAR, EPSILON_PULSAR); // update pulsar parameters for star2 diff --git a/src/BaseBinaryStar.h b/src/BaseBinaryStar.h index 8c6cfa304..c0318033b 100644 --- a/src/BaseBinaryStar.h +++ b/src/BaseBinaryStar.h @@ -83,6 +83,7 @@ class BaseBinaryStar { m_MassTransferTrackerHistory = p_Star.m_MassTransferTrackerHistory; + m_Omega = p_Star.m_Omega; m_OrbitalVelocityPreSN = p_Star.m_OrbitalVelocityPreSN; m_RLOFDetails = p_Star.m_RLOFDetails; @@ -109,6 +110,7 @@ class BaseBinaryStar { m_TimePrev = p_Star.m_TimePrev; m_TimeToCoalescence = p_Star.m_TimeToCoalescence; + m_TotalAngularMomentumPrev = p_Star.m_TotalAngularMomentumPrev; m_TotalAngularMomentum = p_Star.m_TotalAngularMomentum; m_TotalEnergy = p_Star.m_TotalEnergy; @@ -216,10 +218,11 @@ class BaseBinaryStar { bool MassesEquilibratedAtBirth() const { return m_Flags.massesEquilibratedAtBirth; } MT_TRACKING MassTransferTrackerHistory() const { return m_MassTransferTrackerHistory; } bool MergesInHubbleTime() const { return m_Flags.mergesInHubbleTime; } + double Omega() const { return m_Omega; } bool OptimisticCommonEnvelope() const { return m_CEDetails.optimisticCE; } double OrbitalAngularVelocity() const { return std::sqrt(G1 * (m_Star1->Mass() + m_Star2->Mass()) / (m_SemiMajorAxis * m_SemiMajorAxis * m_SemiMajorAxis)); } // rads/year double OrbitalVelocityPreSN() const { return m_OrbitalVelocityPreSN; } - double Periastron() const { return m_SemiMajorAxis * (1.0-m_Eccentricity); } + double Periastron() const { return m_SemiMajorAxis * (1.0 - m_Eccentricity); } double PeriastronRsol() const { return Periastron() * AU_TO_RSOL; } double Radius1PostCEE() const { return m_Star1->RadiusPostCEE(); } double Radius2PostCEE() const { return m_Star2->RadiusPostCEE(); } @@ -348,6 +351,7 @@ class BaseBinaryStar { MT_TRACKING m_MassTransferTrackerHistory; + double m_Omega; // Orbital frequency double m_OrbitalVelocityPreSN; BinaryRLOFDetailsT m_RLOFDetails; // RLOF details @@ -374,6 +378,7 @@ class BaseBinaryStar { double m_DCOFormationTime; // Time of DCO formation double m_TotalAngularMomentum; + double m_TotalAngularMomentumPrev; double m_TotalEnergy; @@ -543,18 +548,16 @@ class BaseBinaryStar { //Functor for the boost root finder to determine how much mass needs to be lost from a donor without an envelope in order to fit inside the Roche lobe template - struct RadiusEqualsRocheLobeFunctor - { - RadiusEqualsRocheLobeFunctor(BaseBinaryStar *p_Binary, BinaryConstituentStar *p_Donor, BinaryConstituentStar *p_Accretor, ERROR *p_Error, double p_FractionAccreted) - { + struct RadiusEqualsRocheLobeFunctor { + RadiusEqualsRocheLobeFunctor(BaseBinaryStar *p_Binary, BinaryConstituentStar *p_Donor, BinaryConstituentStar *p_Accretor, ERROR *p_Error, double p_FractionAccreted) { m_Binary = p_Binary; m_Donor = p_Donor; m_Accretor = p_Accretor; m_Error = p_Error; m_FractionAccreted = p_FractionAccreted; } - T operator()(double const& dM) - { + T operator()(double const& dM) { + if (dM >= m_Donor->Mass()) { // Can't remove more than the donor's mass *m_Error = ERROR::TOO_MANY_RLOF_ITERATIONS; return m_Donor->Radius(); @@ -565,9 +568,9 @@ class BaseBinaryStar { BinaryConstituentStar* donorCopy = new BinaryConstituentStar(*m_Donor); double semiMajorAxis = m_Binary->CalculateMassTransferOrbit(donorCopy->Mass(), -dM , *m_Accretor, m_FractionAccreted); - double RLRadius = semiMajorAxis * (1 - m_Binary->Eccentricity()) * CalculateRocheLobeRadius_Static(donorMass - dM, accretorMass + (m_Binary->FractionAccreted() * dM)) * AU_TO_RSOL; + double RLRadius = semiMajorAxis * (1.0 - m_Binary->Eccentricity()) * CalculateRocheLobeRadius_Static(donorMass - dM, accretorMass + (m_Binary->FractionAccreted() * dM)) * AU_TO_RSOL; - (void)donorCopy->UpdateAttributes(-dM, -dM*donorCopy->Mass0()/donorCopy->Mass()); + (void)donorCopy->UpdateAttributes(-dM, -dM * donorCopy->Mass0() / donorCopy->Mass()); // Modify donor Mass0 and Age for MS (including HeMS) and HG stars donorCopy->UpdateInitialMass(); // update initial mass (MS, HG & HeMS) @@ -579,20 +582,20 @@ class BaseBinaryStar { delete donorCopy; donorCopy = nullptr; - return (RLRadius-thisRadiusAfterMassLoss); + return (RLRadius - thisRadiusAfterMassLoss); } private: - BaseBinaryStar *m_Binary; + BaseBinaryStar *m_Binary; BinaryConstituentStar *m_Donor; BinaryConstituentStar *m_Accretor; - ERROR *m_Error; - double m_FractionAccreted; + ERROR *m_Error; + double m_FractionAccreted; }; //Root solver to determine how much mass needs to be lost from a donor without an envelope in order to fit inside the Roche lobe - double MassLossToFitInsideRocheLobe(BaseBinaryStar * p_Binary, BinaryConstituentStar * p_Donor, BinaryConstituentStar * p_Accretor, double p_FractionAccreted) - { + double MassLossToFitInsideRocheLobe(BaseBinaryStar *p_Binary, BinaryConstituentStar *p_Donor, BinaryConstituentStar *p_Accretor, double p_FractionAccreted) { + using namespace std; // Help ADL of std functions. using namespace boost::math::tools; // For bracket_and_solve_root. @@ -612,7 +615,7 @@ class BaseBinaryStar { // iterations just thrash around. eps_tolerance tol(get_digits); // Set the tolerance. - std::pair root; + std::pair root(0.0, 0.0); try { ERROR error = ERROR::NONE; root = bracket_and_solve_root(RadiusEqualsRocheLobeFunctor(p_Binary, p_Donor, p_Accretor, &error, p_FractionAccreted), guess, factor, is_rising, tol, it); @@ -620,12 +623,73 @@ class BaseBinaryStar { } catch(exception& e) { SHOW_ERROR(ERROR::TOO_MANY_RLOF_ITERATIONS, e.what()); // Catch generic boost root finding error - m_Donor->Radius(); } - SHOW_WARN_IF(it>=maxit, ERROR::TOO_MANY_RLOF_ITERATIONS); + SHOW_WARN_IF(it >= maxit, ERROR::TOO_MANY_RLOF_ITERATIONS); - return root.first + (root.second - root.first)/2; // Midway between brackets is our result, if necessary we could return the result as an interval here. + return root.first + (root.second - root.first) / 2.0; // Midway between brackets is our result, if necessary we could return the result as an interval here. + } + + + /* + * Root solver to determine rotational frequency after synchronisation for tides + * + * Uses boost::math::tools::bracket_and_solve_root() + * + * + * double OmegaAfterSynchronisation(const double p_M1, const double p_M2, const double p_I1, const double p_I2, const double p_Omega) + * + * @param [IN] p_M1 Mass of star 1 + * @param [IN] p_M2 Mass of star 2 + * @param [IN] p_I1 Moment of intertia of star 1 + * @param [IN] p_I2 Moment of intertia of star 1 + * @param [IN] p_Ltot Total angular momentum for binary + * @param [IN] p_Guess Initial guess for value of root + * @return Root found: will be -1.0 if no real root found + */ + double OmegaAfterSynchronisation(const double p_M1, const double p_M2, const double p_I1, const double p_I2, const double p_Ltot, const double p_Guess) { + + const boost::uintmax_t maxit = TIDES_OMEGA_MAX_ITERATIONS; // maximum iterations + boost::uintmax_t it = maxit; // initially max iterations, but updated with actual count + int digits = std::numeric_limits::digits; // maximum possible binary digits accuracy + int get_digits = digits - 5; // we have to have a non-zero interval at each step + + // maximum accuracy is digits - 1. But we also have to allow for inaccuracy + // in the functor, otherwise the last few iterations just thrash around. + boost::math::tools::eps_tolerance tol(get_digits); // tolerance + + // define functor + // function: ax + bx^(1/3) + c = 0 + double a = p_I1 + p_I2; + double b = PPOW(G1, 2.0 / 3.0) * p_M1 * p_M2 / std::cbrt(p_M1 + p_M2); + double c = -p_Ltot; + + auto func = [this, a, b, c](double x) -> double { return (a * x) + (b / std::cbrt(x)) + c; }; // functor + + // find root + double factor = TIDES_OMEGA_SEARCH_FACTOR; // size of search steps + bool is_rising = func(p_Guess) > func(p_Guess * factor) ? false : true; // so bracket_and_solve_root() knows whether to increase or decrease guess per iteration + + std::pair root(-1.0, -1.0); // initialise root + try { + root = boost::math::tools::bracket_and_solve_root(func, p_Guess, factor, is_rising, tol, it); // iterate to find root + } + catch(std::exception& e) { // catch generic boost root finding error + root.first = -1.0; // set error return + root.second = -1.0; + if (it < maxit) { // too many iterations? + SHOW_ERROR(ERROR::ROOT_FINDER_FAILED, e.what()); // no - some other error - show it + } + } + + if (it >= maxit) { // too many iterations? + root.first = -1.0; // yes - set error return + root.second = -1.0; + SHOW_WARN(ERROR::TOO_MANY_OMEGA_ITERATIONS); // show warning + } + + return root.first + (root.second - root.first) / 2.0; // midway between brackets (could return brackets...) } + }; #endif // __BaseBinaryStar_h__ diff --git a/src/BaseStar.cpp b/src/BaseStar.cpp index 1c8489d01..616d05c8a 100755 --- a/src/BaseStar.cpp +++ b/src/BaseStar.cpp @@ -111,6 +111,7 @@ BaseStar::BaseStar(const unsigned long int p_RandomSeed, m_TZAMS = CalculateTemperatureOnPhase_Static(m_LZAMS, m_RZAMS); m_OmegaCHE = CalculateOmegaCHE(m_MZAMS, m_Metallicity); + m_OmegaZAMS = p_RotationalFrequency >= 0.0 // valid rotational frequency passed in? ? p_RotationalFrequency // yes - use it : CalculateZAMSAngularFrequency(m_MZAMS, m_RZAMS); // no - calculate it @@ -140,6 +141,7 @@ BaseStar::BaseStar(const unsigned long int p_RandomSeed, m_DominantMassLossRate = MASS_LOSS_TYPE::NONE; m_Omega = m_OmegaZAMS; + m_AngularMomentum = DEFAULT_INITIAL_DOUBLE_VALUE; m_MinimumLuminosityOnPhase = DEFAULT_INITIAL_DOUBLE_VALUE; m_LBVphaseFlag = false; @@ -151,7 +153,7 @@ BaseStar::BaseStar(const unsigned long int p_RandomSeed, m_RadiusPrev = m_RZAMS; m_DtPrev = DEFAULT_INITIAL_DOUBLE_VALUE; m_OmegaPrev = m_OmegaZAMS; - + // Lambdas m_Lambdas.dewi = DEFAULT_INITIAL_DOUBLE_VALUE; m_Lambdas.fixed = DEFAULT_INITIAL_DOUBLE_VALUE; @@ -2444,7 +2446,7 @@ double BaseStar::CalculateRotationalVelocity(double p_MZAMS) const { break; default: // unknown rorational velocity prescription - SHOW_WARN(ERROR::UNKNOWN_VROT_PRESCRIPTION, "Using default vRot = 0.0"); // show warning + SHOW_WARN(ERROR::UNKNOWN_VROT_PRESCRIPTION, "Using default vRot = 0.0"); // show warning } return vRot; } @@ -2457,7 +2459,7 @@ double BaseStar::CalculateRotationalVelocity(double p_MZAMS) const { * Hurley et al. 2000, eq 108 * * - * double CalculateRotationalAngularFrequency(const double p_MZAMS, const double p_RZAMS) + * double CalculateZAMSAngularFrequency(const double p_MZAMS, const double p_RZAMS) * * @param [IN] p_MZAMS Zero age main sequence mass in Msol * @param [IN] p_RZAMS Zero age main sequence radius in Rsol @@ -2465,7 +2467,7 @@ double BaseStar::CalculateRotationalVelocity(double p_MZAMS) const { */ double BaseStar::CalculateZAMSAngularFrequency(const double p_MZAMS, const double p_RZAMS) const { double vRot = CalculateRotationalVelocity(p_MZAMS); - return utils::Compare(vRot, 0.0) == 0 ? 0.0 : 45.35 * vRot / p_RZAMS; // Hurley et al. 2000, eq 108 JR: todo: added check for vRot = 0 + return utils::Compare(vRot, 0.0) == 0 ? 0.0 : 45.35 * vRot / p_RZAMS; // Hurley et al. 2000, eq 108 } @@ -2516,7 +2518,7 @@ double BaseStar::CalculateOmegaCHE(const double p_MZAMS, const double p_Metallic } // calculate omegaCHE(M, Z) - return (1.0 / ((0.09 * log(p_Metallicity / 0.004)) + 1.0) * omegaZ004) * SECONDS_IN_YEAR; // in rads/yr + return (1.0 / ((0.09 * log(p_Metallicity / 0.004)) + 1.0) * omegaZ004) * SECONDS_IN_YEAR; // in rads/yr #undef massCutoffs } @@ -3242,10 +3244,10 @@ void BaseStar::UpdateAttributesAndAgeOneTimestepPreamble(const double p_DeltaMas // record some current values before they are (possibly) changed by evolution if (p_DeltaTime > 0.0) { // don't use utils::Compare() here - m_StellarTypePrev = m_StellarType; - m_DtPrev = m_Dt; - m_MassPrev = m_Mass; - m_RadiusPrev = m_Radius; + m_StellarTypePrev = m_StellarType; + m_DtPrev = m_Dt; + m_MassPrev = m_Mass; + m_RadiusPrev = m_Radius; } // the GBParams and Timescale calculations need to be done @@ -3438,21 +3440,21 @@ STELLAR_TYPE BaseStar::EvolveOnPhase() { if (ShouldEvolveOnPhase()) { // Evolve timestep on phase - m_Tau = CalculateTauOnPhase(); + m_Tau = CalculateTauOnPhase(); - m_COCoreMass = CalculateCOCoreMassOnPhase(); - m_CoreMass = CalculateCoreMassOnPhase(); - m_HeCoreMass = CalculateHeCoreMassOnPhase(); + m_COCoreMass = CalculateCOCoreMassOnPhase(); + m_CoreMass = CalculateCoreMassOnPhase(); + m_HeCoreMass = CalculateHeCoreMassOnPhase(); - m_Luminosity = CalculateLuminosityOnPhase(); + m_Luminosity = CalculateLuminosityOnPhase(); std::tie(m_Radius, stellarType) = CalculateRadiusAndStellarTypeOnPhase(); // Radius and possibly new stellar type - m_Mu = CalculatePerturbationMuOnPhase(); + m_Mu = CalculatePerturbationMuOnPhase(); PerturbLuminosityAndRadiusOnPhase(); - m_Temperature = CalculateTemperatureOnPhase(); + m_Temperature = CalculateTemperatureOnPhase(); STELLAR_TYPE thisStellarType = ResolveEnvelopeLoss(); // Resolve envelope loss if it occurs - possibly new stellar type if (thisStellarType != m_StellarType) { // thisStellarType overrides stellarType (from CalculateRadiusAndStellarTypeOnPhase()) diff --git a/src/BaseStar.h b/src/BaseStar.h index 45dfe0a8c..a9304cee1 100644 --- a/src/BaseStar.h +++ b/src/BaseStar.h @@ -195,7 +195,7 @@ class BaseStar { double CalculateOmegaCHE(const double p_MZAMS, const double p_Metallicity) const; - double CalculateRadialChange() const { return (utils::Compare(m_RadiusPrev,0)<=0)? 0 : std::abs(m_Radius - m_RadiusPrev) / m_RadiusPrev; } // Return fractional radial change (if previous radius is negative or zero, return 0 to avoid NaN + double CalculateRadialChange() const { return (utils::Compare(m_RadiusPrev,0)<=0)? 0 : std::abs(m_Radius - m_RadiusPrev) / m_RadiusPrev; } // Return fractional radial change (if previous radius is negative or zero, return 0 to avoid NaN double CalculateRadialExpansionTimescale() const { return CalculateRadialExpansionTimescale_Static(m_StellarType, m_StellarTypePrev, m_Radius, m_RadiusPrev, m_DtPrev); } // Use class member variables @@ -214,13 +214,13 @@ class BaseStar { double CalculateTimestep(); double CalculateZetaAdiabatic(); - virtual double CalculateZetaConstantsByEnvelope(ZETA_PRESCRIPTION p_ZetaPrescription) { return 0.0; } // Use inheritance hierarchy + virtual double CalculateZetaConstantsByEnvelope(ZETA_PRESCRIPTION p_ZetaPrescription) { return 0.0; } // Use inheritance hierarchy void ClearCurrentSNEvent() { m_SupernovaDetails.events.current = SN_EVENT::NONE; } // Clear supernova event/state for current timestep virtual ENVELOPE DetermineEnvelopeType() const { return ENVELOPE::REMNANT; } // Default is REMNANT - but should never be called - void HaltWinds() { m_Mdot = 0.0; } // Disable wind mass loss in current time step (e.g., if star is a donor or accretor in a RLOF episode) + void HaltWinds() { m_Mdot = 0.0; } // Disable wind mass loss in current time step (e.g., if star is a donor or accretor in a RLOF episode) virtual ACCRETION_REGIME DetermineAccretionRegime(const bool p_HeRich, const double p_DonorThermalMassLossRate) { return ACCRETION_REGIME::NONE; } // Placeholder, use inheritance for WDs @@ -242,12 +242,12 @@ class BaseStar { void SetStellarTypePrev(const STELLAR_TYPE p_StellarTypePrev) { m_StellarTypePrev = p_StellarTypePrev; } - bool ShouldEnvelopeBeExpelledByPulsations() const { return false; } // Default is that there is no envelope expulsion by pulsations + bool ShouldEnvelopeBeExpelledByPulsations() const { return false; } // Default is that there is no envelope expulsion by pulsations void StashSupernovaDetails(const STELLAR_TYPE p_StellarType, const SSE_SN_RECORD_TYPE p_RecordType = SSE_SN_RECORD_TYPE::DEFAULT) { LOGGING->StashSSESupernovaDetails(this, p_StellarType, p_RecordType); } - virtual void UpdateAgeAfterMassLoss() { } // Default is NO-OP + virtual void UpdateAgeAfterMassLoss() { } // Default is NO-OP STELLAR_TYPE UpdateAttributesAndAgeOneTimestep(const double p_DeltaMass, const double p_DeltaMass0, @@ -313,6 +313,7 @@ class BaseStar { double m_RZAMS; // ZAMS Radius double m_TZAMS; // ZAMS Temperature + // Effective Zero Age Main Sequence double m_LZAMS0; // Effective ZAMS Luminosity double m_RZAMS0; // Effective ZAMS Radius @@ -334,7 +335,8 @@ class BaseStar { double m_Mdot; // Current mass loss rate (Msol per ?) MASS_LOSS_TYPE m_DominantMassLossRate; // Current dominant mass loss rate double m_Mu; // Current small envelope parameter mu - double m_Omega; // Current angular frequency (yr-1) + double m_Omega; // Current angular frequency (yr^-1) + double m_AngularMomentum; // Current angular momentum double m_Radius; // Current radius (Rsol) double m_Tau; // Relative time double m_Temperature; // Current temperature (Tsol) @@ -343,7 +345,7 @@ class BaseStar { // Previous timestep variables double m_DtPrev; // Previous timestep double m_MassPrev; // Previous mass (Msol) - double m_OmegaPrev; // Previous angular frequency (yr-1) + double m_OmegaPrev; // Previous angular frequency (yr^-1) double m_RadiusPrev; // Previous radius (Rsol) STELLAR_TYPE m_StellarTypePrev; // Stellar type at previous timestep @@ -417,8 +419,6 @@ class BaseStar { DBL_VECTOR &p_RConstants, DBL_VECTOR &p_GammaConstants); - virtual void CalculateAndSetPulsarParameters() { } // NO-OP for most stellar types - double CalculateBindingEnergy(const double p_CoreMass, const double p_EnvMass, const double p_Radius, const double p_Lambda) const; void CalculateBnCoefficients(DBL_VECTOR &p_BnCoefficients); diff --git a/src/COWD.h b/src/COWD.h index dcebf257d..c0f799d7c 100755 --- a/src/COWD.h +++ b/src/COWD.h @@ -26,7 +26,6 @@ class COWD: virtual public BaseStar, public WhiteDwarfs { // member functions - void CalculateAngularMomentum() const { } // NO-OP static double CalculateLuminosityOnPhase_Static(const double p_Mass, const double p_Time, diff --git a/src/HG.h b/src/HG.h index 663a75792..53a92a201 100755 --- a/src/HG.h +++ b/src/HG.h @@ -127,20 +127,20 @@ class HG: virtual public BaseStar, public GiantBranch { { Mass0YieldsDesiredCoreMassFunctor(HG *p_Star, double p_DesiredCoreMass, ERROR *p_Error) { - m_Star = p_Star; - m_DesiredCoreMass = p_DesiredCoreMass; - m_Error = p_Error; + m_Star = p_Star; + m_DesiredCoreMass = p_DesiredCoreMass; + m_Error = p_Error; } T operator()(double const& guessMass0) { - HG * copy = new HG(*m_Star, false); + HG *copy = new HG(*m_Star, false); copy->UpdateAttributesAndAgeOneTimestep(0.0, guessMass0 - copy->Mass0(), 0.0, true); - double coreMassEstimate=copy->CalculateCoreMassOnPhase(guessMass0, copy->Age()); + double coreMassEstimate = copy->CalculateCoreMassOnPhase(guessMass0, copy->Age()); delete copy; copy = nullptr; return (coreMassEstimate - m_DesiredCoreMass); } private: - HG *m_Star; + HG *m_Star; double m_DesiredCoreMass; ERROR *m_Error; }; @@ -155,7 +155,7 @@ class HG: virtual public BaseStar, public GiantBranch { double guess = p_Star->Mass(); // Rough guess at solution double factor = ADAPTIVE_MASS0_SEARCH_FACTOR; // Size of search steps - const boost::uintmax_t maxit = ADAPTIVE_MASS0_MAX_ITERATIONS; // Limit to maximum iterations. + const boost::uintmax_t maxit = ADAPTIVE_MASS0_MAX_ITERATIONS; // Limit to maximum iterations. boost::uintmax_t it = maxit; // Initially our chosen max iterations, but updated with actual. bool is_rising = true; // So if result with guess is too low, then try increasing guess. int digits = std::numeric_limits::digits; // Maximum possible binary digits accuracy for type T. @@ -175,11 +175,11 @@ class HG: virtual public BaseStar, public GiantBranch { if (error != ERROR::NONE) SHOW_WARN(error); } catch(exception& e) { - SHOW_ERROR(ERROR::TOO_MANY_MASS0_ITERATIONS, e.what()); // Catch generic boost root finding error + SHOW_ERROR(ERROR::TOO_MANY_MASS0_ITERATIONS, e.what()); // Catch generic boost root finding error } SHOW_WARN_IF(it>=maxit, ERROR::TOO_MANY_MASS0_ITERATIONS); - return root.first + (root.second - root.first)/2; // Midway between brackets is our result, if necessary we could return the result as an interval here. + return root.first + (root.second - root.first) / 2.0; // Midway between brackets is our result, if necessary we could return the result as an interval here. } }; diff --git a/src/HeHG.h b/src/HeHG.h index 25fe27531..389e007a3 100755 --- a/src/HeHG.h +++ b/src/HeHG.h @@ -72,7 +72,7 @@ class HeHG: virtual public BaseStar, public HeMS { double CalculateHeCoreMassOnPhase() const { return m_Mass; } // NO-OP double CalculateLambdaNanjingStarTrack(const double p_Mass, const double p_Metallicity) const; - double CalculateLambdaNanjingEnhanced(const int p_MassInd, const int p_Zind) const { return CalculateLambdaNanjingStarTrack(0.0, 0.0); } // 0.0 are dummy values that are not used + double CalculateLambdaNanjingEnhanced(const int p_MassInd, const int p_Zind) const { return CalculateLambdaNanjingStarTrack(0.0, 0.0); } // 0.0 are dummy values that are not used double CalculateLuminosityOnPhase() const; double CalculateLuminosityAtPhaseEnd() const { return m_Luminosity; } // NO-OP @@ -104,7 +104,7 @@ class HeHG: virtual public BaseStar, public HeMS { void CalculateTimescales(const double p_Mass, DBL_VECTOR &p_Timescales); void CalculateTimescales() { CalculateTimescales(m_Mass0, m_Timescales); } // Use class member variables - double CalculateZetaConstantsByEnvelope(ZETA_PRESCRIPTION p_ZetaPrescription) { return GiantBranch::CalculateZetaConstantsByEnvelope(p_ZetaPrescription); } // Calculate Zetas as for other giant stars (HeMS stars were an exception) + double CalculateZetaConstantsByEnvelope(ZETA_PRESCRIPTION p_ZetaPrescription) { return GiantBranch::CalculateZetaConstantsByEnvelope(p_ZetaPrescription); } // Calculate Zetas as for other giant stars (HeMS stars were an exception) double ChooseTimestep(const double p_Time) const; @@ -122,7 +122,7 @@ class HeHG: virtual public BaseStar, public HeMS { void ResolveHeliumFlash() { } // NO-OP STELLAR_TYPE ResolveSkippedPhase() { return m_StellarType; } // NO-OP - bool ShouldEnvelopeBeExpelledByPulsations() const { return CHeB::ShouldEnvelopeBeExpelledByPulsations(); } // Envelope of convective star with luminosity to mass ratio beyond threshold should be expelled + bool ShouldEnvelopeBeExpelledByPulsations() const { return CHeB::ShouldEnvelopeBeExpelledByPulsations(); } // Envelope of convective star with luminosity to mass ratio beyond threshold should be expelled bool ShouldEvolveOnPhase() const; bool ShouldSkipPhase() const { return false; } // Never skip HeMS phase diff --git a/src/MainSequence.cpp b/src/MainSequence.cpp index abdfa7218..4d5a04a10 100644 --- a/src/MainSequence.cpp +++ b/src/MainSequence.cpp @@ -598,7 +598,7 @@ double MainSequence::CalculateGyrationRadius() const { if ((utils::Compare(log10M, 0.2) > 0)) CUpper = -1.5; // log10(M) > 0.2 else if ((utils::Compare(log10M, 0.0) > 0)) CUpper = -2.5 + (5.0 * log10M); // 0.2 <= log10(M) > 0.0 (de Mink doesn't include '=' - we assume it to be here (and for log10(M) <= 0.0)) - double k0 = cLower+ std::min(0.21, std::max(0.09 - (0.27 * log10M), 0.037 + (0.033 * log10M))); // gyration radius squared for ZAMS stars + double k0 = cLower + std::min(0.21, std::max(0.09 - (0.27 * log10M), 0.037 + (0.033 * log10M)));// gyration radius squared for ZAMS stars double radiusRatio = m_Radius / m_RZAMS; diff --git a/src/NS.cpp b/src/NS.cpp index 0a04aaf84..0b45843d3 100755 --- a/src/NS.cpp +++ b/src/NS.cpp @@ -1,5 +1,4 @@ #include "Rand.h" - #include "NS.h" @@ -8,6 +7,8 @@ * * Hurley et al. 2000, eq 93 * + * Called (indirectly) from GiantBranch, so must be static. + * * * double CalculateLuminosityOnPhase_Static(const double p_Mass, const double p_Time) * @@ -17,7 +18,7 @@ */ double NS::CalculateLuminosityOnPhase_Static(const double p_Mass, const double p_Time) { double t = std::max(p_Time, 0.1); - return 0.02 * PPOW(p_Mass, 2.0/3.0) / (t * t); + return 0.02 * PPOW(p_Mass, 2.0 / 3.0) / (t * t); } @@ -25,10 +26,8 @@ double NS::CalculateLuminosityOnPhase_Static(const double p_Mass, const double p * Choose timestep for Pulsar Evolution * * Pulsars evolve very fast when they are first born, and evolve slower as they age. - * Hence, timestep is chosen to be small when pulsar is young, - * and is slowly increased as the pulsar ages. - * - * Can change it to a choice that suits your simulation. + * Hence, timestep is chosen to be small when pulsar is young, and is slowly increased + * as the pulsar ages. * * double ChooseTimestep(const double p_Time) * @@ -36,7 +35,8 @@ double NS::CalculateLuminosityOnPhase_Static(const double p_Mass, const double p * @return Suggested timestep (dt) */ double NS::ChooseTimestep(const double p_Time) const { - double result; + double result = 500.0; // default value + if (p_Time < 0.01) { result = 0.001; } @@ -50,13 +50,11 @@ double NS::ChooseTimestep(const double p_Time) const { result = 1.0; } else if (p_Time < 500.0) { - double slope = log10(500.0) / (log10(500.0) - 1.0); + double slope = 1.58859191006; // 1.58859191006 = log10(500.0) / (log10(500.0) - 1.0) double log10_step = slope * (log10(p_Time) - 1.0); result = PPOW(10.0, log10_step); } - else { - result = 500.0; - } + return result; } @@ -74,13 +72,13 @@ double NS::CalculateRadiusOnPhaseInKM_Static(const double p_Mass) { double radius; - switch (OPTIONS->NeutronStarEquationOfState()) { // which equation-of-state? + switch (OPTIONS->NeutronStarEquationOfState()) { // which equation-of-state? - case NS_EOS::SSE: // SSE + case NS_EOS::SSE: // SSE radius = 10.0; break; - case NS_EOS::ARP3: { // ARP3 + case NS_EOS::ARP3: { // ARP3 // We don't extrapolate so masses outside table just set to extreme values @@ -105,13 +103,14 @@ double NS::CalculateRadiusOnPhaseInKM_Static(const double p_Mass) { } } break; - default: // unknown equation-of-state - SHOW_WARN_STATIC(ERROR::UNKNOWN_NS_EOS, // show warning + default: // unknown equation-of-state + SHOW_WARN_STATIC(ERROR::UNKNOWN_NS_EOS, // show warning "Using default NS radius = 10.0", OBJECT_TYPE::BASE_STAR, STELLAR_TYPE::NEUTRON_STAR); radius = 10.0; } + return radius; } @@ -119,16 +118,18 @@ double NS::CalculateRadiusOnPhaseInKM_Static(const double p_Mass) { /* * Calculate core collapse Supernova parameters * - * + * Called from GiantBranch, so must be static. + * + * * DBL_DBL_DBL CalculateCoreCollapseSNParams_Static(const double p_Mass) * * @param [IN] p_Mass Mass in Msol * @return Tuple containing Luminosity, Radius and Temperature of Neutron Star */ DBL_DBL_DBL NS::CalculateCoreCollapseSNParams_Static(const double p_Mass) { - double luminosity = CalculateLuminosityOnPhase_Static(p_Mass, 0.0); // Luminosity of Neutron Star as it cools - double radius = CalculateRadiusOnPhase_Static(p_Mass); // Radius of Neutron Star in Rsol - double temperature = BaseStar::CalculateTemperatureOnPhase_Static(luminosity, radius); // Temperature of NS + double luminosity = CalculateLuminosityOnPhase_Static(p_Mass, 0.0); // Luminosity of Neutron Star as it cools + double radius = CalculateRadiusOnPhase_Static(p_Mass); // Radius of Neutron Star in Rsol + double temperature = BaseStar::CalculateTemperatureOnPhase_Static(luminosity, radius); // Temperature of NS return std::make_tuple(luminosity, radius, temperature); } @@ -138,30 +139,30 @@ DBL_DBL_DBL NS::CalculateCoreCollapseSNParams_Static(const double p_Mass) { * Calculate the spin period of a Pulsar at birth according to selected distribution (by commandline option) * * - * double CalculatePulsarBirthSpinPeriod_Static() + * double CalculatePulsarBirthSpinPeriod() * * @return Birth spin period of Pulsar in ms */ -double NS::CalculatePulsarBirthSpinPeriod_Static() { +double NS::CalculatePulsarBirthSpinPeriod() { double pSpin; - switch (OPTIONS->PulsarBirthSpinPeriodDistribution()) { // which distribution? + switch (OPTIONS->PulsarBirthSpinPeriodDistribution()) { // which distribution? - case PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION::ZERO: // ZERO + case PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION::ZERO: // ZERO pSpin = 0.0; break; - case PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION::FIXED: // FIXED constant value as used in default model in Oslowski et al 2011 https://arxiv.org/abs/0903.3538 - SHOW_WARN_STATIC(ERROR::UNSUPPORTED_PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION, // show warning + case PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION::FIXED: // FIXED constant value as used in default model in Oslowski et al 2011 https://arxiv.org/abs/0903.3538 + SHOW_WARN_STATIC(ERROR::UNSUPPORTED_PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION, // show warning "Using spin = 0.0", OBJECT_TYPE::BASE_STAR, STELLAR_TYPE::NEUTRON_STAR); pSpin = 0.0; break; - case PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION::UNIFORM: { // UNIFORM distribution between minimum and maximum value as in Oslowski et al 2011 https://arxiv.org/abs/0903.3538 (default Pmin = and Pmax = ) - // and also Kiel et al 2008 https://arxiv.org/abs/0805.0059 (default Pmin = 10 ms and Pmax 100 ms, section 3.4) + case PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION::UNIFORM: { // UNIFORM distribution between minimum and maximum value as in Oslowski et al 2011 https://arxiv.org/abs/0903.3538 (default Pmin = and Pmax = ) + // and also Kiel et al 2008 https://arxiv.org/abs/0805.0059 (default Pmin = 10 ms and Pmax 100 ms, section 3.4) double maximum = OPTIONS->PulsarBirthSpinPeriodDistributionMax(); double minimum = OPTIONS->PulsarBirthSpinPeriodDistributionMin(); @@ -169,7 +170,7 @@ double NS::CalculatePulsarBirthSpinPeriod_Static() { pSpin = minimum + (RAND->Random() * (maximum - minimum)); } break; - case PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION::NORMAL: { // NORMAL distribution from Faucher-Giguere and Kaspi 2006 https://arxiv.org/abs/astro-ph/0512585 + case PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION::NORMAL: { // NORMAL distribution from Faucher-Giguere and Kaspi 2006 https://arxiv.org/abs/astro-ph/0512585 // Values hard-coded for now, can make them options if necessary // pulsarBirthSpinPeriodDistributionFaucherGiguereKaspi2006Mean = 300.0; @@ -182,8 +183,8 @@ double NS::CalculatePulsarBirthSpinPeriod_Static() { } break; - default: // unknown distribution - SHOW_WARN_STATIC(ERROR::UNKNOWN_PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION, // show warning + default: // unknown distribution + SHOW_WARN_STATIC(ERROR::UNKNOWN_PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION, // show warning "Using spin = 0.0", OBJECT_TYPE::BASE_STAR, STELLAR_TYPE::NEUTRON_STAR); @@ -199,29 +200,29 @@ double NS::CalculatePulsarBirthSpinPeriod_Static() { * according to selected distribution (by commandline option) * * - * double CalculatePulsarBirthMagneticField_Static() + * double CalculatePulsarBirthMagneticField() * * @return log10 of the birth magnetic field in G */ -double NS::CalculatePulsarBirthMagneticField_Static() { +double NS::CalculatePulsarBirthMagneticField() { double log10B; - switch (OPTIONS->PulsarBirthMagneticFieldDistribution()) { // which distribution? + switch (OPTIONS->PulsarBirthMagneticFieldDistribution()) { // which distribution? - case PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION::ZERO: // ZERO + case PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION::ZERO: // ZERO log10B = 0.0; break; - case PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION::FIXED: // FIXED - set to a fixed constant value - SHOW_WARN_STATIC(ERROR::UNSUPPORTED_PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION, // show warning + case PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION::FIXED: // FIXED - set to a fixed constant value + SHOW_WARN_STATIC(ERROR::UNSUPPORTED_PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION, // show warning "Using 0.0", OBJECT_TYPE::BASE_STAR, STELLAR_TYPE::NEUTRON_STAR); log10B = 0.0; break; - case PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION::FLATINLOG: { // FLAT IN LOG distribution from Oslowski et al 2011 https://arxiv.org/abs/0903.3538 (log10B0min = , log10B0max = ) + case PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION::FLATINLOG: { // FLAT IN LOG distribution from Oslowski et al 2011 https://arxiv.org/abs/0903.3538 (log10B0min = , log10B0max = ) double maximum = OPTIONS->PulsarBirthMagneticFieldDistributionMax(); double minimum = OPTIONS->PulsarBirthMagneticFieldDistributionMin(); @@ -230,7 +231,7 @@ double NS::CalculatePulsarBirthMagneticField_Static() { } break; - case PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION::UNIFORM: { // UNIFORM flat distribution used in Kiel et al 2008 https://arxiv.org/abs/0805.0059 (log10B0min = 11, log10B0max = 13.5 see section 3.4 and Table 1.) + case PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION::UNIFORM: { // UNIFORM flat distribution used in Kiel et al 2008 https://arxiv.org/abs/0805.0059 (log10B0min = 11, log10B0max = 13.5 see section 3.4 and Table 1.) double maximum = PPOW(10.0, OPTIONS->PulsarBirthMagneticFieldDistributionMax()); @@ -239,7 +240,7 @@ double NS::CalculatePulsarBirthMagneticField_Static() { log10B = log10(minimum + (RAND->Random() * (maximum - minimum))); } break; - case PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION::LOGNORMAL: { // LOG NORMAL distribution from Faucher-Giguere and Kaspi 2006 https://arxiv.org/abs/astro-ph/0512585 + case PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION::LOGNORMAL: { // LOG NORMAL distribution from Faucher-Giguere and Kaspi 2006 https://arxiv.org/abs/astro-ph/0512585 // Values hard-coded for now, can make them options if necessary // pulsarBirthMagneticFieldDistributionFaucherGiguereKaspi2006Mean = 12.65 @@ -251,8 +252,8 @@ double NS::CalculatePulsarBirthMagneticField_Static() { log10B = RAND->RandomGaussian(sigma) + mean; } break; - default: // unknown distribution - SHOW_WARN_STATIC(ERROR::UNKNOWN_PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION, // show warning + default: // unknown distribution + SHOW_WARN_STATIC(ERROR::UNKNOWN_PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION, // show warning "Using 0.0", OBJECT_TYPE::BASE_STAR, STELLAR_TYPE::NEUTRON_STAR); @@ -265,27 +266,29 @@ double NS::CalculatePulsarBirthMagneticField_Static() { /* * Calculate the moment of inertia for a Neutron Star using a model independent relation between - * the moment of inertia, mass and radius of a neutron star + * the moment of inertia, mass and radius of a neutron star - return MoI in CGS. + * + * Uses m_Mass and m_Radius to calculate moment of inertia. * * Raithel et al. 2016, eq 8 in https://arxiv.org/abs/1603.06594 * https://tap.arizona.edu/sites/tap.arizona.edu/files/Raithel_2015_manuscript.pdf * * - * double CalculateMomentOfInertia_Static(const double p_Mass, const double p_Radius) + * double CalculateMomentOfInertiaCGS() * - * @param [IN] p_Mass Mass in Msol - * @param [IN] p_Radius Radius in km * @return Moment of inertia in g cm^2 */ -double NS::CalculateMomentOfInertia_Static(const double p_Mass, const double p_Radius) { +double NS::CalculateMomentOfInertiaCGS() const { // pow() is slow - use multiplication - double m_r = p_Mass / p_Radius; + + double radius = m_Radius * RSOL_TO_KM; + double m_r = m_Mass / radius; double m_r_4 = m_r * m_r * m_r * m_r; - double r_km = p_Radius * KM_TO_CM; - double r_km_2 = r_km * r_km; + double r_cm = m_Radius * KM_TO_CM; + double r_cm_2 = r_cm * r_cm; - return 0.237 * p_Mass * MSOL_TO_G * r_km_2 * (1.0 + (4.2 * m_r) + 90.0 * m_r_4); + return 0.237 * m_Mass * MSOL_TO_G * r_cm_2 * (1.0 + (4.2 * m_r) + 90.0 * m_r_4); } @@ -297,7 +300,7 @@ double NS::CalculateMomentOfInertia_Static(const double p_Mass, const double p_R * This is changed to the form of calculating spindown with P and Pdot, then convert to OmegaDot and to be recorded in the output file. * Evolution of the inclination between pulsar magnetic and rotational axes will be considered in a future version. * - * double CalculateSpinDownRate_Static(const double p_Omega, const double p_MomentOfInteria, const double p_MagField, const double p_Radius) + * double CalculateSpinDownRate(const double p_Omega, const double p_MomentOfInteria, const double p_MagField, const double p_Radius) * * @param [IN] p_Omega Pulsar spin frequency. * @param [IN] p_MomentOfInteria Moment of Interia of the Neutron Star in kg m^2 @@ -305,44 +308,55 @@ double NS::CalculateMomentOfInertia_Static(const double p_Mass, const double p_R * @param [IN] p_Radius Radius of the Neutron Star in metres * @return Spin down rate (spin frequency derivative) of an isolated Neutron Star in s^(-2) */ -double NS::CalculateSpinDownRate_Static(const double p_Omega, const double p_MomentOfInteria, const double p_MagField, const double p_Radius) { +double NS::CalculateSpinDownRate(const double p_Omega, const double p_MomentOfInteria, const double p_MagField, const double p_Radius) const { // pow() is slow - use multiplication - double Period = 2 * M_PI / p_Omega ; //Convert frequency to period - double cgs_Radius = p_Radius * KM_TO_CM; // radius in cm - double radius_6 = cgs_Radius * cgs_Radius * cgs_Radius * cgs_Radius * cgs_Radius * cgs_Radius; - double cgs_MagField = p_MagField * TESLA_TO_GAUSS; //B field in G - double magField_2 = cgs_MagField * cgs_MagField; - constexpr double _8_PI_2 = 8.0 * M_PI * M_PI; - constexpr double _3_C_3 = 3.0 * (C * 100.0) * (C * 100.0) * (C * 100.0) ; - double pDotTop = _8_PI_2 * radius_6 * magField_2; - double pDotBottom = _3_C_3 * p_MomentOfInteria * Period ; - double pDot = pDotTop / pDotBottom ; //Period derivative - return( -1 * pDot * p_Omega / Period); //convert period derivative to frequency derivative, which is what is recorded in the output. - + double period = _2_PI / p_Omega; // convert frequency to period + double cgsRadius = p_Radius * KM_TO_CM; // radius in cm + double radius_6 = cgsRadius * cgsRadius * cgsRadius * cgsRadius * cgsRadius * cgsRadius; + double cgsMagField = p_MagField * TESLA_TO_GAUSS; // B field in G + double magField_2 = cgsMagField * cgsMagField; + constexpr double _8_PI_2 = 8.0 * PI_2; + constexpr double _3_C_3 = 3.0 * (C * 100.0) * (C * 100.0) * (C * 100.0); + double pDotTop = _8_PI_2 * radius_6 * magField_2; + double pDotBottom = _3_C_3 * p_MomentOfInteria * period; + double pDot = pDotTop / pDotBottom; // period derivative + + return(-pDot * p_Omega / period); // convert period derivative to frequency derivative, which is what is recorded in the output } /* * Calculates and sets pulsar parameters at birth of pulsar * + * Modifies the following class member variables: + * + * m_AngularMomentum + * m_MomentOfInertia_CGS + * m_PulsarDetails.birthPeriod + * m_PulsarDetails.birthSpinDownRate + * m_PulsarDetails.magneticField + * m_PulsarDetails.spinDownRate + * m_PulsarDetails.spinFrequency + * + * * void CalculateAndSetPulsarParameters() */ void NS::CalculateAndSetPulsarParameters() { - m_PulsarDetails.magneticField = PPOW(10.0, CalculatePulsarBirthMagneticField_Static()) * GAUSS_TO_TESLA ; // magnetic field in Gauss -> convert to Tesla - m_PulsarDetails.spinPeriod = CalculatePulsarBirthSpinPeriod_Static(); // spin period in ms + m_PulsarDetails.magneticField = PPOW(10.0, CalculatePulsarBirthMagneticField()) * GAUSS_TO_TESLA; // magnetic field in Gauss -> convert to Tesla + m_PulsarDetails.spinPeriod = CalculatePulsarBirthSpinPeriod(); // spin period in ms m_PulsarDetails.spinFrequency = _2_PI / (m_PulsarDetails.spinPeriod * SECONDS_IN_MS); - m_PulsarDetails.birthPeriod = m_PulsarDetails.spinPeriod / 1000.0; //convert from ms to s + m_PulsarDetails.birthPeriod = m_PulsarDetails.spinPeriod / 1000.0; // convert from ms to s - m_MomentOfInertia = CalculateMomentOfInertia_Static(m_Mass, m_Radius * RSOL_TO_KM) ; // in CGS g cm^2 - // Note we convert neutronStarMomentOfInertia from CGS to SI here + m_MomentOfInertia_CGS = CalculateMomentOfInertiaCGS(); // in CGS g cm^2 + + // Note we convert neutronStarMomentOfInertia from CGS to SI here constexpr double factor = G_TO_KG * CM_TO_M * CM_TO_M; - m_PulsarDetails.spinDownRate = CalculateSpinDownRate_Static(m_PulsarDetails.spinFrequency, m_MomentOfInertia, m_PulsarDetails.magneticField, m_Radius * RSOL_TO_KM); + m_PulsarDetails.spinDownRate = CalculateSpinDownRate(m_PulsarDetails.spinFrequency, m_MomentOfInertia_CGS, m_PulsarDetails.magneticField, m_Radius * RSOL_TO_KM); m_PulsarDetails.birthSpinDownRate = m_PulsarDetails.spinDownRate; - m_AngularMomentum = _2_PI * m_MomentOfInertia / (m_PulsarDetails.spinPeriod * SECONDS_IN_MS) * factor; // in kg m^2 sec^-1 - + m_AngularMomentum = _2_PI * m_MomentOfInertia_CGS / (m_PulsarDetails.spinPeriod * SECONDS_IN_MS) * factor; // in kg m^2 sec^-1 } @@ -365,38 +379,42 @@ void NS::CalculateAndSetPulsarParameters() { */ void NS::SpinDownIsolatedPulsar(const double p_Stepsize) { - double NSRadius_IN_CM = m_Radius * RSOL_TO_KM * KM_TO_CM ; - double NSRadius_3 = NSRadius_IN_CM * NSRadius_IN_CM * NSRadius_IN_CM; - double NSRadius_6 = NSRadius_3 * NSRadius_3; - constexpr double _8_PI_2 = 8.0 * M_PI * M_PI; - constexpr double _3_C_3 = 3.0 * C * C * C * 1000000.0; + double NSradius_IN_CM = m_Radius * RSOL_TO_KM * KM_TO_CM; + double NSradius_3 = NSradius_IN_CM * NSradius_IN_CM * NSradius_IN_CM; + double NSradius_6 = NSradius_3 * NSradius_3; + constexpr double _8_PI_2 = 8.0 * PI_2; + constexpr double _3_C_3 = 3.0 * C * C * C * 1000000.0; - double initialMagField = m_PulsarDetails.magneticField; // (in T) - double initialMagField_G = initialMagField * TESLA_TO_GAUSS; - double initialSpinPeriod = 2.0 * M_PI / m_PulsarDetails.spinFrequency; - double magFieldLowerLimit = PPOW(10.0, OPTIONS->PulsarLog10MinimumMagneticField()) * GAUSS_TO_TESLA; - double magFieldLowerLimit_G = magFieldLowerLimit * TESLA_TO_GAUSS; - double momentOfInertia = m_MomentOfInertia; - double tau = OPTIONS->PulsarMagneticFieldDecayTimescale() * MYR_TO_YEAR * SECONDS_IN_YEAR; - - // calculate isolated decay of the magnetic field for a neutron star see Equation 6 in arXiv:0903.3538v2 - m_PulsarDetails.magneticField = magFieldLowerLimit + (initialMagField - magFieldLowerLimit) * exp(-p_Stepsize / tau); //Update pulsar magnetic field in SI. - // calculate the spin down rate for isolated neutron stars, see Equation 6 in arxiv:1912.02415. The rest of the calculations are carried out in cgs. - double constant_2 = (_8_PI_2 * NSRadius_6) / (_3_C_3 * momentOfInertia); - double term1 = magFieldLowerLimit_G * magFieldLowerLimit_G * p_Stepsize; - double term2 = tau * magFieldLowerLimit_G * ( m_PulsarDetails.magneticField * TESLA_TO_GAUSS - initialMagField_G); - double term3 = (tau / 2.0) * (TESLA_TO_GAUSS * TESLA_TO_GAUSS * (m_PulsarDetails.magneticField * m_PulsarDetails.magneticField) - (initialMagField_G * initialMagField_G)); - double PSquared = 2 * constant_2 * (term1 - term2 - term3) + (initialSpinPeriod * initialSpinPeriod); + double initialMagField = m_PulsarDetails.magneticField; // (in T) + double initialMagField_G = initialMagField * TESLA_TO_GAUSS; + double initialSpinPeriod = PI_2 / m_PulsarDetails.spinFrequency; + double magFieldLowerLimit = PPOW(10.0, OPTIONS->PulsarLog10MinimumMagneticField()) * GAUSS_TO_TESLA; + double magFieldLowerLimit_G = magFieldLowerLimit * TESLA_TO_GAUSS; + double tau = OPTIONS->PulsarMagneticFieldDecayTimescale() * MYR_TO_YEAR * SECONDS_IN_YEAR; + + // calculate isolated decay of the magnetic field for a neutron star + // see Equation 6 in arXiv:0903.3538v2 + m_PulsarDetails.magneticField = magFieldLowerLimit + (initialMagField - magFieldLowerLimit) * exp(-p_Stepsize / tau); // update pulsar magnetic field in SI. - double P_f = std::sqrt(PSquared); - m_PulsarDetails.spinFrequency = 2.0 * M_PI / P_f; // pulsar spin frequency + // calculate the spin down rate for isolated neutron stars + // see Equation 6 in arxiv:1912.02415 + // The rest of the calculations are carried out in cgs. + double constant2 = (_8_PI_2 * NSradius_6) / (_3_C_3 * m_MomentOfInertia_CGS); + double term1 = magFieldLowerLimit_G * magFieldLowerLimit_G * p_Stepsize; + double term2 = tau * magFieldLowerLimit_G * ( m_PulsarDetails.magneticField * TESLA_TO_GAUSS - initialMagField_G); + double term3 = (tau / 2.0) * (TESLA_TO_GAUSS * TESLA_TO_GAUSS * (m_PulsarDetails.magneticField * m_PulsarDetails.magneticField) - (initialMagField_G * initialMagField_G)); + double Psquared = 2 * constant2 * (term1 - term2 - term3) + (initialSpinPeriod * initialSpinPeriod); + + double P_f = std::sqrt(Psquared); + m_PulsarDetails.spinFrequency = _2_PI / P_f; // pulsar spin frequency - // calculate the spin down rate for isolated neutron stars, see Equation 4 in arXiv:0903.3538v2 (Our version is in cgs) - double pDotTop = constant_2 * TESLA_TO_GAUSS * TESLA_TO_GAUSS * m_PulsarDetails.magneticField * m_PulsarDetails.magneticField; - double pDot = pDotTop / P_f; - m_PulsarDetails.spinDownRate = -2.0 * M_PI * pDot / (P_f * P_f); + // calculate the spin down rate for isolated neutron stars + // see Equation 4 in arXiv:0903.3538v2 (Our version is in cgs) + double pDotTop = constant2 * TESLA_TO_GAUSS * TESLA_TO_GAUSS * m_PulsarDetails.magneticField * m_PulsarDetails.magneticField; + double pDot = pDotTop / P_f; + m_PulsarDetails.spinDownRate = -_2_PI * pDot / (P_f * P_f); - m_AngularMomentum = m_PulsarDetails.spinFrequency * momentOfInertia; // angular momentum of star + m_AngularMomentum = m_PulsarDetails.spinFrequency * m_MomentOfInertia_CGS; // angular momentum of star } @@ -426,51 +444,55 @@ void NS::SpinDownIsolatedPulsar(const double p_Stepsize) { */ void NS::UpdateMagneticFieldAndSpin(const bool p_CommonEnvelope, const bool p_RecycledNS, const double p_Stepsize, const double p_MassGainPerTimeStep, const double p_Epsilon) { - constexpr double unitsMoI = G_TO_KG * CM_TO_M * CM_TO_M; - double initialMagField = m_PulsarDetails.magneticField; // (in T) + constexpr double unitsMoI = G_TO_KG * CM_TO_M * CM_TO_M; // converts CGS -> SI + + double initialMagField = m_PulsarDetails.magneticField; // (in T) double magFieldLowerLimit = PPOW(10.0, OPTIONS->PulsarLog10MinimumMagneticField()) * GAUSS_TO_TESLA; double kappa = OPTIONS->PulsarMagneticFieldDecayMassscale() * MSOL_TO_KG; if ((!p_RecycledNS && !p_CommonEnvelope) || (!p_RecycledNS && utils::Compare(p_MassGainPerTimeStep, 0.0) == 0 )) { - //These are the ''classical'' isolated pulsars + // these are the ''classical'' isolated pulsars SpinDownIsolatedPulsar(p_Stepsize); } - else if ((m_PulsarDetails.spinFrequency < 2.0 * M_PI * 1000.0) && (p_RecycledNS || p_CommonEnvelope) && utils::Compare(p_MassGainPerTimeStep, 0.0) > 0) { - //This part of the code does pulsar recycling through acretion - //Recycling happens for pulsar with spin period larger than 1 ms and in a binary system with mass transfer - //The pulsar being recycled is either in a common envolope, or should have started the recycling process in previous time steps. - double mass_kg = m_Mass * MSOL_TO_KG; //in kg - double r_m = m_Radius * RSOL_TO_KM * 1000.0; //in meters + else if (utils::Compare(m_PulsarDetails.spinFrequency, _2_PI * 1000.0) < 0 && + (p_RecycledNS || p_CommonEnvelope) && utils::Compare(p_MassGainPerTimeStep, 0.0) > 0) { - double MOI_SI = m_MomentOfInertia * unitsMoI; - double angularMomentum_SI = m_AngularMomentum * unitsMoI; - - double newPulsarMagneticField = (initialMagField - magFieldLowerLimit) * exp(-1 * p_MassGainPerTimeStep / 1000.0 / kappa) + magFieldLowerLimit ; + // his part of the code does pulsar recycling through acretion + // recycling happens for pulsar with spin period larger than 1 ms and in a binary system with mass transfer + // the pulsar being recycled is either in a common envolope, or should have started the recycling process in previous time steps. + double mass_kg = m_Mass * MSOL_TO_KG; // in kg + double r_m = m_Radius * RSOL_TO_KM * 1000.0; // in metres - - //Calculate the Alfven radius for an accreting neutron star, see Equation 8 in arXiv:0903.3538v2 - double mDot = p_MassGainPerTimeStep / 1000.0 / p_Stepsize ; - double R_M_6 = r_m * r_m * r_m * r_m * r_m * r_m; - double B_4 = newPulsarMagneticField * newPulsarMagneticField * newPulsarMagneticField * newPulsarMagneticField; - double R_a_top = 8.0 * R_M_6 * R_M_6 * B_4; - double R_a_bot = mass_kg * mDot * mDot * G; - double alfvenRadius = PPOW(R_a_top / R_a_bot, 1.0/7.0); + double MoI_SI = m_MomentOfInertia_CGS * unitsMoI; + double angularMomentum_SI = m_AngularMomentum * unitsMoI; + + double newPulsarMagneticField = (initialMagField - magFieldLowerLimit) * exp(-p_MassGainPerTimeStep / 1000.0 / kappa) + magFieldLowerLimit; - // calculate the difference in the keplerian angular velocity and surface angular velocity of the neutron star in m - see Equation 2 in 1994MNRAS.269..455J - double keplerianVelocityAtAlfvenRadius = std::sqrt(2.0 * G * mass_kg / alfvenRadius); - double keplerianAngularVelocityAtAlfvenRadius = 4.0 * M_PI * keplerianVelocityAtAlfvenRadius / alfvenRadius; - double velocityDifference = keplerianAngularVelocityAtAlfvenRadius - m_PulsarDetails.spinFrequency; - - // calculate the change in angular momentum due to accretion, see Equation 12 in arXiv:0805.0059/ Equation 8 in arxiv:1912.02415 - double Jdot = p_Epsilon * velocityDifference * alfvenRadius * alfvenRadius * mDot ; + // calculate the Alfven radius for an accreting neutron star + // see Equation 8 in arXiv:0903.3538v2 + double mDot = p_MassGainPerTimeStep / 1000.0 / p_Stepsize; + double R_M_6 = r_m * r_m * r_m * r_m * r_m * r_m; + double B_4 = newPulsarMagneticField * newPulsarMagneticField * newPulsarMagneticField * newPulsarMagneticField; + double R_a_top = 8.0 * R_M_6 * R_M_6 * B_4; + double R_a_bot = mass_kg * mDot * mDot * G; + double alfvenRadius = PPOW(R_a_top / R_a_bot, 1.0 / 7.0); - angularMomentum_SI = angularMomentum_SI + Jdot * p_Stepsize ; + // calculate the difference in the keplerian angular velocity and surface angular velocity of the neutron star in m + // see Equation 2 in 1994MNRAS.269..455J + double keplerianVelocityAtAlfvenRadius = std::sqrt(2.0 * G * mass_kg / alfvenRadius); + double keplerianAngularVelocityAtAlfvenRadius = 4.0 * M_PI * keplerianVelocityAtAlfvenRadius / alfvenRadius; + double velocityDifference = keplerianAngularVelocityAtAlfvenRadius - m_PulsarDetails.spinFrequency; + + // calculate the change in angular momentum due to accretion + // see Equation 12 in arXiv:0805.0059/ Equation 8 in arxiv:1912.02415 + double Jdot = p_Epsilon * velocityDifference * alfvenRadius * alfvenRadius * mDot; + angularMomentum_SI = angularMomentum_SI + Jdot * p_Stepsize; - if (angularMomentum_SI / MOI_SI > 0) { - m_PulsarDetails.magneticField = newPulsarMagneticField ; - m_PulsarDetails.spinFrequency = angularMomentum_SI / MOI_SI; - m_PulsarDetails.spinDownRate = Jdot / MOI_SI; - m_AngularMomentum = angularMomentum_SI / unitsMoI; + if (utils::Compare(angularMomentum_SI / MoI_SI, 0.0) > 0) { + m_PulsarDetails.magneticField = newPulsarMagneticField; + m_PulsarDetails.spinFrequency = angularMomentum_SI / MoI_SI; + m_PulsarDetails.spinDownRate = Jdot / MoI_SI; + m_AngularMomentum = angularMomentum_SI / unitsMoI; } else { SpinDownIsolatedPulsar(p_Stepsize); diff --git a/src/NS.h b/src/NS.h index fbcdb18e5..ca46ca4f6 100755 --- a/src/NS.h +++ b/src/NS.h @@ -33,7 +33,7 @@ class NS: virtual public BaseStar, public Remnants { static double CalculateLuminosityOnPhase_Static(const double p_Mass, const double p_Time); - static double CalculatePulsarBirthSpinPeriod_Static(); + static double CalculatePulsarBirthSpinPeriod(); static double CalculateRadiusOnPhaseInKM_Static(const double p_Mass); // Radius on phase in km static double CalculateRadiusOnPhase_Static(const double p_Mass) { return CalculateRadiusOnPhaseInKM_Static(p_Mass) * KM_TO_RSOL; } // Radius on phase in Rsol @@ -42,56 +42,55 @@ class NS: virtual public BaseStar, public Remnants { protected: - + void Initialise() { m_StellarType = STELLAR_TYPE::NEUTRON_STAR; // Set stellar type CalculateTimescales(); // Initialise timescales //Set internal properties to zero to avoid meaningless values - m_Age = 0.0; - m_COCoreMass = 0.0; - m_HeCoreMass = 0.0; - m_CoreMass = 0.0; - m_Mass0 = 0.0; + m_Age = 0.0; + m_COCoreMass = 0.0; + m_HeCoreMass = 0.0; + m_CoreMass = 0.0; + m_Mass0 = 0.0; - m_Radius = NS::CalculateRadiusOnPhase_Static(m_Mass); // Set the NS radius, in Rsol - m_Luminosity = NS::CalculateLuminosityOnPhase_Static(m_Mass, m_Age); // Set the NS luminosity + m_Radius = CalculateRadiusOnPhase_Static(m_Mass); // Set the NS radius, in Rsol + m_Luminosity = CalculateLuminosityOnPhase_Static(m_Mass, m_Age); // Set the NS luminosity CalculateAndSetPulsarParameters(); } - double m_MomentOfInertia; // in CGS g cm^2 - double m_AngularMomentum; // Current angular momentum in (Msol AU^2 yr-1) - + double m_MomentOfInertia_CGS; // MoI in CGS - only required in NS class + + // member functions - alphabetically - void CalculateAndSetPulsarParameters(); + void CalculateAndSetPulsarParameters(); - double CalculateLuminosityOnPhase() const { return CalculateLuminosityOnPhase_Static(m_Mass, m_Age); } // Use class member variables + double CalculateLuminosityOnPhase() const { return CalculateLuminosityOnPhase_Static(m_Mass, m_Age); } // Use class member variables - double CalculateMassLossRate() { return 0.0; } // Ensure that NSs don't lose mass in winds + double CalculateMassLossRate() { return 0.0; } // Ensure that NSs don't lose mass in winds - static double CalculateMomentOfInertia_Static(const double p_Mass, const double p_Radius); + double CalculateMomentOfInertiaCGS() const; // MoI in CGS + double CalculateMomentOfInertia(const double p_RemnantRadius = 0.0) const { return CalculateMomentOfInertiaCGS() / MSOL_TO_G / RSOL_TO_CM / RSOL_TO_CM; } // MoI (default is solar units) + double CalculateMomentOfInertiaAU(const double p_RemnantRadius = 0.0) const { return CalculateMomentOfInertia() * RSOL_TO_AU * RSOL_TO_AU; } - static double CalculatePulsarBirthMagneticField_Static(); + double CalculatePulsarBirthMagneticField(); - double CalculateRadiusOnPhase() const { return CalculateRadiusOnPhase_Static(m_Mass); } // Use class member variables - returns radius in Rsol + double CalculateRadiusOnPhase() const { return CalculateRadiusOnPhase_Static(m_Mass); } // Use class member variables - returns radius in Rsol - static double CalculateSpinDownRate_Static(const double p_Omega, - const double p_MomentOfInteria, - const double p_MagField, - const double p_Radius); + double CalculateSpinDownRate(const double p_Omega, const double p_MomentOfInteria, const double p_MagField, const double p_Radius) const; - double ChooseTimestep(const double p_Time) const; + double ChooseTimestep(const double p_Time) const; - STELLAR_TYPE EvolveToNextPhase() { return STELLAR_TYPE::BLACK_HOLE; } + STELLAR_TYPE EvolveToNextPhase() { return STELLAR_TYPE::BLACK_HOLE; } - bool ShouldEvolveOnPhase() const { return (m_Mass <= OPTIONS->MaximumNeutronStarMass()); } // Evolve as a neutron star unless mass > maximum neutron star mass (e.g. through accretion) - void SpinDownIsolatedPulsar(const double p_Stepsize); - void UpdateMagneticFieldAndSpin(const bool p_CommonEnvelope, - const bool p_RecycledNS, - const double p_Stepsize, - const double p_MassGainPerTimeStep, - const double p_Epsilon); + bool ShouldEvolveOnPhase() const { return (m_Mass <= OPTIONS->MaximumNeutronStarMass()); } // Evolve as a neutron star unless mass > maximum neutron star mass (e.g. through accretion) + void SpinDownIsolatedPulsar(const double p_Stepsize); + void UpdateMagneticFieldAndSpin(const bool p_CommonEnvelope, + const bool p_RecycledNS, + const double p_Stepsize, + const double p_MassGainPerTimeStep, + const double p_Epsilon); }; diff --git a/src/Options.cpp b/src/Options.cpp index 5ddabfad5..c32f27991 100644 --- a/src/Options.cpp +++ b/src/Options.cpp @@ -474,6 +474,10 @@ void Options::OptionValues::Initialise() { m_CommonEnvelopeRecombinationEnergyDensity = 1.5E13; + // Tides + m_EnableTides = false; // default is no tides + + // Zetas m_StellarZetaPrescription.type = ZETA_PRESCRIPTION::SOBERMAN; m_StellarZetaPrescription.typeString = ZETA_PRESCRIPTION_LABEL.at(m_StellarZetaPrescription.type); @@ -751,6 +755,12 @@ bool Options::AddOptions(OptionValues *p_Options, po::options_description *p_Opt ("Print detailed output to file (default = " + std::string(p_Options->m_DetailedOutput ? "TRUE" : "FALSE") + ")").c_str() ) + ( + "enable-tides", + po::value(&p_Options->m_EnableTides)->default_value(p_Options->m_EnableTides)->implicit_value(true), + ("Enable tides (default = " + std::string(p_Options->m_EnableTides ? "TRUE" : "FALSE") + ")").c_str() + ) + ( "enable-warnings", po::value(&p_Options->m_EnableWarnings)->default_value(p_Options->m_EnableWarnings)->implicit_value(true), diff --git a/src/Options.h b/src/Options.h index d3157df48..37e2b84a7 100755 --- a/src/Options.h +++ b/src/Options.h @@ -333,6 +333,7 @@ class Options { "eccentricity-distribution", "eccentricity-max", "eccentricity-min", + "enable-tides", "evolve-double-white-dwarfs", "evolve-pulsars", "evolve-unbound-systems", @@ -437,6 +438,7 @@ class Options { "detailed-output", "eccentricity-distribution", + "enable-tides", "enable-warnings", "envelope-state-prescription", "errors-to-file", @@ -896,6 +898,10 @@ class Options { double m_CommonEnvelopeRecombinationEnergyDensity; // Factor using to calculate the binding energy depending on the mass of the envelope. (default = 1.5x10^13 erg/g) + // Tides + bool m_EnableTides; // Whether to enable tides (default = False) + + // Zetas ENUM_OPT m_StellarZetaPrescription; // Prescription to use for calculating stellar zetas (default = SOBERMAN) @@ -1206,6 +1212,7 @@ class Options { bool DebugToFile() const { return m_CmdLine.optionValues.m_DebugToFile; } bool DetailedOutput() const { return m_CmdLine.optionValues.m_DetailedOutput; } + bool EnableTides() const { return OPT_VALUE("enable-tides", m_EnableTides, true); } bool EnableWarnings() const { return m_CmdLine.optionValues.m_EnableWarnings; } bool ErrorsToFile() const { return m_CmdLine.optionValues.m_ErrorsToFile; } double Eccentricity() const { return OPT_VALUE("eccentricity", m_Eccentricity, true); } diff --git a/src/changelog.h b/src/changelog.h index 0a77a1b38..7c086b25d 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1057,8 +1057,10 @@ // - Fixed a few typos, a little code cleanup. // 02.39.01 LC - Sep 01, 2023 - Defect repair: // - Fix for issue #945 - made HeSD SN types a sub-class of SNIA types. +// 02.40.00 JR - Oct 25, 2023 - Enhancement, a little cleanup: +// - Added naive tides implementation. Functionality enabled with new option `--enable-tides`. Default is no tides. -const std::string VERSION_STRING = "02.39.01"; +const std::string VERSION_STRING = "02.40.00"; # endif // __changelog_h__ diff --git a/src/constants.h b/src/constants.h index 3b80565eb..9d6e79599 100755 --- a/src/constants.h +++ b/src/constants.h @@ -177,10 +177,10 @@ extern OBJECT_ID globalObjectId; // // I've added _2_PI and SQRT_M_2_PI below -#undef COMPARE_WITH_TOLERANCE // define/undef this to compare floats with/without tolerance (see FLOAT_TOLERANCE_ABSOLUTE, FLOAT_TOLERANCE_RELATIVE and Compare() function) +#undef COMPARE_GLOBAL_TOLERANCE // define/undef this to compare floats with/without tolerance (see FLOAT_TOLERANCE_ABSOLUTE, FLOAT_TOLERANCE_RELATIVE and Compare() function) -constexpr double FLOAT_TOLERANCE_ABSOLUTE = 0.0000005; // Absolute tolerance for floating-point comparisons if COMPARE_WITH_TOLERANCE is defined -constexpr double FLOAT_TOLERANCE_RELATIVE = 0.0000005; // Relative tolerance for floating-point comparisons if COMPARE_WITH_TOLERANCE is defined +constexpr double FLOAT_TOLERANCE_ABSOLUTE = 0.0000005; // Absolute tolerance for floating-point comparisons if COMPARE_GLOBAL_TOLERANCE is defined +constexpr double FLOAT_TOLERANCE_RELATIVE = 0.0000005; // Relative tolerance for floating-point comparisons if COMPARE_GLOBAL_TOLERANCE is defined // initialisation constants @@ -232,6 +232,7 @@ constexpr double GAUSS_TO_TESLA = 1.0 / TESLA_TO_GAUSS; // constants constexpr double _2_PI = M_PI * 2; // 2PI +constexpr double PI_2 = M_PI * M_PI; // PI squared constexpr double SQRT_M_2_PI = 0.79788456080286536; // sqrt(2/PI) constexpr double DEGREE = M_PI / 180.0; // 1 degree in radians @@ -243,7 +244,7 @@ constexpr double HUBBLE_TIME = 1 / H0SI; constexpr double G = 6.67E-11; // Gravitational constant in m^3 kg^-1 s^-2 (more accurately known as G M_sol) constexpr double G_CGS = 6.6743E-8; // Gravitational constant in cm^3 g^-1 s^-2 -constexpr double G1 = 4.0 * M_PI * M_PI; // Gravitational constant in AU^3 Msol^-1 yr^-2 +constexpr double G1 = 4.0 * PI_2; // Gravitational constant in AU^3 Msol^-1 yr^-2 constexpr double G_SN = G * 1.0E-9 / KG_TO_MSOL; // Gravitational constant in km^3 Msol^-1 s^-2, for use in the ResolveSupernova() function constexpr double G_SOLAR_YEAR = 3.14E7; // Gravitational constant in Lsol Rsol yr Msol^-2 for calculating photon tiring limit @@ -316,16 +317,19 @@ constexpr double EPSILON_PULSAR = 1.0; constexpr double MIN_HMXRB_STAR_TO_ROCHE_LOBE_RADIUS_RATIO = 0.8; // Minimum value of stellar radius | Roche Lobe radius for visible HMXRBs -constexpr double ADAPTIVE_RLOF_FRACTION_DONOR_GUESS = 0.001; // Fraction of donor mass to use as guess in MassLossToFitInsideRocheLobe() -constexpr int ADAPTIVE_RLOF_MAX_ITERATIONS = 50; // Maximum number of iterations in MassLossToFitInsideRocheLobe() -constexpr double ADAPTIVE_RLOF_SEARCH_FACTOR = 2.0; // Search factor in MassLossToFitInsideRocheLobe() -constexpr int ADAPTIVE_MASS0_MAX_ITERATIONS = 50; // Maximum number of iterations in Mass0ToMatchDesiredCoreMass() -constexpr double ADAPTIVE_MASS0_SEARCH_FACTOR = 2.0; // Search factor in Mass0ToMatchDesiredCoreMass() +constexpr double ADAPTIVE_RLOF_FRACTION_DONOR_GUESS = 0.001; // Fraction of donor mass to use as guess in BaseBinaryStar::MassLossToFitInsideRocheLobe() +constexpr int ADAPTIVE_RLOF_MAX_ITERATIONS = 50; // Maximum number of iterations in BaseBinaryStar::MassLossToFitInsideRocheLobe() +constexpr double ADAPTIVE_RLOF_SEARCH_FACTOR = 2.0; // Search factor in BaseBinaryStar::MassLossToFitInsideRocheLobe() +constexpr int ADAPTIVE_MASS0_MAX_ITERATIONS = 50; // Maximum number of iterations in HG::Mass0ToMatchDesiredCoreMass() +constexpr double ADAPTIVE_MASS0_SEARCH_FACTOR = 2.0; // Search factor in HG::Mass0ToMatchDesiredCoreMass() constexpr double FARMER_PPISN_UPP_LIM_LIN_REGIME = 38.0; // Maximum CO core mass to result in the linear remnant mass regime of the FARMER PPISN prescription constexpr double FARMER_PPISN_UPP_LIM_QUAD_REGIME = 60.0; // Maximum CO core mass to result in the quadratic remnant mass regime of the FARMER PPISN prescription constexpr double FARMER_PPISN_UPP_LIM_INSTABILLITY = 140.0; // Maximum CO core mass to result in PI (upper edge of PISN gap) from FARMER PPISN prescription constexpr double STARTRACK_PPISN_HE_CORE_MASS = 45.0; // Helium core mass remaining following PPISN as assumed in StarTrack (Belczynski et al. 2017 https://arxiv.org/abs/1607.03116) +constexpr int TIDES_OMEGA_MAX_ITERATIONS = 1000; // Maximum number of iterations in BaseBinaryStar::OmegaAfterCircularisation() +constexpr double TIDES_OMEGA_SEARCH_FACTOR = 1.1; // Search factor in BaseBinaryStar::OmegaAfterCircularisation() + // logging constants @@ -590,11 +594,13 @@ enum class ERROR: int { RADIUS_NOT_POSITIVE, // radius is <= 0.0 - invalid RADIUS_NOT_POSITIVE_ONCE, // radius is <= 0.0 - invalid RESOLVE_SUPERNOVA_IMPROPERLY_CALLED, // ResolveSupernova() called, but m_Supernova->IsSNevent() is false + ROOT_FINDER_FAILED, // root finder threw an exception STELLAR_EVOLUTION_STOPPED, // evolution of current star stopped STELLAR_SIMULATION_STOPPED, // stellar simulation stopped SUGGEST_HELP, // suggest using --help TIMESTEP_BELOW_MINIMUM, // timestep too small - below minimum TOO_MANY_MASS0_ITERATIONS, // too many iterations in MASS0 root finder + TOO_MANY_OMEGA_ITERATIONS, // too many iterations in OMEGA root finder TOO_MANY_RLOF_ITERATIONS, // too many iterations in RLOF root finder UNEXPECTED_END_OF_FILE, // unexpected end of file UNEXPECTED_LOG_FILE_TYPE, // unexpected log file type @@ -730,11 +736,13 @@ const COMPASUnorderedMap> ERROR_CATA { ERROR::RADIUS_NOT_POSITIVE, { ERROR_SCOPE::ALWAYS, "Radius <= 0.0" }}, { ERROR::RADIUS_NOT_POSITIVE_ONCE, { ERROR_SCOPE::FIRST_IN_FUNCTION, "Radius <= 0.0" }}, { ERROR::RESOLVE_SUPERNOVA_IMPROPERLY_CALLED, { ERROR_SCOPE::ALWAYS, "ResolveSupernova() called, but m_Supernova->IsSNevent() is false" }}, + { ERROR::ROOT_FINDER_FAILED, { ERROR_SCOPE::ALWAYS, "Exception encountered in root finder" }}, { ERROR::STELLAR_EVOLUTION_STOPPED, { ERROR_SCOPE::ALWAYS, "Evolution of current star stopped" }}, { ERROR::STELLAR_SIMULATION_STOPPED, { ERROR_SCOPE::ALWAYS, "Stellar simulation stopped" }}, { ERROR::SUGGEST_HELP, { ERROR_SCOPE::ALWAYS, "Use option '-h' (or '--help') to see (descriptions of) available options" }}, { ERROR::TIMESTEP_BELOW_MINIMUM, { ERROR_SCOPE::ALWAYS, "Timestep below minimum - timestep taken" }}, { ERROR::TOO_MANY_MASS0_ITERATIONS, { ERROR_SCOPE::ALWAYS, "Reached maximum number of iterations when looking for effective initial mass Mass_0 to match desired stellar core of HG star following case A mass transfer" }}, + { ERROR::TOO_MANY_OMEGA_ITERATIONS, { ERROR_SCOPE::ALWAYS, "Reached maximum number of iterations when looking for omega when circularising and synchronising for tides" }}, { ERROR::TOO_MANY_RLOF_ITERATIONS, { ERROR_SCOPE::ALWAYS, "Reached maximum number of iterations when fitting star inside Roche Lobe in RLOF" }}, { ERROR::UNEXPECTED_END_OF_FILE, { ERROR_SCOPE::ALWAYS, "Unexpected end of file" }}, { ERROR::UNEXPECTED_LOG_FILE_TYPE, { ERROR_SCOPE::ALWAYS, "Unexpected log file type" }}, @@ -1795,6 +1803,7 @@ const COMPASUnorderedMap PROPERTY_TYPE_LABEL = { MDOT, \ MEAN_ANOMALY, \ METALLICITY, \ + MOMENT_OF_INERTIA, \ MZAMS, \ NUCLEAR_TIMESCALE, \ NUCLEAR_TIMESCALE_POST_COMMON_ENVELOPE, \ @@ -1949,6 +1958,7 @@ const COMPASUnorderedMap STAR_PROPERTY_LABEL = { { STAR_PROPERTY::MDOT, "MDOT" }, { STAR_PROPERTY::MEAN_ANOMALY, "MEAN_ANOMALY" }, { STAR_PROPERTY::METALLICITY, "METALLICITY" }, + { STAR_PROPERTY::MOMENT_OF_INERTIA, "MOMENT_OF_INERTIA"}, { STAR_PROPERTY::MZAMS, "MZAMS" }, { STAR_PROPERTY::NUCLEAR_TIMESCALE, "NUCLEAR_TIMESCALE" }, { STAR_PROPERTY::NUCLEAR_TIMESCALE_POST_COMMON_ENVELOPE, "NUCLEAR_TIMESCALE_POST_COMMON_ENVELOPE" }, @@ -2761,7 +2771,7 @@ typedef std::tuple PROPERTY_DETAIL // STAR_PROPERTIES_LABEL const std::map ANY_STAR_PROPERTY_DETAIL = { { ANY_STAR_PROPERTY::AGE, { TYPENAME::DOUBLE, "Age", "Myr", 16, 8 }}, - { ANY_STAR_PROPERTY::ANGULAR_MOMENTUM, { TYPENAME::DOUBLE, "Ang_Momentum", "Msol*AU^2*yr^-1", 14, 6 }}, + { ANY_STAR_PROPERTY::ANGULAR_MOMENTUM, { TYPENAME::DOUBLE, "Ang_Momentum", "Msol AU^2 yr^-1", 14, 6 }}, { ANY_STAR_PROPERTY::BINDING_ENERGY_AT_COMMON_ENVELOPE, { TYPENAME::DOUBLE, "Binding_Energy@CE", "erg", 14, 6 }}, { ANY_STAR_PROPERTY::BINDING_ENERGY_FIXED, { TYPENAME::DOUBLE, "BE_Fixed", "erg", 14, 6 }}, { ANY_STAR_PROPERTY::BINDING_ENERGY_NANJING, { TYPENAME::DOUBLE, "BE_Nanjing", "erg", 14, 6 }}, @@ -2835,6 +2845,7 @@ const std::map ANY_STAR_PROPERTY_DETAIL = { { ANY_STAR_PROPERTY::MASS_TRANSFER_DONOR_HISTORY, { TYPENAME::STRING, "MT_Donor_Hist", "-", 16, 1 }}, { ANY_STAR_PROPERTY::MDOT, { TYPENAME::DOUBLE, "Mdot", "Msol yr^-1", 14, 6 }}, { ANY_STAR_PROPERTY::METALLICITY, { TYPENAME::DOUBLE, "Metallicity@ZAMS", "-", 14, 6 }}, + { ANY_STAR_PROPERTY::MOMENT_OF_INERTIA, { TYPENAME::DOUBLE, "Moment_Of_Inertia", "Msol Rsol^2", 14, 6 }}, { ANY_STAR_PROPERTY::MZAMS, { TYPENAME::DOUBLE, "Mass@ZAMS", "Msol", 14, 6 }}, { ANY_STAR_PROPERTY::NUCLEAR_TIMESCALE, { TYPENAME::DOUBLE, "Tau_Nuclear", "Myr", 16, 8 }}, { ANY_STAR_PROPERTY::NUCLEAR_TIMESCALE_POST_COMMON_ENVELOPE, { TYPENAME::DOUBLE, "Tau_Nuclear>CE", "Myr", 16, 8 }}, @@ -3006,8 +3017,8 @@ const std::map BINARY_PROPERTY_DETAIL = { { BINARY_PROPERTY::SYSTEMIC_SPEED, { TYPENAME::DOUBLE, "SystemicSpeed", "kms^-1", 14, 6 }}, { BINARY_PROPERTY::TIME, { TYPENAME::DOUBLE, "Time", "Myr", 16, 8 }}, { BINARY_PROPERTY::TIME_TO_COALESCENCE, { TYPENAME::DOUBLE, "Coalescence_Time", "Myr", 16, 8 }}, - { BINARY_PROPERTY::TOTAL_ANGULAR_MOMENTUM, { TYPENAME::DOUBLE, "Ang_Momentum_Total", "Msol*AU^2*yr^-1", 14, 6 }}, - { BINARY_PROPERTY::TOTAL_ENERGY, { TYPENAME::DOUBLE, "Energy_Total", "Msol*AU^2*yr^-2", 14, 6 }}, + { BINARY_PROPERTY::TOTAL_ANGULAR_MOMENTUM, { TYPENAME::DOUBLE, "Ang_Momentum_Total", "Msol AU^2 yr^-1", 14, 6 }}, + { BINARY_PROPERTY::TOTAL_ENERGY, { TYPENAME::DOUBLE, "Energy_Total", "Msol AU^2 yr^-2", 14, 6 }}, { BINARY_PROPERTY::UNBOUND, { TYPENAME::BOOL, "Unbound", "State", 0, 0 }}, { BINARY_PROPERTY::ZETA_LOBE, { TYPENAME::DOUBLE, "Zeta_Lobe", "-", 14, 6 }}, { BINARY_PROPERTY::ZETA_STAR, { TYPENAME::DOUBLE, "Zeta_Star", "-", 14, 6 }} diff --git a/src/utils.cpp b/src/utils.cpp index 05051061c..326ab4a47 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -168,26 +168,53 @@ namespace utils { /* * Compare floating-point numbers with tolerance * - * Absolute and relative tolerance can be different - see constants.h - * Set relative tolerance = 0.0 to always use absolute - * Set absolute tolerance = 0.0 to always use relative - * Set both to zero for no tolerance - or #undef COMPARE_WITH_TOLERANCE for performance - * + * For comparisons using the global tolerance values (FLOAT_TOLERANCE_ABSOLUTE, FLOAT_TOLERANCE_RELATIVE): + * - Absolute and relative tolerance can be different - see constants.h + * - Set relative tolerance = 0.0 to always use absolute + * - Set absolute tolerance = 0.0 to always use relative + * - Set both to zero for no tolerance - or #undef COMPARE_GLOBAL_TOLERANCE for performance + * + * If p_Tolerance is > 0.0 it will be used in preference to the global tolerance values + * If p_Tolerance is > 0.0, then p_Absolute determines if p_tolerance should be treated as an absolute + * torelace (p_Absolute = true), or a relative tolerance (p_Absolete = false). + * * * int Compare(const double p_X, const double p_Y) * * @param [IN] p_X Floating-point value to be compared * @param [IN] p_Y Floating-point value to be compared + * @param [IN] p_Tolerance Floating-point tolerance value - if > 0.0 supersedes global tolerance + * @param [IN] p_Absolute Boolean indicatin whether p_Tolerance should be treated as absolute tolerance (true) or relative tolerance (false) * @return Integer indicating result of comparison: * -1 indicates p_X is less than p_Y * 0 indicates equality * 1 indicates p_X is greater than p_Y */ - int Compare(const double p_X, const double p_Y) { - #ifdef COMPARE_WITH_TOLERANCE - return (std::abs(p_X - p_Y) <= std::max(FLOAT_TOLERANCE_ABSOLUTE, FLOAT_TOLERANCE_RELATIVE * std::max(std::abs(p_X), fabs(p_Y)))) ? 0 : (p_X < p_Y ? -1 : 1); + int Compare(const double p_X, const double p_Y, const double p_Tolerance, const bool p_Absolute) { + #ifdef COMPARE_GLOBAL_TOLERANCE + if (p_Tolerance > 0.0) { // use tolerance passed? + if (p_Absolute) { // yes - absolute tolerance? + return (std::abs(p_X - p_Y) <= p_Tolerance) ? 0 : (p_X < p_Y ? -1 : 1); // yes + } + else { // no - relative tolerance + return (std::abs(p_X - p_Y) <= p_Tolerance * std::max(std::abs(p_X), fabs(p_Y))) ? 0 : (p_X < p_Y ? -1 : 1); + } + } + else { // use global tolerance + return (std::abs(p_X - p_Y) <= std::max(FLOAT_TOLERANCE_ABSOLUTE, FLOAT_TOLERANCE_RELATIVE * std::max(std::abs(p_X), fabs(p_Y)))) ? 0 : (p_X < p_Y ? -1 : 1); + } #else - return (p_X == p_Y) ? 0 : (p_X < p_Y ? -1 : 1); + if (p_Tolerance > 0.0) { // use tolerance passed? + if (p_Absolute) { // yes - absolute tolerance? + return (std::abs(p_X - p_Y) <= p_Tolerance) ? 0 : (p_X < p_Y ? -1 : 1); // yes + } + else { // no - relative tolerance + return (std::abs(p_X - p_Y) <= p_Tolerance * std::max(std::abs(p_X), fabs(p_Y))) ? 0 : (p_X < p_Y ? -1 : 1); + } + } + else { // no tolerance + return (p_X == p_Y) ? 0 : (p_X < p_Y ? -1 : 1); + } #endif } @@ -205,10 +232,9 @@ namespace utils { */ double ConvertPeriodInDaysToSemiMajorAxisInAU(const double p_Mass1, const double p_Mass2, const double p_Period) { - double a_cubed_SI_top = G * ((p_Mass1 * MSOL_TO_KG) + (p_Mass2 * MSOL_TO_KG)) * p_Period * p_Period * SECONDS_IN_DAY * SECONDS_IN_DAY; - double a_cubed_SI_bottom = 4.0 * M_PI * M_PI; - double a_cubed_SI = a_cubed_SI_top / a_cubed_SI_bottom; - double a_SI = std::cbrt(a_cubed_SI); + double a_cubed_SI_top = G * ((p_Mass1 * MSOL_TO_KG) + (p_Mass2 * MSOL_TO_KG)) * p_Period * p_Period * SECONDS_IN_DAY * SECONDS_IN_DAY; + double a_cubed_SI = a_cubed_SI_top / G1; + double a_SI = std::cbrt(a_cubed_SI); return a_SI / AU; } @@ -883,7 +909,7 @@ namespace utils { // Sampling function taken from binpop.f in NBODY6 do { - eccentricity = 0.23 * std::sqrt(-2.0 * log(RAND->Random())) * cos(2.0 * M_PI * RAND->Random()) + 0.38; + eccentricity = 0.23 * std::sqrt(-2.0 * log(RAND->Random())) * cos(_2_PI * RAND->Random()) + 0.38; } while (eccentricity < p_Min || eccentricity > p_Max); // JR: don't use utils::Compare() here break; @@ -892,7 +918,7 @@ namespace utils { // Sampling function taken from binpop.f in NBODY6 do { - eccentricity = 0.15 * std::sqrt(-2.0 * log(RAND->Random())) * cos(2.0 * M_PI * RAND->Random()) + 0.3; + eccentricity = 0.15 * std::sqrt(-2.0 * log(RAND->Random())) * cos(_2_PI * RAND->Random()) + 0.3; } while (eccentricity < p_Min or eccentricity > p_Max); // JR: don't use utils::Compare() here break; @@ -1100,7 +1126,7 @@ namespace utils { case MASS_RATIO_DISTRIBUTION::DUQUENNOYMAYOR1991: // mass ratio distribution from Duquennoy & Mayor (1991) (http://adsabs.harvard.edu/abs/1991A%26A...248..485D) do { // JR: todo: catch non-convergence? - q = 0.42 * std::sqrt(-2.0 * log(RAND->Random())) * cos(2.0 * M_PI * RAND->Random()) + 0.23; + q = 0.42 * std::sqrt(-2.0 * log(RAND->Random())) * cos(_2_PI * RAND->Random()) + 0.23; } while (q < p_Min || q > p_Max); // JR: don't use utils::Compare() here break; @@ -1238,7 +1264,7 @@ namespace utils { // Make sure that the drawn semi-major axis is in the range specified by the user do { // JR: todo: catch for non-convergence? - double periodInDays = PPOW(10.0, 2.3 * std::sqrt(-2.0 * log(RAND->Random())) * cos(2.0 * M_PI * RAND->Random()) + 4.8); + double periodInDays = PPOW(10.0, 2.3 * std::sqrt(-2.0 * log(RAND->Random())) * cos(_2_PI * RAND->Random()) + 4.8); semiMajorAxis = utils::ConvertPeriodInDaysToSemiMajorAxisInAU(p_Mass1, p_Mass2, periodInDays); // convert period in days to semi-major axis in AU } while (semiMajorAxis < p_AdistMin || semiMajorAxis > p_AdistMax); // JR: don't use utils::Compare() here break; diff --git a/src/utils.h b/src/utils.h index 591b04b2b..459e78e26 100755 --- a/src/utils.h +++ b/src/utils.h @@ -21,7 +21,7 @@ namespace utils { std::string CentreJustify(const std::string p_Str, std::size_t p_Width); - int Compare(const double p_X, const double p_Y); + int Compare(const double p_X, const double p_Y, const double p_Tolerance = -1.0, const bool p_Absolute = true); double ConvertPeriodInDaysToSemiMajorAxisInAU(const double p_Mass1, const double p_Mass2, const double p_Period); diff --git a/src/yaml.h b/src/yaml.h index 3ec9e417d..1f12b9464 100644 --- a/src/yaml.h +++ b/src/yaml.h @@ -113,6 +113,9 @@ namespace yaml { " --common-envelope-lambda-nanjing-use-rejuvenated-mass", " --revised-energy-formalism-nandez-ivanova", "", + " ### TIDES", + " --enable-tides", + "", " ### SUPERNOVAE, KICKS AND REMNANTS", " --allow-non-stripped-ECSN", " --pair-instability-supernovae",