diff --git a/online-docs/pages/User guide/Program options/program-options-list-defaults.rst b/online-docs/pages/User guide/Program options/program-options-list-defaults.rst index eb210cc76..484f8d196 100644 --- a/online-docs/pages/User guide/Program options/program-options-list-defaults.rst +++ b/online-docs/pages/User guide/Program options/program-options-list-defaults.rst @@ -845,6 +845,11 @@ Default = FIXED_MASS Amount of mass lost in neutrinos during BH formation (either as fraction or in solar masses, depending on the value of ``--neutrino-mass-loss-bh-formation``). |br| Default = 0.1 +**--neutron-star-accretion-in-ce** |br| +Neutron star accretion behavior in common envelope. |br| +Options: { ZERO, DISK, SURFACE } |br| +Default = ZERO + **--neutron-star-equation-of-state** |br| Neutron star equation of state. |br| Options: { SSE, ARP3 } |br| diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 4b1165e64..dc30d8b52 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -388,11 +388,11 @@ void BaseBinaryStar::SetRemainingValues() { m_Flags.mergesInHubbleTime = false; m_Unbound = false; - m_SystemicVelocity = Vector3d(); - m_NormalizedOrbitalAngularMomentumVector = Vector3d(); - m_ThetaE = DEFAULT_INITIAL_DOUBLE_VALUE; - m_PhiE = DEFAULT_INITIAL_DOUBLE_VALUE; - m_PsiE = DEFAULT_INITIAL_DOUBLE_VALUE; + m_SystemicVelocity = Vector3d(); + m_NormalizedOrbitalAngularMomentumVector = 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; @@ -1363,7 +1363,7 @@ bool BaseBinaryStar::ResolveSupernova() { 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 + double orbitalAngularMomentum = orbitalAngularMomentumVector.mag; // km^2 s^-1 - PostSN specific orbital angular momentum m_NormalizedOrbitalAngularMomentumVector = orbitalAngularMomentumVector/orbitalAngularMomentum; // set unit vector here to make printing out the inclination vector easier Vector3d eccentricityVector = cross(relativeVelocityVector, orbitalAngularMomentumVector) / @@ -2166,8 +2166,8 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) { } } - // Check for recycled pulsars. Not considering CEE as a way of recycling NSs. - if (!m_CEDetails.CEEnow && m_Accretor->IsOneOf({ STELLAR_TYPE::NEUTRON_STAR })) { // accretor is a neutron star + // Check for recycled pulsars. + if ((!m_CEDetails.CEEnow || OPTIONS->NeutronStarAccretionInCE() != NS_ACCRETION_IN_CE::ZERO) && m_Accretor->IsOneOf({ STELLAR_TYPE::NEUTRON_STAR })) { m_Donor->SetRLOFOntoNS(); // donor donated mass to a neutron star m_Accretor->SetRecycledNS(); // accretor is (was) a recycled NS } diff --git a/src/Mode b/src/Mode new file mode 100644 index 000000000..fbbdd8f76 --- /dev/null +++ b/src/Mode @@ -0,0 +1,495 @@ +BaseBinaryStar.cpp: bool sampled = OPTIONS->OptionSpecified("initial-mass-1") == 0 || +BaseBinaryStar.cpp: OPTIONS->OptionSpecified("initial-mass-2") == 0 || +BaseBinaryStar.cpp: (OPTIONS->OptionSpecified("metallicity") == 0 && OPTIONS->MetallicityDistribution() != METALLICITY_DISTRIBUTION::ZSOLAR) || +BaseBinaryStar.cpp: (OPTIONS->OptionSpecified("semi-major-axis") == 0 && OPTIONS->OptionSpecified("orbital-period") == 0) || +BaseBinaryStar.cpp: (OPTIONS->OptionSpecified("eccentricity") == 0 && OPTIONS->EccentricityDistribution() != ECCENTRICITY_DISTRIBUTION::ZERO); +BaseBinaryStar.cpp: kickParameters1.magnitudeRandomSpecified = OPTIONS->OptionSpecified("kick-magnitude-random-1") == 1; +BaseBinaryStar.cpp: kickParameters1.magnitudeRandom = OPTIONS->KickMagnitudeRandom1(); +BaseBinaryStar.cpp: kickParameters1.magnitudeSpecified = OPTIONS->OptionSpecified("kick-magnitude-1") == 1; +BaseBinaryStar.cpp: kickParameters1.magnitude = OPTIONS->KickMagnitude1(); +BaseBinaryStar.cpp: kickParameters1.phiSpecified = OPTIONS->OptionSpecified("kick-phi-1") == 1; +BaseBinaryStar.cpp: kickParameters1.phi = OPTIONS->SN_Phi1(); +BaseBinaryStar.cpp: kickParameters1.thetaSpecified = OPTIONS->OptionSpecified("kick-theta-1") == 1; +BaseBinaryStar.cpp: kickParameters1.theta = OPTIONS->SN_Theta1(); +BaseBinaryStar.cpp: kickParameters1.meanAnomalySpecified = OPTIONS->OptionSpecified("kick-mean-anomaly-1") == 1; +BaseBinaryStar.cpp: kickParameters1.meanAnomaly = OPTIONS->SN_MeanAnomaly1(); +BaseBinaryStar.cpp: kickParameters2.magnitudeRandomSpecified = OPTIONS->OptionSpecified("kick-magnitude-random-2") == 1; +BaseBinaryStar.cpp: kickParameters2.magnitudeRandom = OPTIONS->KickMagnitudeRandom2(); +BaseBinaryStar.cpp: kickParameters2.magnitudeSpecified = OPTIONS->OptionSpecified("kick-magnitude-2") == 1; +BaseBinaryStar.cpp: kickParameters2.magnitude = OPTIONS->KickMagnitude2(); +BaseBinaryStar.cpp: kickParameters2.phiSpecified = OPTIONS->OptionSpecified("kick-phi-2") == 1; +BaseBinaryStar.cpp: kickParameters2.phi = OPTIONS->SN_Phi2(); +BaseBinaryStar.cpp: kickParameters2.thetaSpecified = OPTIONS->OptionSpecified("kick-theta-2") == 1; +BaseBinaryStar.cpp: kickParameters2.theta = OPTIONS->SN_Theta2(); +BaseBinaryStar.cpp: kickParameters2.meanAnomalySpecified = OPTIONS->OptionSpecified("kick-mean-anomaly-2") == 1; +BaseBinaryStar.cpp: kickParameters2.meanAnomaly = OPTIONS->SN_MeanAnomaly2(); +BaseBinaryStar.cpp: double mass1 = OPTIONS->OptionSpecified("initial-mass-1") == 1 // user specified primary mass? +BaseBinaryStar.cpp: ? OPTIONS->InitialMass1() // yes, use it +BaseBinaryStar.cpp: : utils::SampleInitialMass(OPTIONS->InitialMassFunction(), +BaseBinaryStar.cpp: OPTIONS->InitialMassFunctionMax(), +BaseBinaryStar.cpp: OPTIONS->InitialMassFunctionMin(), +BaseBinaryStar.cpp: OPTIONS->InitialMassFunctionPower()); // no - asmple it +BaseBinaryStar.cpp: if (OPTIONS->OptionSpecified("initial-mass-2") == 1) { // user specified secondary mass? +BaseBinaryStar.cpp: mass2 = OPTIONS->InitialMass2(); // yes, use it +BaseBinaryStar.cpp: double q = OPTIONS->OptionSpecified("mass-ratio") == 1 // user specified mass ratio? +BaseBinaryStar.cpp: ? OPTIONS->MassRatio() // yes, use it +BaseBinaryStar.cpp: : utils::SampleMassRatio(OPTIONS->MassRatioDistribution(), +BaseBinaryStar.cpp: OPTIONS->MassRatioDistributionMax(), +BaseBinaryStar.cpp: OPTIONS->MassRatioDistributionMin()); // no - sample it +BaseBinaryStar.cpp: double metallicity = OPTIONS->OptionSpecified("metallicity") == 1 // user specified metallicity? +BaseBinaryStar.cpp: ? OPTIONS->Metallicity() // yes, use it +BaseBinaryStar.cpp: : utils::SampleMetallicity(OPTIONS->MetallicityDistribution(), +BaseBinaryStar.cpp: OPTIONS->MetallicityDistributionMax(), +BaseBinaryStar.cpp: OPTIONS->MetallicityDistributionMin()); // no, sample it +BaseBinaryStar.cpp: if (OPTIONS->OptionSpecified("semi-major-axis") == 1) { // user specified semi-major axis? +BaseBinaryStar.cpp: m_SemiMajorAxis = OPTIONS->SemiMajorAxis(); // yes, use it +BaseBinaryStar.cpp: if (OPTIONS->OptionSpecified("orbital-period") == 1) { // user specified orbital period? +BaseBinaryStar.cpp: m_SemiMajorAxis = utils::ConvertPeriodInDaysToSemiMajorAxisInAU(mass1, mass2, OPTIONS->OrbitalPeriod()); // yes - calculate semi-major axis from period +BaseBinaryStar.cpp: if (OPTIONS->OptionSpecified("semi-major-axis-distribution") == 1 || // user specified semi-major axis distribution, or +BaseBinaryStar.cpp: OPTIONS->OptionSpecified("orbital-period-distribution" ) == 0) { // user did not specify oprbital period distribution +BaseBinaryStar.cpp: m_SemiMajorAxis = utils::SampleSemiMajorAxis(OPTIONS->SemiMajorAxisDistribution(), +BaseBinaryStar.cpp: OPTIONS->SemiMajorAxisDistributionMax(), +BaseBinaryStar.cpp: OPTIONS->SemiMajorAxisDistributionMin(), +BaseBinaryStar.cpp: OPTIONS->SemiMajorAxisDistributionPower(), +BaseBinaryStar.cpp: OPTIONS->OrbitalPeriodDistributionMax(), +BaseBinaryStar.cpp: OPTIONS->OrbitalPeriodDistributionMin(), +BaseBinaryStar.cpp: double orbitalPeriod = utils::SampleOrbitalPeriod(OPTIONS->OrbitalPeriodDistribution(), +BaseBinaryStar.cpp: OPTIONS->OrbitalPeriodDistributionMax(), +BaseBinaryStar.cpp: OPTIONS->OrbitalPeriodDistributionMin()); +BaseBinaryStar.cpp: m_Eccentricity = OPTIONS->OptionSpecified("eccentricity") == 1 // user specified eccentricity? +BaseBinaryStar.cpp: ? OPTIONS->Eccentricity() // yes, use it +BaseBinaryStar.cpp: : utils::SampleEccentricity(OPTIONS->EccentricityDistribution(), +BaseBinaryStar.cpp: OPTIONS->EccentricityDistributionMax(), +BaseBinaryStar.cpp: OPTIONS->EccentricityDistributionMin()); // no, sample it +BaseBinaryStar.cpp: m_Star1 = OPTIONS->OptionSpecified("rotational-frequency-1") == 1 // user specified primary rotational frequency? +BaseBinaryStar.cpp: ? new BinaryConstituentStar(m_RandomSeed, mass1, metallicity, kickParameters1, OPTIONS->RotationalFrequency1() * SECONDS_IN_YEAR) // yes - use it (convert from Hz to cycles per year - see BaseStar::CalculateZAMSAngularFrequency()) +BaseBinaryStar.cpp: m_Star2 = OPTIONS->OptionSpecified("rotational-frequency-2") == 1 // user specified secondary rotational frequency? +BaseBinaryStar.cpp: ? new BinaryConstituentStar(m_RandomSeed, mass2, metallicity, kickParameters2, OPTIONS->RotationalFrequency2() * SECONDS_IN_YEAR) // yes - use it (convert from Hz to cycles per year - see BaseStar::CalculateZAMSAngularFrequency()) +BaseBinaryStar.cpp: if (rlof && OPTIONS->AllowRLOFAtBirth()) { // over-contact binaries at birth allowed? +BaseBinaryStar.cpp: m_Star1 = OPTIONS->OptionSpecified("rotational-frequency-1") == 1 // user specified primary rotational frequency? +BaseBinaryStar.cpp: ? new BinaryConstituentStar(m_RandomSeed, mass1, metallicity, kickParameters1, OPTIONS->RotationalFrequency1() * SECONDS_IN_YEAR) // yes - use it (convert from Hz to cycles per year - see BaseStar::CalculateZAMSAngularFrequency()) +BaseBinaryStar.cpp: m_Star2 = OPTIONS->OptionSpecified("rotational-frequency-2") == 1 // user specified secondary rotational frequency? +BaseBinaryStar.cpp: ? new BinaryConstituentStar(m_RandomSeed, mass2, metallicity, kickParameters2, OPTIONS->RotationalFrequency2() * SECONDS_IN_YEAR) // yes - use it (convert from Hz to cycles per year - see BaseStar::CalculateZAMSAngularFrequency()) +BaseBinaryStar.cpp: secondarySmallerThanMinimumMass = utils::Compare(mass2, OPTIONS->MinimumMassSecondary()) < 0; +BaseBinaryStar.cpp: bool ok = !((!OPTIONS->AllowRLOFAtBirth() && rlof) || (!OPTIONS->AllowTouchingAtBirth() && merger) || secondarySmallerThanMinimumMass); +BaseBinaryStar.cpp: if (OPTIONS->PopulationDataPrinting()) { // user wants to see details of binary? +BaseBinaryStar.cpp: if (OPTIONS->CHEMode() != CHE_MODE::NONE) { // CHE enabled? +BaseBinaryStar.cpp: m_JLoss = OPTIONS->MassTransferJloss(); +BaseBinaryStar.cpp: m_FractionAccreted = OPTIONS->MassTransferFractionAccreted(); +BaseBinaryStar.cpp: std::tie(ok, value) = OPTIONS->OptionValue(p_Property); // get the value +BaseBinaryStar.cpp: if (!OPTIONS->RLOFPrinting()) return ok; // do not print if printing option off +BaseBinaryStar.cpp: if (OPTIONS->HMXRBinaries()) { +BaseBinaryStar.cpp: if (!OPTIONS->BeBinaries()) return true; // do not print if printing option off +BaseBinaryStar.cpp: if (!OPTIONS->RLOFPrinting()) return; // nothing to do +BaseBinaryStar.cpp: if (!OPTIONS->BeBinaries() || !IsBeBinary()) return; // nothing to do; +BaseBinaryStar.cpp: double alphaCE = OPTIONS->CommonEnvelopeAlpha(); // CE efficiency parameter +BaseBinaryStar.cpp: if (OPTIONS->AllowMainSequenceStarToSurviveCommonEnvelope()) { // allow main sequence stars to survive CEE? +BaseBinaryStar.cpp: if (OPTIONS->CommonEnvelopeFormalism() == CE_FORMALISM::ENERGY) { +BaseBinaryStar.cpp: else if ( OPTIONS->CommonEnvelopeFormalism() == CE_FORMALISM::TWO_STAGE ) { +BaseBinaryStar.cpp: if (!OPTIONS->AllowRadiativeEnvelopeStarToSurviveCommonEnvelope()) { // stellar merger +BaseBinaryStar.cpp: if (OPTIONS->CHEMode() != CHE_MODE::NONE) m_Star1->SetOmega(omega); +BaseBinaryStar.cpp: if (OPTIONS->CHEMode() != CHE_MODE::NONE) m_Star2->SetOmega(omega); +BaseBinaryStar.cpp: if (m_RLOFDetails.immediateRLOFPostCEE == true && !OPTIONS->AllowImmediateRLOFpostCEToSurviveCommonEnvelope()) { // is there immediate post-CE RLOF which is not allowed? +BaseBinaryStar.cpp: if (!(m_Star1->IsOneOf(MAIN_SEQUENCE) && m_Star2->IsOneOf(MAIN_SEQUENCE) && OPTIONS->EvolveMainSequenceMergers())) +BaseBinaryStar.cpp: switch (OPTIONS->MassTransferAngularMomentumLossPrescription()) { // which precription? +BaseBinaryStar.cpp: ? OPTIONS->MassTransferJlossMacLeodLinearFractionDegen() +BaseBinaryStar.cpp: : OPTIONS->MassTransferJlossMacLeodLinearFractionNonDegen(); +BaseBinaryStar.cpp: case MT_ANGULAR_MOMENTUM_LOSS_PRESCRIPTION::ARBITRARY : gamma = OPTIONS->MassTransferJloss(); break; +BaseBinaryStar.cpp: if (OPTIONS->UseMassTransfer() && m_MassTransfer) { +BaseBinaryStar.cpp: if (OPTIONS->UseMassLoss()) { // mass loss enabled? +BaseBinaryStar.cpp: if (!OPTIONS->UseMassTransfer()) return; // mass transfer not enabled - nothing to do +BaseBinaryStar.cpp: if (OPTIONS->CHEMode() != CHE_MODE::NONE && HasTwoOf({STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS}) && HasStarsTouching()) { // CHE enabled and both stars CH? +BaseBinaryStar.cpp: if (OPTIONS->MassTransferAngularMomentumLossPrescription() != MT_ANGULAR_MOMENTUM_LOSS_PRESCRIPTION::ARBITRARY) { // arbitrary angular momentum loss prescription? +BaseBinaryStar.cpp: bool caseBBAlwaysStable = OPTIONS->CaseBBStabilityPrescription() == CASE_BB_STABILITY_PRESCRIPTION::ALWAYS_STABLE; +BaseBinaryStar.cpp: bool caseBBAlwaysUnstable = OPTIONS->CaseBBStabilityPrescription() == CASE_BB_STABILITY_PRESCRIPTION::ALWAYS_UNSTABLE; +BaseBinaryStar.cpp: bool caseBBAlwaysUnstableOntoNSBH = OPTIONS->CaseBBStabilityPrescription() == CASE_BB_STABILITY_PRESCRIPTION::ALWAYS_STABLE_ONTO_NSBH; +BaseBinaryStar.cpp: else if (OPTIONS->QCritPrescription() != QCRIT_PRESCRIPTION::NONE) { // Determine stability based on critical mass ratios +BaseBinaryStar.cpp: if ((!m_CEDetails.CEEnow || OPTIONS->NeutronStarAccretionInCE() != NS_ACCRETION_IN_CE::ZERO) && m_Accretor->IsOneOf({ STELLAR_TYPE::NEUTRON_STAR })) { +BaseBinaryStar.cpp: if (OPTIONS->CHEMode() != CHE_MODE::NONE && HasTwoOf({STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS})) { // CHE enabled and both stars CH? +BaseBinaryStar.cpp: if (OPTIONS->CirculariseBinaryDuringMassTransfer()) { // circularise binary to the periapsis separation? +BaseBinaryStar.cpp: m_SemiMajorAxis *= OPTIONS->AngularMomentumConservationDuringCircularisation() // yes - conserve angular momentum? +BaseBinaryStar.cpp: if (OPTIONS->CHEMode() != CHE_MODE::NONE) m_Star1->SetOmega(omega); +BaseBinaryStar.cpp: if (OPTIONS->CHEMode() != CHE_MODE::NONE) m_Star2->SetOmega(omega); +BaseBinaryStar.cpp: !(OPTIONS->CHEMode() != CHE_MODE::NONE && HasTwoOf({STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS}))) { // yes - avoid CEE if CH+CH +BaseBinaryStar.cpp: if (!m_Unbound && OPTIONS->TidesPrescription() != TIDES_PRESCRIPTION::NONE) { +BaseBinaryStar.cpp: if (OPTIONS->TidesPrescription() == TIDES_PRESCRIPTION::KAPIL2024) { +BaseBinaryStar.cpp: if (OPTIONS->PopulationDataPrinting()) { +BaseBinaryStar.cpp: if (!OPTIONS->TimestepsFileName().empty()) { // have timesteps filename? +BaseBinaryStar.cpp: std::tie(error, timesteps) = utils::ReadTimesteps(OPTIONS->TimestepsFileName()); // read timesteps from file +BaseBinaryStar.cpp: dt = std::min(m_Star1->CalculateTimestep(), m_Star2->CalculateTimestep()) * OPTIONS->TimestepMultiplier() / 1000.0; // calculate timestep - make first step small +BaseBinaryStar.cpp: else if (!OPTIONS->EvolveDoubleWhiteDwarfs() && IsWDandWD()) { // double WD and their evolution is not enabled? +BaseBinaryStar.cpp: else if (HasOneOf({ STELLAR_TYPE::MASSLESS_REMNANT }) && !OPTIONS->EvolveMainSequenceMergers()) { // at least one massless remnant and not evolving MS merger products? +BaseBinaryStar.cpp: if (m_Star1->IsOneOf(MAIN_SEQUENCE) && m_Star2->IsOneOf(MAIN_SEQUENCE) && OPTIONS->EvolveMainSequenceMergers()) // yes - both MS and evolving MS merger products? +BaseBinaryStar.cpp: else if (IsUnbound() && !OPTIONS->EvolveUnboundSystems()) { // binary is unbound and we don't want unbound systems? +BaseBinaryStar.cpp: if (OPTIONS->RLOFPrinting()) StashRLOFProperties(MASS_TRANSFER_TIMING::PRE_MT); // stash properties immediately pre-Mass Transfer +BaseBinaryStar.cpp: if (m_Star1->IsOneOf(MAIN_SEQUENCE) && m_Star2->IsOneOf(MAIN_SEQUENCE) && OPTIONS->EvolveMainSequenceMergers()) // yes - both MS and evolving MS merger products? +BaseBinaryStar.cpp: if (!OPTIONS->EvolveUnboundSystems() || IsDCO()) { // should we evolve unbound systems? +BaseBinaryStar.cpp: if (!(OPTIONS->EvolvePulsars() && HasOneOf({ STELLAR_TYPE::NEUTRON_STAR }))) { // evolve pulsar? +BaseBinaryStar.cpp: else if (m_Time > OPTIONS->MaxEvolutionTime()) { // evolution time exceeds maximum? +BaseBinaryStar.cpp: else if (!OPTIONS->EvolveDoubleWhiteDwarfs() && IsWDandWD()) { // double WD and their evolution is not enabled? +BaseBinaryStar.cpp: else if (HasOneOf({ STELLAR_TYPE::MASSLESS_REMNANT }) && !OPTIONS->EvolveMainSequenceMergers()) { // at least one massless remnant and not evolving MS merger products? +BaseBinaryStar.cpp: if (m_Star1->IsOneOf(MAIN_SEQUENCE) && m_Star2->IsOneOf(MAIN_SEQUENCE) && OPTIONS->EvolveMainSequenceMergers()) // yes - both MS and evolving MS merger products? +BaseBinaryStar.cpp: else if (IsUnbound() && !OPTIONS->EvolveUnboundSystems()) { // binary is unbound and we don't want unbound systems? +BaseBinaryStar.cpp: if (stepNum >= OPTIONS->MaxNumberOfTimestepIterations()) evolutionStatus = EVOLUTION_STATUS::STEPS_UP; // number of timesteps for evolution exceeds maximum +BaseBinaryStar.cpp: dt = std::min(m_Star1->CalculateTimestep(), m_Star2->CalculateTimestep()) * OPTIONS->TimestepMultiplier(); // calculate new timestep +BaseBinaryStar.h: return OPTIONS->SwitchLog() ? // switch logging enabled? +BaseBinaryStar.h: return OPTIONS->DetailedOutput() ? LOGGING->LogBSEDetailedOutput(this, p_Id, p_RecordType) : true; +BaseBinaryStar.h: return OPTIONS->EvolvePulsars() ? LOGGING->LogBSEPulsarEvolutionParameters(this, p_RecordType) : true; +BaseBinaryStar.h: return (OPTIONS->RocketKickMagnitude1() > 0) || (OPTIONS->RocketKickMagnitude2() > 0); +BaseStar.cpp: m_BaryonicMassOfMaximumNeutronStarMass = (0.075 * OPTIONS->MaximumNeutronStarMass() * OPTIONS->MaximumNeutronStarMass()) + OPTIONS->MaximumNeutronStarMass(); +BaseStar.cpp: std::tie(ok, value) = OPTIONS->OptionValue(p_Property); +BaseStar.cpp: if (OPTIONS->CommonEnvelopeLambdaNanjingUseRejuvenatedMass()) mass = m_Mass0; // Use rejuvenated mass to calculate lambda instead of true birth mass +BaseStar.cpp: if (OPTIONS->CommonEnvelopeLambdaNanjingEnhanced()) { // If using enhanced Nanjing lambdas +BaseStar.cpp: if (OPTIONS->CommonEnvelopeLambdaNanjingInterpolateInMass()) { +BaseStar.cpp: if (OPTIONS->CommonEnvelopeLambdaNanjingInterpolateInMetallicity()) { +BaseStar.cpp: if (OPTIONS->CommonEnvelopeLambdaNanjingInterpolateInMetallicity()) { +BaseStar.cpp: ZETA_PRESCRIPTION zetaPrescription = OPTIONS->StellarZetaPrescription(); +BaseStar.cpp: QCRIT_PRESCRIPTION qCritPrescription = OPTIONS->QCritPrescription(); +BaseStar.cpp: m_Lambdas.fixed = OPTIONS->CommonEnvelopeLambda(); +BaseStar.cpp: m_Lambdas.kruckow = CalculateLambdaKruckow(m_Radius, OPTIONS->CommonEnvelopeSlopeKruckow()); +BaseStar.cpp: switch (OPTIONS->MassTransferRejuvenationPrescription()) { // which prescription +BaseStar.cpp: return OPTIONS->LuminousBlueVariableFactor() * 1.0E-4; +BaseStar.cpp: rate = OPTIONS->WolfRayetFactor() * 1.0E-13 * PPOW(m_Luminosity, 1.5) * PPOW(m_Metallicity / ZSOL, 0.86) * (1.0 - p_Mu); +BaseStar.cpp: Mdot = PPOW(10.0, logMdot) * OPTIONS->WolfRayetFactor(); +BaseStar.cpp: double LBVRate = CalculateMassLossRateLBV(OPTIONS->LuminousBlueVariablePrescription()); // start with LBV winds (can be, and is often, 0.0) +BaseStar.cpp: OPTIONS->LuminousBlueVariablePrescription() == LBV_PRESCRIPTION::HURLEY_ADD ) { // check whether we should add other winds to the LBV winds (always for HURLEY_ADD prescription, only if not in LBV regime for others) +BaseStar.cpp: otherWindsRate = CalculateMassLossRateHurley() * OPTIONS->CoolWindMassLossMultiplier(); // Apply cool wind mass loss multiplier +BaseStar.cpp: double LBVRate = CalculateMassLossRateLBV(OPTIONS->LuminousBlueVariablePrescription()); // start with LBV winds (can be, and is often, 0.0) +BaseStar.cpp: OPTIONS->LuminousBlueVariablePrescription() == LBV_PRESCRIPTION::HURLEY_ADD ) { // check whether we should add other winds to the LBV winds (always for HURLEY_ADD prescription, only if not in LBV regime for others) +BaseStar.cpp: otherWindsRate = CalculateMassLossRateRSG(OPTIONS->RSGMassLoss()); +BaseStar.cpp: otherWindsRate = CalculateMassLossRateHurley() * OPTIONS->CoolWindMassLossMultiplier(); // apply cool wind mass loss multiplier +BaseStar.cpp: otherWindsRate = CalculateMassLossRateVMS(OPTIONS->VMSMassLoss()); +BaseStar.cpp: otherWindsRate = CalculateMassLossRateOB(OPTIONS->OBMassLoss()); +BaseStar.cpp: * Calls relevant mass loss function based on mass loss prescription given in program options (OPTIONS->massLossPrescription) +BaseStar.cpp: if (OPTIONS->UseMassLoss()) { +BaseStar.cpp: switch (OPTIONS->MassLossPrescription()) { // which prescription? +BaseStar.cpp: mDot = mDot * OPTIONS->OverallWindMassLossMultiplier(); // apply overall wind mass loss multiplier +BaseStar.cpp: if (OPTIONS->UseMassLoss()) { // only if using mass loss (program option) +BaseStar.cpp: if (OPTIONS->CheckPhotonTiringLimit()) { +BaseStar.cpp: if (OPTIONS->UseMassLoss()) { +BaseStar.cpp: switch (OPTIONS->MassTransferAccretionEfficiencyPrescription()) { +BaseStar.cpp: acceptanceRate = min(OPTIONS->MassTransferCParameter() * p_AccretorMassRate, p_DonorMassRate); +BaseStar.cpp: fractionAccreted = OPTIONS->MassTransferFractionAccreted(); +BaseStar.cpp: switch( OPTIONS->MassTransferThermallyLimitedVariation() ) { +BaseStar.cpp: switch (OPTIONS->RotationalVelocityDistribution()) { // which prescription? +BaseStar.cpp: if (utils::Compare(p_RemnantMass, OPTIONS->MaximumNeutronStarMass()) < 0) { +BaseStar.cpp: muKick = max(OPTIONS->MullerMandelKickMultiplierNS() * (p_COCoreMass - p_RemnantMass) / p_RemnantMass, 0.0); +BaseStar.cpp: muKick = max(OPTIONS->MullerMandelKickMultiplierBH() * (p_COCoreMass - p_RemnantMass) / p_RemnantMass, 0.0); +BaseStar.cpp: remnantKick = muKick * (1.0 + gsl_cdf_gaussian_Pinv(rand, OPTIONS->MullerMandelSigmaKick())); +BaseStar.cpp: switch (OPTIONS->KickMagnitudeDistribution()) { // which distribution +BaseStar.cpp: kickMagnitude = DrawKickMagnitudeDistributionFlat(OPTIONS->KickMagnitudeDistributionMaximum(), p_Rand); +BaseStar.cpp: return kickMagnitude / OPTIONS->KickScalingFactor(); +BaseStar.cpp: sigma = OPTIONS->KickMagnitudeDistributionSigmaForECSN(); +BaseStar.cpp: sigma = OPTIONS->KickMagnitudeDistributionSigmaForUSSN(); +BaseStar.cpp: if( utils::SNEventType(m_SupernovaDetails.events.current) == SN_EVENT::PPISN && !OPTIONS->NatalKickForPPISN() ) { +BaseStar.cpp: sigma = OPTIONS->KickMagnitudeDistributionSigmaCCSN_NS(); +BaseStar.cpp: sigma = OPTIONS->KickMagnitudeDistributionSigmaCCSN_BH(); +BaseStar.cpp: if( OPTIONS->RadialChangeFraction()!=0 && radialExpansionTimescale > 0.0 ) // if radial expansion timescale was computed +BaseStar.cpp: dt = min(dt, OPTIONS->RadialChangeFraction()*radialExpansionTimescale); +BaseStar.cpp: if( OPTIONS->MassChangeFraction()!=0 && massChangeTimescale > 0.0 ) // if mass change timescale was computed +BaseStar.cpp: dt = min(dt, OPTIONS->MassChangeFraction()*massChangeTimescale); +BaseStar.h: return OPTIONS->DetailedOutput() ? LOGGING->LogSSEDetailedOutput(this, p_Id, p_RecordType) : true; // Write record to SSE Detailed Output log file +BaseStar.h: return OPTIONS->SwitchLog() ? (LOGGING->ObjectSwitchingPersistence() == OBJECT_PERSISTENCE::PERMANENT ? LOGGING->LogSSESwitchLog(this) : true) : true; // Write record to SSE Switchlog log file +BaseStar.h: virtual double CalculateEddingtonCriticalRate() const { return 2.08E-3 / 1.7 * m_Radius * MYR_TO_YEAR * OPTIONS->EddingtonAccretionFactor() ; } // Hurley+, 2002, Eq. (67) +BH.cpp: switch (OPTIONS->NeutrinoMassLossAssumptionBH()) { // which assumption? +BH.cpp: gravitationalMass = p_BaryonicMass * (1.0 - OPTIONS->NeutrinoMassLossValueBH()); +BH.cpp: gravitationalMass = p_BaryonicMass - OPTIONS->NeutrinoMassLossValueBH(); +BH.cpp: switch (OPTIONS->BlackHoleKicks()) { // which BH kicks option specified? +BinaryConstituentStar.cpp: switch (OPTIONS->CommonEnvelopeMassAccretionPrescription()) { // which prescription? +BinaryConstituentStar.cpp: deltaMass = OPTIONS->CommonEnvelopeMassAccretionConstant(); // use program option +BinaryConstituentStar.cpp: deltaMass = RAND->Random(OPTIONS->CommonEnvelopeMassAccretionMin(), OPTIONS->CommonEnvelopeMassAccretionMax()); // uniform random distribution - Oslowski+ (2011) +BinaryConstituentStar.cpp: deltaMass = std::min(OPTIONS->CommonEnvelopeMassAccretionMax(), std::max(OPTIONS->CommonEnvelopeMassAccretionMin(), m * p_CompanionRadius + c)); +BinaryConstituentStar.cpp: switch (OPTIONS->CommonEnvelopeLambdaPrescription()) { // which common envelope lambda prescription? +BinaryConstituentStar.cpp: m_CEDetails.lambda *= OPTIONS->CommonEnvelopeLambdaMultiplier(); // multiply by constant (program option, default = 1.0) +BinaryConstituentStar.cpp: if ( OPTIONS->CommonEnvelopeFormalism() == CE_FORMALISM::TWO_STAGE ) +changelog.h:// - removed redundant OPTIONS->MassTransferCriticalMassRatioHeliumGiant() from qcritflag if statement in BaseBinaryStar::CalculateMassTransfer() +changelog.h:// - fixed OPTIONS->FixedMetallicity() - always returned true, now returns actual value +changelog.h:// - fixed OPTIONS->OutputPathString() - was always returning raw option instead of fully qualified path +changelog.h:// (a) m_JLoss = OPTIONS->MassTransferJloss(); +changelog.h:// (b) m_FractionAccreted = OPTIONS->MassTransferFractionAccreted(); +changelog.h:// - added OPTIONS->ZetaAdiabaticArbitrary() - option existed, but Options code had no function to retrieve value +changelog.h:// - added OPTIONS->MassTransferFractionAccreted() to options - erroneously not ported from legacy code +changelog.h:// - added ProgramOptionDetails() to Options.cpp and OPTIONS->OptionsDetails() in preparation for change in output functionality +changelog.h:// - fixed issue #162 OPTIONS->UseFixedUK() always returns FALSE. Now returns TRUE if user supplies a fixed kick velocity via --fix-dimensionless-kick-velocity command line option +changelog.h:// - OPTIONS->UseFixedUK() returns TRUE when user supplies -ve value via --fix-dimensionless-kick-velocity. Now return TRUE iff the user supplies a value >=0 via --fix-dimensionless-kick-velocity +changelog.h:// - Removed m_LBVfactor variable from BaseStar - use OPTIONS->LuminousBlueVariableFactor() +changelog.h:// - Removed m_LBVfactor variable from BaseStar - use OPTIONS->WolfRayetFactor() +changelog.h:// - Removed variable 'alpha' from BinaryCEDetails struct - use OPTIONS->CommonEnvelopeAlpha() +changelog.h:// - removed the fixed constant MULLERMANDEL_MAXNS; instead, OPTIONS->MaximumNeutronStarMass() is used for consistency (see issue #1114) +CHeB.cpp: return (OPTIONS->CommonEnvelopeAlphaThermal() * lambdaBG[0]) + ((1.0 - OPTIONS->CommonEnvelopeAlphaThermal()) * lambdaBG[1]); +CHeB.cpp: return (OPTIONS->CommonEnvelopeAlphaThermal() * lambdaBG[0]) + ((1.0 - OPTIONS->CommonEnvelopeAlphaThermal()) * lambdaBG[1]); +CHeB.cpp: switch (OPTIONS->EnvelopeStatePrescription()) { // which envelope prescription? +CHeB.cpp: envelope = utils::Compare(Temperature() * TSOL, OPTIONS->ConvectiveEnvelopeTemperatureThreshold()) > 0 ? ENVELOPE::RADIATIVE : ENVELOPE::CONVECTIVE; // Envelope is radiative if temperature exceeds fixed threshold, otherwise convective +CHeB.h: bool ShouldEnvelopeBeExpelledByPulsations() const { return ( OPTIONS->ExpelConvectiveEnvelopeAboveLuminosityThreshold() && DetermineEnvelopeType() == ENVELOPE::CONVECTIVE && utils::Compare( log10(m_Luminosity/m_Mass), OPTIONS->LuminosityToMassThreshold() ) >= 0 ) ; } // Envelope of convective star with luminosity to mass ratio beyond threshold should be expelled +CH.h: bool ShouldEvolveOnPhase() const { return m_Age < m_Timescales[static_cast(TIMESCALE::tMS)] && (OPTIONS->OptimisticCHE() || m_Omega >= m_OmegaCHE); } // Evolve on CHE phase if age in MS timescale and spinning at least as fast as CHE threshold +EAGB.cpp: return (OPTIONS->CommonEnvelopeAlphaThermal() * lambdaBG[0]) + ((1.0 - OPTIONS->CommonEnvelopeAlphaThermal()) * lambdaBG[1]); +EAGB.cpp: return (OPTIONS->CommonEnvelopeAlphaThermal() * lambdaBG[0]) + ((1.0 - OPTIONS->CommonEnvelopeAlphaThermal()) * lambdaBG[1]); +Errors.cpp: if (p_Prefix == WARNING_PREFIX && !OPTIONS->EnableWarnings()) return false; // do nothing +GiantBranch.cpp: qCrit = OPTIONS->MassTransferCriticalMassRatioGiantDegenerateAccretor(); +GiantBranch.cpp: qCrit = OPTIONS->MassTransferCriticalMassRatioGiantNonDegenerateAccretor(); +GiantBranch.cpp: zeta = OPTIONS->ZetaAdiabaticArbitrary(); +GiantBranch.cpp: zeta = OPTIONS->ZetaRadiativeEnvelopeGiant(); +GiantBranch.cpp: if (utils::Compare(p_COCoreMass, MULLERMANDEL_M1) < 0 || utils::Compare(p_HeCoreMass, OPTIONS->MaximumNeutronStarMass()) <= 0 ) { +GiantBranch.cpp: while (remnantMassMaximumNeutronStarMass() || remnantMass > p_HeCoreMass) { +GiantBranch.cpp: while (remnantMass < MULLERMANDEL_MINNS || remnantMass > OPTIONS->MaximumNeutronStarMass() || remnantMass > p_HeCoreMass) { +GiantBranch.cpp: while (remnantMass < MULLERMANDEL_MINNS || remnantMass > OPTIONS->MaximumNeutronStarMass() || remnantMass > p_HeCoreMass) { +GiantBranch.cpp: while (remnantMass < MULLERMANDEL_MINNS || remnantMass > OPTIONS->MaximumNeutronStarMass() || remnantMass > p_HeCoreMass) { +GiantBranch.cpp: switch (OPTIONS->FryerSupernovaEngine()) { // which SN_ENGINE? +GiantBranch.cpp: baryonicRemnantMass = 1.2 + 0.05 * OPTIONS->Fryer22fmix() + 0.01 * pow( (p_COCoreMass/OPTIONS->Fryer22fmix()), 2.0) + exp( OPTIONS->Fryer22fmix() * (p_COCoreMass - OPTIONS->Fryer22Mcrit()) ) ; // equation 5. +GiantBranch.cpp: switch (OPTIONS->FryerSupernovaEngine()) { // which SN_ENGINE? +GiantBranch.cpp: switch (OPTIONS->RemnantMassPrescription()) { // which prescription? +GiantBranch.cpp: if (OPTIONS->RemnantMassPrescription() == REMNANT_MASS_PRESCRIPTION::MULLER2016) { +GiantBranch.cpp: else if (OPTIONS->RemnantMassPrescription() == REMNANT_MASS_PRESCRIPTION::MULLERMANDEL) { +GiantBranch.cpp: if (utils::Compare(m_Mass, OPTIONS->MaximumNeutronStarMass() ) > 0) +GiantBranch.cpp: else if (OPTIONS->RemnantMassPrescription() == REMNANT_MASS_PRESCRIPTION::HURLEY2000) { +GiantBranch.cpp: else if (utils::Compare(m_Mass, OPTIONS->MaximumNeutronStarMass()) > 0) { +GiantBranch.cpp: if (!m_MassTransferDonorHistory.empty() || (OPTIONS->AllowNonStrippedECSN())) { // If progenitor has never been a MT donor, is it allowed to ECSN? +GiantBranch.cpp: switch (OPTIONS->PulsationalPairInstabilityPrescription()) { // which prescription? +GiantBranch.cpp: if ( OPTIONS->UsePulsationalPairInstability() && +GiantBranch.cpp: utils::Compare(m_HeCoreMass, OPTIONS->PulsationalPairInstabilityLowerLimit()) >= 0 && +GiantBranch.cpp: utils::Compare(m_HeCoreMass, OPTIONS->PulsationalPairInstabilityUpperLimit()) <= 0) { // Pulsational Pair Instability Supernova +GiantBranch.cpp: else if ( OPTIONS->UsePairInstabilitySupernovae() && +GiantBranch.cpp: utils::Compare(m_HeCoreMass, OPTIONS->PairInstabilityLowerLimit()) >= 0 && +GiantBranch.cpp: utils::Compare(m_HeCoreMass, OPTIONS->PairInstabilityUpperLimit()) <= 0) { // Pair Instability Supernova +GiantBranch.cpp: if (OPTIONS->EvolutionMode() == EVOLUTION_MODE::SSE && m_ObjectPersistence == OBJECT_PERSISTENCE::PERMANENT) { +HeGB.cpp: ? OPTIONS->MassTransferCriticalMassRatioHeliumGiantDegenerateAccretor() // degenerate accretor +HeGB.cpp: : OPTIONS->MassTransferCriticalMassRatioHeliumGiantNonDegenerateAccretor(); // non-degenerate accretor +HeHG.cpp: switch (OPTIONS->MassTransferRejuvenationPrescription()) { // which rejuvenation prescription? +HeHG.cpp: switch (OPTIONS->EnvelopeStatePrescription()) { // which envelope prescription? +HeHG.cpp: envelope = utils::Compare(Temperature() * TSOL, OPTIONS->ConvectiveEnvelopeTemperatureThreshold()) > 0 ? ENVELOPE::RADIATIVE : ENVELOPE::CONVECTIVE; // Envelope is radiative if temperature exceeds fixed threshold, otherwise convective +HeHG.cpp: ? OPTIONS->MassTransferCriticalMassRatioHeliumHGDegenerateAccretor() // degenerate accretor +HeHG.cpp: : OPTIONS->MassTransferCriticalMassRatioHeliumHGNonDegenerateAccretor(); // non-degenerate accretor +HeHG.cpp: stellarType = (utils::Compare(m_COCoreMass, OPTIONS->MCBUR1() ) < 0) ? STELLAR_TYPE::CARBON_OXYGEN_WHITE_DWARF : STELLAR_TYPE::OXYGEN_NEON_WHITE_DWARF; +HeMS.cpp: switch (OPTIONS->MassTransferRejuvenationPrescription()) { +HeMS.cpp: if (OPTIONS->WRMassLoss() == WR_MASS_LOSS::SANDERVINK2023) { +HeMS.cpp: else if (OPTIONS->WRMassLoss() == WR_MASS_LOSS::SHENAR2019) { +HeMS.cpp: ? OPTIONS->MassTransferCriticalMassRatioHeliumMSDegenerateAccretor() // degenerate accretor +HeMS.cpp: : OPTIONS->MassTransferCriticalMassRatioHeliumMSNonDegenerateAccretor(); // non-degenerate accretor +HeMS.h: double CalculateZetaConstantsByEnvelope(ZETA_PRESCRIPTION p_ZetaPrescription) { return OPTIONS->ZetaMainSequence(); } // A HeMS star is treated as any other MS star for Zeta calculation purposes +HG.cpp: return (OPTIONS->CommonEnvelopeAlphaThermal() * lambdaBG[0]) + ((1.0 - OPTIONS->CommonEnvelopeAlphaThermal()) * lambdaBG[1]); +HG.cpp: return (OPTIONS->CommonEnvelopeAlphaThermal() * lambdaBG[0]) + ((1.0 - OPTIONS->CommonEnvelopeAlphaThermal()) * lambdaBG[1]); +HG.cpp: switch (OPTIONS->MassTransferRejuvenationPrescription()) { // which rejuvenation prescription? +HG.cpp: ? OPTIONS->MassTransferCriticalMassRatioHGDegenerateAccretor() // degenerate accretor +HG.cpp: : OPTIONS->MassTransferCriticalMassRatioHGNonDegenerateAccretor(); // non-degenerate accretor +HG.cpp: switch (OPTIONS->EnvelopeStatePrescription()) { // which envelope prescription? +HG.cpp: envelope = utils::Compare(Temperature() * TSOL, OPTIONS->ConvectiveEnvelopeTemperatureThreshold()) > 0 ? ENVELOPE::RADIATIVE : ENVELOPE::CONVECTIVE; // Envelope is radiative if temperature exceeds fixed threshold, otherwise convective +Log.cpp: size_t IOBufSize = OPTIONS->HDF5BufferSize() * chunkSize; // IO buffer size +Log.cpp: m_OptionDetails = OPTIONS->CmdLineOptionsDetails(); // get commandline option details +Log.cpp: string fileExt = "." + LOGFILETYPEFileExt.at(OPTIONS->LogfileType()); // file extension for HDF5 files +Log.cpp: if (OPTIONS->StoreInputFiles()) { // user wants input files stored in output container? +Log.cpp: if (!OPTIONS->GridFilename().empty()) { // user specified a grid file? +Log.cpp: boost::filesystem::path srcPath(OPTIONS->GridFilename()); // grid file fully-qualified name +Log.cpp: (void)boost::filesystem::copy_file(OPTIONS->GridFilename(), dstFn, BOOST_OVERWRITE_EXISTING); // copy grid file - overwrite any existing file (shouldn't be one, but just in case we want this one) +Log.cpp: Squawk("ERROR: Unable to copy grid file " + OPTIONS->GridFilename() + " to output container " + dstPath); // announce error +Log.cpp: if (m_Enabled && !OPTIONS->LogfileDefinitionsFilename().empty()) { // user specified a logfile-definitions file? +Log.cpp: boost::filesystem::path srcPath(OPTIONS->LogfileDefinitionsFilename()); // logfile-definitions file fully-qualified name +Log.cpp: (void)boost::filesystem::copy_file(OPTIONS->LogfileDefinitionsFilename(), dstFn, BOOST_OVERWRITE_EXISTING); // copy logfile-definitions file - overwrite any existing file (shouldn't be one, but just in case we want this one) +Log.cpp: Squawk("ERROR: Unable to copy logfile-definitions file " + OPTIONS->LogfileDefinitionsFilename() + " to output container " + dstPath); // announce error +Log.cpp: unsigned long int actualRandomSeed = OPTIONS->FixedRandomSeedCmdLine() ? OPTIONS->RandomSeedCmdLine() : RAND->DefaultSeed(); // actual random seed used +Log.cpp: if (objectsRequested >= 0 && (int)OPTIONS->nObjectsToEvolve() == objectsRequested) derivation = "USER_SUPPLIED"; // should be right most of the time (not critical) +Log.cpp: if (OPTIONS->PrintBoolAsString()) +Log.cpp: if (OPTIONS->EvolutionMode() == EVOLUTION_MODE::SSE) +Log.cpp: if (OPTIONS->EvolutionMode() == EVOLUTION_MODE::SSE) { +Log.cpp: string fileext = LOGFILETYPEFileExt.at(OPTIONS->LogfileType()); // file extension +Log.cpp: details = std::make_tuple(TYPENAME::STRING, OPTIONS->NotesHdrs(p_Idx), "-", 0, 1); // yes - construct details +Log.cpp: if ((OPTIONS->AddOptionsToSysParms() == ADD_OPTIONS_TO_SYSPARMS::ALWAYS) || // always add option columns? +Log.cpp: ((OPTIONS->AddOptionsToSysParms() == ADD_OPTIONS_TO_SYSPARMS::GRID) && // add for grids? +Log.cpp: (!OPTIONS->GridFilename().empty() || OPTIONS->CommandLineGrid()))) { // have grid file or ranges/sets? +Log.cpp: if ((OPTIONS->AddOptionsToSysParms() == ADD_OPTIONS_TO_SYSPARMS::ALWAYS) || // always add option columns? +Log.cpp: ((OPTIONS->AddOptionsToSysParms() == ADD_OPTIONS_TO_SYSPARMS::GRID) && // add for grids? +Log.cpp: (!OPTIONS->GridFilename().empty() || OPTIONS->CommandLineGrid()))) { // have grid file or ranges/sets? +Log.cpp: if (OPTIONS->PrintBoolAsString()) { // print bool values as strings "TRUE" or "FALSE"? +Log.cpp: fileDetails.filename = OPTIONS->LogfileBeBinaries(); +Log.cpp: fileDetails.recordTypes = OPTIONS->LogfileBeBinariesRecordTypes(); +Log.cpp: fileDetails.filename = OPTIONS->LogfileCommonEnvelopes(); +Log.cpp: fileDetails.recordTypes = OPTIONS->LogfileCommonEnvelopesRecordTypes(); +Log.cpp: fileDetails.filename = OPTIONS->LogfileDoubleCompactObjects(); +Log.cpp: fileDetails.recordTypes = OPTIONS->LogfileDoubleCompactObjectsRecordTypes(); +Log.cpp: fileDetails.filename = OPTIONS->LogfilePulsarEvolution(); +Log.cpp: fileDetails.recordTypes = OPTIONS->LogfilePulsarEvolutionRecordTypes(); +Log.cpp: fileDetails.filename = OPTIONS->LogfileRLOFParameters(); +Log.cpp: fileDetails.recordTypes = OPTIONS->LogfileRLOFParametersRecordTypes(); +Log.cpp: fileDetails.filename = OPTIONS->LogfileSupernovae(); +Log.cpp: fileDetails.recordTypes = OPTIONS->LogfileSupernovaeRecordTypes(); +Log.cpp: fileDetails.filename = OPTIONS->LogfileSwitchLog(); +Log.cpp: fileDetails.filename = OPTIONS->LogfileSystemParameters(); +Log.cpp: fileDetails.recordTypes = OPTIONS->LogfileSystemParametersRecordTypes(); +Log.cpp: if ((OPTIONS->AddOptionsToSysParms() == ADD_OPTIONS_TO_SYSPARMS::ALWAYS) || // always add option columns? +Log.cpp: ((OPTIONS->AddOptionsToSysParms() == ADD_OPTIONS_TO_SYSPARMS::GRID) && // add for grids? +Log.cpp: (!OPTIONS->GridFilename().empty() || OPTIONS->CommandLineGrid()))) { // have grid file or ranges/sets? +Log.cpp: fileDetails.filename = OPTIONS->LogfileSupernovae(); +Log.cpp: fileDetails.recordTypes = OPTIONS->LogfileSupernovaeRecordTypes(); +Log.cpp: fileDetails.filename = OPTIONS->LogfileSwitchLog(); +Log.cpp: fileDetails.filename = OPTIONS->LogfileSystemParameters(); +Log.cpp: fileDetails.recordTypes = OPTIONS->LogfileSystemParametersRecordTypes(); +Log.cpp: if ((OPTIONS->AddOptionsToSysParms() == ADD_OPTIONS_TO_SYSPARMS::ALWAYS) || // always add option columns? +Log.cpp: ((OPTIONS->AddOptionsToSysParms() == ADD_OPTIONS_TO_SYSPARMS::GRID) && // add for grids? +Log.cpp: (!OPTIONS->GridFilename().empty() || OPTIONS->CommandLineGrid()))) { // have grid file or ranges/sets? +Log.cpp: fileDetails.filename = DETAILED_OUTPUT_DIRECTORY_NAME + "/" + OPTIONS->LogfileDetailedOutput(); // logfile filename with directory +Log.cpp: fileDetails.recordTypes = OPTIONS->LogfileDetailedOutputRecordTypes(); // record types +Log.cpp: string fileExt = "." + LOGFILETYPEFileExt.at(OPTIONS->LogfileType()); // file extension for HDF5 files +Log.cpp: // program option NOTES is special - the property is actually a vector of strings (OPTIONS->Notes()), +Log.cpp: if (OPTIONS->LogfileType() == LOGFILETYPE::HDF5) { // logging to HDF5 files? +Log.cpp: size_t chunkSize = OPTIONS->nObjectsToEvolve() < HDF5_MINIMUM_CHUNK_SIZE || +Log.cpp: p_Logfile == LOGFILE::BSE_DETAILED_OUTPUT ? HDF5_MINIMUM_CHUNK_SIZE : OPTIONS->HDF5ChunkSize(); // chunk size +Log.cpp: size_t IOBufSize = OPTIONS->HDF5BufferSize() * chunkSize; // IO buffer size +Log.cpp: string filename = OPTIONS->LogfileDefinitionsFilename(); // get user-specified definitions file +Log.cpp: std::vector addNotes = std::vector(OPTIONS->NotesHdrs().size(), false); // annotations (notes) user wants added to the base properties +Log.cpp: std::vector subtractNotes = std::vector(OPTIONS->NotesHdrs().size(), false); // annotations (notes) user wants subtracted from the base properties +Log.cpp: addNotes = std::vector(OPTIONS->NotesHdrs().size(), false); // start with no annotations to be added +Log.cpp: subtractNotes = std::vector(OPTIONS->NotesHdrs().size(), false); // start with no annotations to be subtracted +Log.cpp: if (notesIdx < 1 or notesIdx > (int)OPTIONS->NotesHdrs().size()) { // index in valid range? +Log.h: string fmt = OPTIONS->PrintBoolAsString() ? "%5s" : "%1s"; +Log.h: string vS = OPTIONS->PrintBoolAsString() ? (v ? "TRUE " : "FALSE") : (v ? "1" : "0"); +Log.h: string fmt = OPTIONS->PrintBoolAsString() ? "%5s" : "%1s"; +Log.h: string vS = OPTIONS->PrintBoolAsString() ? (v ? "TRUE " : "FALSE") : (v ? "1" : "0"); +Log.h: // user-specified option '--notes-hdrs' (OPTIONS->NotesHdrs()) +Log.h: std::vector m_BSE_BE_Binaries_Notes = std::vector(OPTIONS->NotesHdrs().size(), false); +Log.h: std::vector m_BSE_CEE_Notes = std::vector(OPTIONS->NotesHdrs().size(), false); +Log.h: std::vector m_BSE_DCO_Notes = std::vector(OPTIONS->NotesHdrs().size(), false); +Log.h: std::vector m_BSE_Detailed_Notes = std::vector(OPTIONS->NotesHdrs().size(), false); +Log.h: std::vector m_BSE_Pulsars_Notes = std::vector(OPTIONS->NotesHdrs().size(), false); +Log.h: std::vector m_BSE_RLOF_Notes = std::vector(OPTIONS->NotesHdrs().size(), false); +Log.h: std::vector m_BSE_SNE_Notes = std::vector(OPTIONS->NotesHdrs().size(), false); +Log.h: std::vector m_BSE_Switch_Notes = std::vector(OPTIONS->NotesHdrs().size(), false); +Log.h: std::vector m_BSE_SysParms_Notes = std::vector(OPTIONS->NotesHdrs().size(), false); +Log.h: std::vector m_SSE_Detailed_Notes = std::vector(OPTIONS->NotesHdrs().size(), false); +Log.h: std::vector m_SSE_SNE_Notes = std::vector(OPTIONS->NotesHdrs().size(), false); +Log.h: std::vector m_SSE_Switch_Notes = std::vector(OPTIONS->NotesHdrs().size(), false); +Log.h: std::vector m_SSE_SysParms_Notes = std::vector(OPTIONS->NotesHdrs().size(), false); +Log.h: bool hdf5 = OPTIONS->LogfileType() == LOGFILETYPE::HDF5; // logging to hdf5 file? +Log.h: switch (OPTIONS->LogfileType()) { +Log.h: value = boost::variant(OPTIONS->Notes(idx)); // yes - get value +Log.h: if (OPTIONS->LogfileType() == LOGFILETYPE::HDF5) // logging to HDF5 file? +Log.h: if (OPTIONS->LogfileType() == LOGFILETYPE::HDF5) { // logging to HDF5 file? +main.cpp: * The signal is raised in the Star::SwitchTo() function if OPTIONS->BSESwitchLog() +main.cpp: if (evolvingBinaryStarValid && OPTIONS->SwitchLog()) { // yes - do we have a valid binary star, and are we logging switches? +main.cpp: bool usingGrid = !OPTIONS->GridFilename().empty(); // using grid file? +main.cpp: // OPTIONS->AdvanceCmdLineOptionValues(), called at the end of the loop, advances the +main.cpp: int gridResult = OPTIONS->ApplyNextGridLine(); // set options according to specified values in grid file +main.cpp: SHOW_ERROR(error, "Accessing grid file '" + OPTIONS->GridFilename() + "'"); // show error +main.cpp: ERROR error = OPTIONS->RewindGridFile(); // ready for next commandline options variation +main.cpp: SHOW_ERROR(error, "Accessing grid file '" + OPTIONS->GridFilename() + "'"); // no - show error (should never happen here - should be picked up at file open) +main.cpp: SHOW_ERROR(ERROR::ERROR, "Accessing grid file '" + OPTIONS->GridFilename() + "'"); // unexpected error - show error +main.cpp: // (by OPTIONS->ApplyNextGridLine()). +main.cpp: // OPTIONS->AdvanceGridLineOptionValues(), called at the end of the loop, advances the +main.cpp: if (OPTIONS->FixedRandomSeedGridLine()) { // user specified a random seed in the grid file for this binary? +main.cpp: randomSeed = OPTIONS->RandomSeedGridLine() + (unsigned long int)gridLineVariation; // random seed +main.cpp: if (OPTIONS->SetRandomSeed(randomSeed, optsOrigin) < 0) { // ok? +main.cpp: else if (OPTIONS->FixedRandomSeedCmdLine()) { // no - user specified a random seed on the commandline? +main.cpp: randomSeed = OPTIONS->RandomSeedCmdLine() + (unsigned long int)index; // random seed +main.cpp: if (OPTIONS->SetRandomSeed(randomSeed, optsOrigin) < 0) { // ok? +main.cpp: if (OPTIONS->SetRandomSeed(randomSeed, optsOrigin) < 0) { // ok? +main.cpp: double initialMass = OPTIONS->OptionSpecified("initial-mass") == 1 // user specified mass? +main.cpp: ? OPTIONS->InitialMass() // yes, use it +main.cpp: : utils::SampleInitialMass(OPTIONS->InitialMassFunction(), // no, sample it +main.cpp: OPTIONS->InitialMassFunctionMax(), +main.cpp: OPTIONS->InitialMassFunctionMin(), +main.cpp: OPTIONS->InitialMassFunctionPower()); +main.cpp: double metallicity = OPTIONS->OptionSpecified("metallicity") == 1 // user specified metallicity? +main.cpp: ? OPTIONS->Metallicity() // yes, use it +main.cpp: : utils::SampleMetallicity(OPTIONS->MetallicityDistribution(), +main.cpp: OPTIONS->MetallicityDistributionMax(), +main.cpp: OPTIONS->MetallicityDistributionMin()); // no, sample it +main.cpp: kickParameters.magnitudeRandomSpecified = OPTIONS->OptionSpecified("kick-magnitude-random") == 1; +main.cpp: kickParameters.magnitudeRandom = OPTIONS->KickMagnitudeRandom(); +main.cpp: kickParameters.magnitudeSpecified = OPTIONS->OptionSpecified("kick-magnitude") == 1; +main.cpp: kickParameters.magnitude = OPTIONS->KickMagnitude(); +main.cpp: star = OPTIONS->OptionSpecified("rotational-frequency") == 1 // user specified rotational frequency? +main.cpp: ? new Star(randomSeed, initialMass, metallicity, kickParameters, OPTIONS->RotationalFrequency() * SECONDS_IN_YEAR) // yes - use it (convert from Hz to cycles per year - see BaseStar::CalculateZAMSAngularFrequency()) +main.cpp: if (!OPTIONS->Quiet()) { // quiet mode? +main.cpp: int optionsStatus = OPTIONS->AdvanceGridLineOptionValues(); // apply next grid file options (ranges/sets) +main.cpp: int optionsStatus = OPTIONS->AdvanceCmdLineOptionValues(); // yes - apply next commandline options (ranges/sets) +main.cpp: if (usingGrid || OPTIONS->CommandLineGrid() || (!usingGrid && index >= OPTIONS->nObjectsToEvolve())) { // created required number of stars? +main.cpp: if (!OPTIONS->Quiet()) { +main.cpp: bool usingGrid = !OPTIONS->GridFilename().empty(); // using grid file? +main.cpp: // OPTIONS->AdvanceCmdLineOptionValues(), called at the end of the loop, advances the +main.cpp: int gridResult = OPTIONS->ApplyNextGridLine(); // yes - set options according to specified values in grid file +main.cpp: SHOW_ERROR(error, "Accessing grid file '" + OPTIONS->GridFilename() + "'"); // show error +main.cpp: ERROR error = OPTIONS->RewindGridFile(); // ready for next commandline options variation +main.cpp: SHOW_ERROR(error, "Accessing grid file '" + OPTIONS->GridFilename() + "'"); // no - show error (should never happen here - should be picked up at file open) +main.cpp: SHOW_ERROR(ERROR::ERROR, "Accessing grid file '" + OPTIONS->GridFilename() + "'"); // unexpected error - show error +main.cpp: // (by OPTIONS->ApplyNextGridLine()). +main.cpp: // OPTIONS->AdvanceGridLineOptionValues(), called at the end of the loop, advances the +main.cpp: if (OPTIONS->FixedRandomSeedGridLine()) { // user specified a random seed in the grid file for this binary? +main.cpp: randomSeed = OPTIONS->RandomSeedGridLine() + (unsigned long int)gridLineVariation; // random seed +main.cpp: if (OPTIONS->SetRandomSeed(randomSeed, optsOrigin) < 0) { // ok? +main.cpp: else if (OPTIONS->FixedRandomSeedCmdLine()) { // no - user specified a random seed on the commandline? +main.cpp: randomSeed = OPTIONS->RandomSeedCmdLine() + (unsigned long int)index + (unsigned long int)gridLineVariation; // random seed +main.cpp: if (OPTIONS->SetRandomSeed(randomSeed, optsOrigin) < 0) { // ok? +main.cpp: if (OPTIONS->SetRandomSeed(randomSeed, optsOrigin) < 0) { // ok? +main.cpp: if (!OPTIONS->Quiet()) { // quiet mode? +main.cpp: if (OPTIONS->CHEMode() == CHE_MODE::NONE) { // CHE enabled? +main.cpp: int optionsStatus = OPTIONS->AdvanceGridLineOptionValues(); // apply next grid file options (ranges/sets) +main.cpp: int optionsStatus = OPTIONS->AdvanceCmdLineOptionValues(); // apply next commandline options (ranges/sets) +main.cpp: if (usingGrid || OPTIONS->CommandLineGrid() || (!usingGrid && index >= OPTIONS->nObjectsToEvolve())) { // created required number of stars? +main.cpp: if (!OPTIONS->Quiet()) { +main.cpp: bool ok = OPTIONS->Initialise(argc, argv); // get the program options from the commandline +main.cpp: if (OPTIONS->RequestedHelp()) { // user requested help? +main.cpp: OPTIONS->ShowHelp(); // show help +main.cpp: else if (OPTIONS->RequestedVersion()) { // user requested version? +main.cpp: else if (!OPTIONS->YAMLfilename().empty()) { // user requested YAML file creation? +main.cpp: yaml::MakeYAMLfile(OPTIONS->YAMLfilename(), OPTIONS->YAMLtemplate()); // create YAML file +main.cpp: LOGGING->Start(OPTIONS->OutputPathString(), // location of logfiles +main.cpp: OPTIONS->OutputContainerName(), // directory to be created for logfiles +main.cpp: OPTIONS->LogfileNamePrefix(), // prefix for logfile names +main.cpp: OPTIONS->LogLevel(), // log level - determines (in part) what is written to log file +main.cpp: OPTIONS->LogClasses(), // log classes - determines (in part) what is written to log file +main.cpp: OPTIONS->DebugLevel(), // debug level - determines (in part) what debug information is displayed +main.cpp: OPTIONS->DebugClasses(), // debug classes - determines (in part) what debug information is displayed +main.cpp: OPTIONS->DebugToFile(), // should debug statements also be written to logfile? +main.cpp: OPTIONS->ErrorsToFile(), // should error messages also be written to logfile? +main.cpp: OPTIONS->LogfileType()); // log file type +main.cpp: if (!OPTIONS->GridFilename().empty()) { // have grid filename? +main.cpp: ERROR error = OPTIONS->OpenGridFile(OPTIONS->GridFilename()); // yes - open grid file +main.cpp: SHOW_ERROR(error, "Accessing grid file '" + OPTIONS->GridFilename() + "'"); // no - show error +main.cpp: if (OPTIONS->EvolutionMode() == EVOLUTION_MODE::SSE) { // SSE? +main.cpp: if (!OPTIONS->GridFilename().empty()) { // have grid filename? +main.cpp: OPTIONS->CloseGridFile(); // yes - close it if it's open +MainSequence.cpp: if (OPTIONS->RetainCoreMassDuringCaseAMassTransfer()) { +MainSequence.h: double CalculateCoreMassAtPhaseEnd() const { return OPTIONS->RetainCoreMassDuringCaseAMassTransfer() ? MinimumCoreMass() : 0.0; } // Accounts for minimal core mass built up prior to mass loss through mass transfer +MainSequence.h: double CalculateZetaConstantsByEnvelope(ZETA_PRESCRIPTION p_ZetaPrescription) { return OPTIONS->ZetaMainSequence(); } +MS_gt_07.cpp: switch (OPTIONS->MassTransferRejuvenationPrescription()) { // which prescription? +MS_gt_07.cpp: ? OPTIONS->MassTransferCriticalMassRatioMSHighMassDegenerateAccretor() // degenerate accretor +MS_gt_07.cpp: : OPTIONS->MassTransferCriticalMassRatioMSHighMassNonDegenerateAccretor(); // non-degenerate accretor +MS_gt_07.cpp: switch (OPTIONS->EnvelopeStatePrescription()) { // which envelope prescription? +MS_gt_07.cpp: envelope = utils::Compare(Temperature() * TSOL, OPTIONS->ConvectiveEnvelopeTemperatureThreshold()) > 0 ? ENVELOPE::RADIATIVE : ENVELOPE::CONVECTIVE; // Envelope is radiative if temperature exceeds fixed threshold, otherwise convective +MS_lte_07.cpp: switch (OPTIONS->MassTransferRejuvenationPrescription()) { // which prescription? +MS_lte_07.cpp: ? OPTIONS->MassTransferCriticalMassRatioMSLowMassDegenerateAccretor() // degenerate accretor +MS_lte_07.cpp: : OPTIONS->MassTransferCriticalMassRatioMSLowMassNonDegenerateAccretor(); // non-degenerate accretor +NS.cpp: switch (OPTIONS->NeutronStarEquationOfState()) { // which equation-of-state? +NS.cpp: switch (OPTIONS->PulsarBirthSpinPeriodDistribution()) { // which distribution? +NS.cpp: double maximum = OPTIONS->PulsarBirthSpinPeriodDistributionMax(); +NS.cpp: double minimum = OPTIONS->PulsarBirthSpinPeriodDistributionMin(); +NS.cpp: switch (OPTIONS->PulsarBirthMagneticFieldDistribution()) { // which distribution? +NS.cpp: double maximum = OPTIONS->PulsarBirthMagneticFieldDistributionMax(); +NS.cpp: double minimum = OPTIONS->PulsarBirthMagneticFieldDistributionMin(); +NS.cpp: double maximum = PPOW(10.0, OPTIONS->PulsarBirthMagneticFieldDistributionMax()); +NS.cpp: double minimum = PPOW(10.0, OPTIONS->PulsarBirthMagneticFieldDistributionMin()); +NS.cpp: double magFieldLowerLimit = PPOW(10.0, OPTIONS->PulsarLog10MinimumMagneticField()) * GAUSS_TO_TESLA; +NS.cpp: double tau = OPTIONS->PulsarMagneticFieldDecayTimescale() * MYR_TO_YEAR * SECONDS_IN_YEAR; +NS.cpp: double magFieldLowerLimit_G = PPOW(10.0, OPTIONS->PulsarLog10MinimumMagneticField()) ; // convert to Gauss +NS.cpp: switch (OPTIONS->NeutronStarAccretionInCE()) { // which mode of CE accretion to use? +NS.cpp: double kappa = OPTIONS->PulsarMagneticFieldDecayMassscale() * MSOL_TO_G; +NS.cpp: double magFieldLowerLimit_G = PPOW(10.0, OPTIONS->PulsarLog10MinimumMagneticField()) ; +NS.h: bool ShouldEvolveOnPhase() const { return (m_Mass <= OPTIONS->MaximumNeutronStarMass()); } // Evolve as a neutron star unless mass > maximum neutron star mass (e.g. through accretion) +Options.cpp: if (OPTIONS->OptionSpecified("grid-lines-to-process") == 1 && // user specified number of grid lines to process? +Options.cpp: m_Gridfile.startLine = OPTIONS->GridStartLine(); // set first line to process (0-based) +Options.cpp: m_Gridfile.linesToProcess = OPTIONS->GridLinesToProcess(); // set number of lines to process (-1 = process to EOF) +Star.cpp: if (OPTIONS->CHEMode() != CHE_MODE::NONE && utils::Compare(m_Star->Omega(), m_Star->OmegaCHE()) >= 0) { // CHE? +Star.cpp: if (utils::IsOneOf(stellarTypePrev, EVOLVABLE_TYPES) && OPTIONS->SwitchLog()) { // star should be evolving from one of the evolvable types (We don't want the initial switch from Star->MS. Not necessary for BSE (handled differently), but no harm) +Star.cpp: if (OPTIONS->EvolutionMode() == EVOLUTION_MODE::BSE) { // BSE? +Star.cpp: if (!OPTIONS->TimestepsFileName().empty()) { // have timesteps filename? +Star.cpp: std::tie(error, timesteps) = utils::ReadTimesteps(OPTIONS->TimestepsFileName()); // read timesteps from file +Star.cpp: if (m_Star->Time() > OPTIONS->MaxEvolutionTime()) { +Star.cpp: else if (stepNum >= OPTIONS->MaxNumberOfTimestepIterations()) { +Star.cpp: dt = m_Star->CalculateTimestep() * OPTIONS->TimestepMultiplier(); // calculate new timestep +TPAGB.cpp: return (OPTIONS->CommonEnvelopeAlphaThermal() * lambdaBG[0]) + ((1.0 - OPTIONS->CommonEnvelopeAlphaThermal()) * lambdaBG[1]); +TPAGB.cpp: return (OPTIONS->CommonEnvelopeAlphaThermal() * lambdaBG[0]) + ((1.0 - OPTIONS->CommonEnvelopeAlphaThermal()) * lambdaBG[1]); +TPAGB.cpp: stellarType = (utils::Compare(gbParams(McBAGB), OPTIONS->MCBUR1() ) < 0) ? STELLAR_TYPE::CARBON_OXYGEN_WHITE_DWARF : STELLAR_TYPE::OXYGEN_NEON_WHITE_DWARF; +TPAGB.h: bool IsSupernova() const { return (utils::Compare(m_COCoreMass, m_GBParams[static_cast(GBP::McSN)]) >= 0 && utils::Compare(CalculateInitialSupernovaMass(), OPTIONS->MCBUR1()) >= 0 && utils::Compare(m_COCoreMass, m_Mass) < 0); } // Going supernova if still has envelope and core mass large enough +typedefs.h: std::vector annotations; // print flags for each annotation specified by the user (e.g. OPTIONS->NotesHdrs()) +typedefs.h: double fixed; // Set to OPTIONS->commonEnvelopeLambda +typedefs.h: double kruckow; // Calculated using m_Radius and OPTIONS->commonEnvelopeSlopeKruckow +typedefs.h: double fixed; // Calculated using lambda = OPTIONS->commonEnvelopeLambda +yaml.cpp: std::vector optionDetails = OPTIONS->CmdLineOptionsDetails(); diff --git a/src/NS.cpp b/src/NS.cpp index c26c7b560..14ca44a64 100755 --- a/src/NS.cpp +++ b/src/NS.cpp @@ -1,6 +1,7 @@ #include "Rand.h" #include "NS.h" +#include "gsl/gsl_poly.h" /* * Calculate the luminosity of a Neutron Star @@ -350,7 +351,6 @@ void NS::CalculateAndSetPulsarParameters() { 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 * SECONDS_IN_MS; // convert from ms to s - m_MomentOfInertia_CGS = CalculateMomentOfInertiaCGS(); // in CGS g cm^2 // Note we convert neutronStarMomentOfInertia from CGS to SI here @@ -396,31 +396,152 @@ void NS::SpinDownIsolatedPulsar(const double p_Stepsize) { // 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 + // calculate the spin frequency 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.0 * constant2 * (term1 - term2 - term3) + (initialSpinPeriod * initialSpinPeriod); + double Psquared = 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 = 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_PulsarDetails.spinDownRate = CalculateSpinDownRate(m_PulsarDetails.spinFrequency,m_MomentOfInertia_CGS , m_PulsarDetails.magneticField, NSradius_IN_CM / KM_TO_CM); m_AngularMomentum_CGS = m_PulsarDetails.spinFrequency * m_MomentOfInertia_CGS; // angular momentum of star in CGS } /* - * Update the magnetic field and spins of neutron stars in the following situations: - * 1). + * Calculate the magnetic field, spin, spin-down rate and angular momentum of a neutron star, + * when it's accreting mass from companion through stable mass transfer. + * The calculations used in this function follow closely Sec. 2.2.1 in arxiv:1912.02415 + * We carry out the calculations in this function using cgs units. + * This function is called in the NS::UpdateMagneticFieldAndSpin() function + * + * Returns, in a tuple, these four quantities of a neutron star: + * + * magnetic field strength + * spin frequency + * spin-down rate + * angular momentum + * + * DBL_DBL_DBL_DBL PulsarAccretion(const double p_MagField, const double p_SpinFrequency, const double p_AngularMomentum, const double p_Stepsize, const double p_MassGainPerTimeStep, const double kappa, const double p_Epsilon) + * @param [IN] p_MagField NS magnetic field strength at the beginning of accretion (in Tesla) + * @param [IN] p_SpinFrequency Spin frequency for the NS at the beginning of accretion (in Hz) + * @param [IN] p_AngularMomentum Angular momentum of the NS at the beginning of accretion (g cm-2 / s) + * @param [IN] p_Stepsize Timestep size for integration (in seconds) + * @param [IN] p_MassGainPerTimeStep Mass transferred from the secondary for each iteration (in g) + * @param [IN] p_Kappa Magnetic field mass decay scale (in g) + * @param [IN] p_Epsilon Efficiency factor allowing for uncertainties of coupling magnetic field and matter. + * @return Tuple containing the updated magnetic field strength, spin frequency, spin-down rate and angular momentum of neutron star + */ +DBL_DBL_DBL_DBL NS::PulsarAccretion(const double p_MagField, const double p_SpinFrequency, const double p_AngularMomentum, const double p_Stepsize, const double p_MassGainPerTimeStep, const double p_Kappa, const double p_Epsilon) { + double initialMagField_G = p_MagField * TESLA_TO_GAUSS; + double magFieldLowerLimit_G = PPOW(10.0, OPTIONS->PulsarLog10MinimumMagneticField()) ; // convert to Gauss + double mass_g = m_Mass * MSOL_TO_G; // in g + double r_cm = m_Radius * RSOL_TO_KM * KM_TO_CM; + double angularMomentum = p_AngularMomentum; + + //Magnetic field decay due to mass transfer. + //Follows Eq. 12 in arxiv:1912.02415 + double newPulsarMagneticField = (initialMagField_G - magFieldLowerLimit_G) * exp(-p_MassGainPerTimeStep / p_Kappa) + magFieldLowerLimit_G; + + // calculate the Alfven radius for an accreting neutron star + // see Equation 10 in arxiv:1912.02415 + double mDot = p_MassGainPerTimeStep / p_Stepsize; + double R_CM_6 = r_cm * r_cm * r_cm * r_cm * r_cm * r_cm; + double p = R_CM_6 * R_CM_6 / (mass_g * mass_g * mDot); + double q = PPOW(p, 1.0/7.0); + double alfvenConst = PPOW(2 * PI_2 / G_CGS, 1.0/7.0) ; + double alfvenRadius = alfvenConst * q * PPOW(initialMagField_G, 4.0/7.0); + double magneticRadius = alfvenRadius / 2.0; + // calculate the difference in the keplerian angular velocity at the magnetic radius + // and surface angular velocity of the neutron star + // magnetic radius is half of alfven radius + // see Equation 2 in 1994MNRAS.269..455J / Equation 9 in arxiv:1912.02415 + double keplerianVelocityAtMagneticRadius = std::sqrt((G_CGS) * mass_g / magneticRadius); + double keplerianAngularVelocityAtMagneticRadius = keplerianVelocityAtMagneticRadius / magneticRadius; + double omegaDifference = keplerianAngularVelocityAtMagneticRadius - p_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 * omegaDifference * magneticRadius * magneticRadius * mDot; + angularMomentum = angularMomentum + Jdot * p_Stepsize; + return std::make_tuple(newPulsarMagneticField * GAUSS_TO_TESLA, angularMomentum / m_MomentOfInertia_CGS, Jdot / m_MomentOfInertia_CGS, angularMomentum); +} + + + +/* + * Calculate the magnetic field, spin, spin-down rate and angular momentum of a neutron star, + * when it's accreting mass from companion through stable mass transfer. + * The calculations used in this function follow closely Sec. 2.2.1 in arxiv:1912.02415 + * We carry out the calculations in this function using cgs units. + * This function is called in the NS::UpdateMagneticFieldAndSpin() function + * + * Returns the updated: + * dJ/dm: change in angular momentum per mass transferred. + * + * DBL_DBL_DBL_DBL PulsarAccretion(const double p_MagField, const double p_SpinFrequency, const double p_AngularMomentum, const double p_Stepsize, const double p_MassGainPerTimeStep, const double kappa, const double p_Epsilon) + * @param [IN] p_MagField NS magnetic field strength at the beginning of accretion (in Gauss) + * @param [IN] p_SpinFrequency Spin frequency for the NS at the beginning of accretion (in Hz) + * @param [IN] p_AngularMomentum Angular momentum of the NS at the beginning of accretion (g cm-2 / s) + * @param [IN] p_Stepsize Timestep size for integration (in seconds) + * @param [IN] p_MassGainPerTimeStep Mass transferred from the secondary for each iteration (in g) + * @param [IN] p_Kappa Magnetic field mass decay scale (in g) + * @param [IN] p_Epsilon Efficiency factor allowing for uncertainties of coupling magnetic field and matter. + * @param [IN] p_MoI Moment of Inertia + * @return Tuple containing the updated magnetic field strength, spin frequency, spin-down rate and angular momentum of neutron star + */ +DBL_DBL_DBL_DBL NS::deltaAngularMomentumByPulsarAccretion(const double p_MassGainPerTimeStep, const double p_Mass, const double p_Radius, const double p_MagField, const double p_SpinFrequency, const double p_AngularMomentum, const double p_Stepsize, const double p_Kappa, const double p_Epsilon, const double p_MoI) { + if (utils::Compare(p_MassGainPerTimeStep, 0.0) <= 0) { + return std::make_tuple(0.0, 0.0, 0.0, 0.0); + } + double initialMagField_G = p_MagField; + double magFieldLowerLimit_G = PPOW(10.0, OPTIONS->PulsarLog10MinimumMagneticField()) ; // convert to Gauss + double mass_g = p_Mass * MSOL_TO_G; // in g + double r_cm = p_Radius * RSOL_TO_KM * KM_TO_CM; + double angularMomentum = p_AngularMomentum; + //std::cout << "in Func: " << p_Kappa <<" " << p_SpinFrequency << " " << mass_g << " " << r_cm << " " << angularMomentum << std::endl; + //Magnetic field decay due to mass transfer. + //Follows Eq. 12 in arxiv:1912.02415 + double newPulsarMagneticField = (initialMagField_G - magFieldLowerLimit_G) * exp(-p_MassGainPerTimeStep / p_Kappa) + magFieldLowerLimit_G; + + // calculate the Alfven radius for an accreting neutron star + // see Equation 10 in arxiv:1912.02415 + double mDot = p_MassGainPerTimeStep / p_Stepsize; + double R_CM_6 = r_cm * r_cm * r_cm * r_cm * r_cm * r_cm; + double p = R_CM_6 * R_CM_6 / (mass_g * mass_g * mDot); + double q = PPOW(p, 1.0/7.0); + double alfvenConst = PPOW(2 * PI_2 / G_CGS, 1.0/7.0) ; + double alfvenRadius = alfvenConst * q * PPOW(initialMagField_G, 4.0/7.0); + double magneticRadius = alfvenRadius / 2.0; + // std::cout<< "MagRad: " << mDot << " " << r_cm << " " << p << " " << mass_g << " " << q << " " << alfvenConst << std::endl; + // calculate the difference in the keplerian angular velocity at the magnetic radius + // and surface angular velocity of the neutron star + // magnetic radius is half of alfven radius + // see Equation 2 in 1994MNRAS.269..455J / Equation 9 in arxiv:1912.02415 + double keplerianVelocityAtMagneticRadius = std::sqrt((G_CGS) * mass_g / magneticRadius); + double keplerianAngularVelocityAtMagneticRadius = keplerianVelocityAtMagneticRadius / magneticRadius; + double omegaDifference = keplerianAngularVelocityAtMagneticRadius - p_SpinFrequency; + // std::cout<< "Denomenators: " << p_Kappa << " " << mass_g << " " << p_Stepsize << " " << magneticRadius << std::endl; + // calculate the change in angular momentum due to accretion + // see Equation 12 in arXiv:0805.0059/ Equation 8 in arxiv:1912.02415 + double deltaJ = p_Epsilon * omegaDifference * magneticRadius * magneticRadius ; + double Jdot = p_Epsilon * omegaDifference * magneticRadius * magneticRadius * mDot; + return std::make_tuple(newPulsarMagneticField, + (angularMomentum + Jdot * p_Stepsize) / p_MoI, + deltaJ * mDot / p_MoI, + deltaJ); +} + + +/* + * Update the magnetic field and spins of neutron stars. * Modifies the following class member variables: * * m_AngularMomentum_CGS @@ -438,66 +559,134 @@ void NS::SpinDownIsolatedPulsar(const double p_Stepsize) { * @param [IN] p_CommonEnvelope Indicates whether there there is a common envelope - true or false * @param [IN] p_RecycledNS Indicates whether this star is/was a recycled neutron star - true or false * @param [IN] p_Stepsize Timestep size for integration (in seconds) - * @param [IN] p_MassGainPerTimeStep Mass loss from the secondary for each iteration (in kg) - * @param [IN] p_Epsilon Uncertainty due to mass loss - * @return Tuple containing the Maximum Mass Acceptance Rate and the Accretion Efficiency Parameter + * @param [IN] p_MassGainPerTimeStep Mass transferred from the secondary for each iteration (in kg) + * @param [IN] p_Epsilon Efficiency factor allowing for uncertainties of coupling magnetic field and matter. */ void NS::UpdateMagneticFieldAndSpin(const bool p_CommonEnvelope, const bool p_RecycledNS, const double p_Stepsize, const double p_MassGainPerTimeStep, const double p_Epsilon) { - - 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 - SpinDownIsolatedPulsar(p_Stepsize); + double kappa = OPTIONS->PulsarMagneticFieldDecayMassscale() * MSOL_TO_G; + double magFieldLowerLimit_G = PPOW(10.0, OPTIONS->PulsarLog10MinimumMagneticField()) ; + double massNS = m_Mass; + double radiusNS = m_Radius; + if ((!p_RecycledNS && !p_CommonEnvelope) || (p_RecycledNS && utils::Compare(p_MassGainPerTimeStep, 0.0) == 0 )) { + // These are the ''classical'' isolated pulsars. They experience spin-down. + SpinDownIsolatedPulsar(p_Stepsize); } - else if (utils::Compare(m_PulsarDetails.spinFrequency, _2_PI * 1000.0) < 0 && - (p_RecycledNS || p_CommonEnvelope) && utils::Compare(p_MassGainPerTimeStep, 0.0) > 0) { - - // This part of the code does pulsar recycling through accretion - // 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 * KM_TO_M; // in metres - - double MoI_SI = m_MomentOfInertia_CGS * CGS_SI; - double angularMomentum_SI = m_AngularMomentum_CGS * CGS_SI; - - double newPulsarMagneticField = (initialMagField - magFieldLowerLimit) * exp(-p_MassGainPerTimeStep / 1000.0 / kappa) + magFieldLowerLimit; - - // 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); + + else if ( utils::Compare(p_MassGainPerTimeStep, 0.0) > 0) { + // Next, we update pulsar that goes through mass transfer. + // Note that we do not need to set hard lower limit on the spin period, + // as propeller effect should be included in the calculation, + // which means pulsars will start spinning down when the AM is no longer transferred along mass transfer. + if ((!p_CommonEnvelope) || (p_CommonEnvelope && OPTIONS->NeutronStarAccretionInCE() == NS_ACCRETION_IN_CE::DISK)){ + // Pulsar evolution during stable mass transfer (RLOF), or if using the 'DISK' option + // in the CE accretion option + + //double newAM = -1.0; //updated angular momentum + double thisB = m_PulsarDetails.magneticField * TESLA_TO_GAUSS ; + double thisF = m_PulsarDetails.spinFrequency ; + double thisFdot = m_PulsarDetails.spinDownRate ; + double thisM = m_Mass ; + double divideTimestepBy = 10.0; + + //while (utils::Compare(newAM, 0.0) < 0) { + double thisMassGain = p_MassGainPerTimeStep / G_TO_KG / divideTimestepBy; + double thisTimestepSize = p_Stepsize / divideTimestepBy; + controlled_stepper_type controlled_stepper; + //boost::numeric::odeint::runge_kutta4 stepper; + + state_type x(5); + x[0] = m_AngularMomentum_CGS; + x[1] = thisB; + x[2] = thisF; + x[3] = thisFdot; + x[4] = thisM; + std::cout<< "the original AM is " << x[0] << std::endl; + std::cout<< massNS << ", "<< radiusNS << ", "<< thisB<< ", " << thisF<< ", " << x[0]<< ", " << thisTimestepSize<< ", "<< kappa<< ", " << p_Epsilon << std::endl; + + // Solve for the angular momentum of the NS after accretion. + // Use boost adaptive ODE solver for speed and accuracy and ensure the result is always positive. + struct DynamicODE + { + double mass, radius, magField, spinFrequency, angularMomentum, stepsize, kkappa, epsilon, momentI; + DynamicODE( double mass0, double radius0, double magField0, double spinFreq0, double AM0, + double stepsize0, double kkappa0, double epsilon0, double momentI0) : mass(mass0), radius(radius0), magField(magField0), + spinFrequency(spinFreq0), angularMomentum(AM0), stepsize(stepsize0), kkappa(kkappa0), epsilon(epsilon0), momentI(momentI0) {} + + void operator()( state_type& x , state_type& dxdt , double p_MassChange ) const { + DBL_DBL_DBL_DBL results = NS::deltaAngularMomentumByPulsarAccretion( + p_MassChange, mass, radius, magField, spinFrequency, angularMomentum, stepsize, kkappa, epsilon, momentI); + dxdt[0] = std::get<3>(results) ; + x[1] = std::get<0>(results); + x[2] = std::get<1>(results); + x[3] = std::get<2>(results); + x[4] = mass + p_MassChange/MSOL_TO_G; + } + }; + std::cout << thisMassGain << " " << p_MassGainPerTimeStep << std::endl; + //DynamicODE ode(massNS, radiusNS, thisB, thisF, x[0], thisTimestepSize, kappa, p_Epsilon); + integrate_adaptive(controlled_stepper, DynamicODE{x[4], radiusNS, x[1], x[2], x[0], thisTimestepSize, + kappa, p_Epsilon, m_MomentOfInertia_CGS}, x, 0.0, p_MassGainPerTimeStep / G_TO_KG, thisMassGain); + // for (double dm = 0.0; dm < p_MassGainPerTimeStep / G_TO_KG; dm += thisMassGain){ + // stepper.do_step( ode, x, dm, thisMassGain); + + // thisB = (thisB - magFieldLowerLimit_G) * exp(-thisMassGain / kappa) + magFieldLowerLimit_G;; + // thisB = thisB ; + // thisF = x[0] / m_MomentOfInertia_CGS; + // thisFdot = (x[0]-newAM)/ thisTimestepSize / m_MomentOfInertia_CGS; + // newAM = x[0]; + // std::cout << thisB << ", " << thisF << ", " << newAM < 0) { - m_PulsarDetails.magneticField = newPulsarMagneticField; - m_PulsarDetails.spinFrequency = angularMomentum_SI / MoI_SI; - m_PulsarDetails.spinDownRate = Jdot / MoI_SI; - m_AngularMomentum_CGS = angularMomentum_SI / CGS_SI; - } - else { + else if (p_CommonEnvelope && OPTIONS->NeutronStarAccretionInCE() == NS_ACCRETION_IN_CE::SURFACE) { + // Mass transfer through CE when accretion happens at the surface of the NS + double initialMagField = m_PulsarDetails.magneticField; + double initialMagField_G = initialMagField * TESLA_TO_GAUSS; + + double mass_g = m_Mass * MSOL_TO_G; // in g + double r_cm = m_Radius * RSOL_TO_KM * KM_TO_CM; // in cm + m_MomentOfInertia_CGS = CalculateMomentOfInertiaCGS() ; + double newPulsarMagneticField = (initialMagField_G - magFieldLowerLimit_G) * exp(-p_MassGainPerTimeStep / G_TO_KG / kappa) + magFieldLowerLimit_G; + double r_cm_3 = r_cm * r_cm * r_cm ; + // angular momentum of the accreted materials as they fall onto the surface of the NS. + double Jacc = m_MomentOfInertia_CGS * PPOW(G_CGS * mass_g /r_cm_3, 0.5) * p_MassGainPerTimeStep / G_TO_KG / mass_g ; + + m_AngularMomentum_CGS = m_AngularMomentum_CGS + Jacc ; + m_PulsarDetails.magneticField = newPulsarMagneticField * GAUSS_TO_TESLA; + m_PulsarDetails.spinFrequency = m_AngularMomentum_CGS / m_MomentOfInertia_CGS ; + m_PulsarDetails.spinDownRate = Jacc / p_Stepsize / m_MomentOfInertia_CGS; + } + + else if (p_CommonEnvelope && OPTIONS->NeutronStarAccretionInCE() == NS_ACCRETION_IN_CE::ZERO) { + // If mass transfer is happening through CE and chosen 'ZERO' in CE accretion mode, + // treat the pulsar as isolated and should be spun down. SpinDownIsolatedPulsar(p_Stepsize); - } + } } - else { - //In all other conditions, treat the pulsar as isolated. - SpinDownIsolatedPulsar(p_Stepsize); - } -} \ No newline at end of file +} diff --git a/src/NS.h b/src/NS.h index 9a8774bdc..67eb792cf 100755 --- a/src/NS.h +++ b/src/NS.h @@ -8,7 +8,7 @@ #include "Remnants.h" #include "BH.h" - +#include class BaseStar; class Remnants; @@ -46,8 +46,17 @@ class NS: virtual public BaseStar, public Remnants { static double CalculateRadiusOnPhase_Static(const double p_Mass) { return CalculateRadiusOnPhaseInKM_Static(p_Mass) * KM_TO_RSOL; } // Radius on phase in Rsol static double CalculateRemnantMass_Static(const double p_COCoreMass) { return 1.17 + (0.09 * p_COCoreMass); } // Hurley et al., eq 92 - - + static DBL_DBL_DBL_DBL deltaAngularMomentumByPulsarAccretion(const double p_MassGainPerTimeStep, + const double p_MagField, + const double p_Mass, + const double p_Radius, + const double p_SpinFrequency, + const double p_AngularMomentum, + const double p_Stepsize, + const double p_Kappa, + const double p_Epsilon, + const double p_MoI + ) ; protected: void Initialise() { @@ -86,12 +95,22 @@ class NS: virtual public BaseStar, public Remnants { 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; } 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); + DBL_DBL_DBL_DBL PulsarAccretion(const double p_MagField, + const double p_SpinPeriod, + const double p_AngularMomentum, + const double p_Stepsize, + const double p_MassGainPerTimeStep, + const double kappa, + const double p_Epsilon); + void UpdateMagneticFieldAndSpin(const bool p_CommonEnvelope, const bool p_RecycledNS, const double p_Stepsize, @@ -100,4 +119,4 @@ class NS: virtual public BaseStar, public Remnants { }; -#endif // __NS_h__ \ No newline at end of file +#endif // __NS_h__ diff --git a/src/Options.cpp b/src/Options.cpp index ceca2aac9..5a6d35da8 100644 --- a/src/Options.cpp +++ b/src/Options.cpp @@ -526,6 +526,9 @@ void Options::OptionValues::Initialise() { m_NeutronStarEquationOfState.type = NS_EOS::SSE; m_NeutronStarEquationOfState.typeString = NS_EOSLabel.at(m_NeutronStarEquationOfState.type); + // Neutron star accretion in common envelope + m_NSAccretionInCE.type = NS_ACCRETION_IN_CE::ZERO; + m_NSAccretionInCE.typeString = NS_ACCRETION_IN_CELabel.at(m_NSAccretionInCE.type) ; // Pulsar birth magnetic field distribution m_PulsarBirthMagneticFieldDistribution.type = PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION::ZERO; @@ -1807,7 +1810,11 @@ bool Options::AddOptions(OptionValues *p_Options, po::options_description *p_Opt po::value(&p_Options->m_NeutronStarEquationOfState.typeString)->default_value(p_Options->m_NeutronStarEquationOfState.typeString), ("Neutron star equation of state to use (" + AllowedOptionValuesFormatted("neutron-star-equation-of-state") + ", default = '" + p_Options->m_NeutronStarEquationOfState.typeString + "')").c_str() ) - + ( + "neutron-star-accretion-in-CE", + po::value(&p_Options->m_NSAccretionInCE.typeString)->default_value(p_Options->m_NSAccretionInCE.typeString), + ("Neutron star accretion in CE to use (" + AllowedOptionValuesFormatted("neutron-star-accretion-in-CE") + ", default = '" + p_Options->m_NSAccretionInCE.typeString + "')").c_str() + ) ( "OB-mass-loss", po::value(&p_Options->m_OBMassLoss.typeString)->default_value(p_Options->m_OBMassLoss.typeString), @@ -2247,6 +2254,11 @@ std::string Options::OptionValues::CheckAndSetOptions() { COMPLAIN_IF(!found, "Unknown Neutron Star Equation of State"); } + if (!DEFAULTED("neutron-star-accretion-in-ce")) { // neutron star accretion in ce + std::tie(found, m_NSAccretionInCE.type) = utils::GetMapKey(m_NSAccretionInCE.typeString, NS_ACCRETION_IN_CELabel, m_NSAccretionInCE.type); + COMPLAIN_IF(!found, "Unknown Neutron Star Accretion in CE"); + } + if (!DEFAULTED("OB-mass-loss")) { // OB (main sequence) loss prescription std::tie(found, m_OBMassLoss.type) = utils::GetMapKey(m_OBMassLoss.typeString, OB_MASS_LOSS_LABEL, m_OBMassLoss.type); COMPLAIN_IF(!found, "Unknown OB Mass Loss Prescription"); @@ -2546,6 +2558,7 @@ std::vector Options::AllowedOptionValues(const std::string p_Option case _("mode") : POPULATE_RET(EVOLUTION_MODE_LABEL); break; case _("neutrino-mass-loss-BH-formation") : POPULATE_RET(NEUTRINO_MASS_LOSS_PRESCRIPTION_LABEL); break; case _("neutron-star-equation-of-state") : POPULATE_RET(NS_EOSLabel); break; + case _("neutron-star-accretion-in-CE") : POPULATE_RET(NS_ACCRETION_IN_CELabel); break; case _("OB-mass-loss") : POPULATE_RET(OB_MASS_LOSS_LABEL); break; case _("orbital-period-distribution") : POPULATE_RET(ORBITAL_PERIOD_DISTRIBUTION_LABEL); break; case _("pulsar-birth-magnetic-field-distribution") : POPULATE_RET(PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION_LABEL); break; @@ -4624,6 +4637,7 @@ COMPAS_VARIABLE Options::OptionValue(const T_ANY_PROPERTY p_Property) const { case PROGRAM_OPTION::NOTES : value = Notes(); break; case PROGRAM_OPTION::NS_EOS : value = static_cast(NeutronStarEquationOfState()); break; + case PROGRAM_OPTION::NS_ACCRETION_IN_CE : value = static_cast(NeutronStarAccretionInCE()); break; case PROGRAM_OPTION::ORBITAL_PERIOD : value = OrbitalPeriod(); break; case PROGRAM_OPTION::ORBITAL_PERIOD_DISTRIBUTION : value = static_cast(OrbitalPeriodDistribution()); break; diff --git a/src/Options.h b/src/Options.h index 040b63952..550170da8 100755 --- a/src/Options.h +++ b/src/Options.h @@ -380,6 +380,8 @@ class Options { "maximum-mass-donor-nandez-ivanova", "minimum-secondary-mass", + "neutron-star-accretion-in-ce", + "orbital-period", "orbital-period-distribution", "orbital-period-max", @@ -510,6 +512,7 @@ class Options { "notes-hdrs", "neutrino-mass-loss-BH-formation", "neutron-star-equation-of-state", + "neutron-star-accretion-in-CE", "OB-mass-loss", "orbital-period-distribution", @@ -965,6 +968,8 @@ class Options { // Neutron star equation of state ENUM_OPT m_NeutronStarEquationOfState; // NS EOS + // Neutron star accretion in common envelope + ENUM_OPT m_NSAccretionInCE; // NS Accretion In CE // Pulsar birth magnetic field distribution string ENUM_OPT m_PulsarBirthMagneticFieldDistribution; // Birth magnetic field distribution for pulsars @@ -1425,6 +1430,7 @@ class Options { double NeutrinoMassLossValueBH() const { return OPT_VALUE("neutrino-mass-loss-BH-formation-value", m_NeutrinoMassLossValueBH, true); } NS_EOS NeutronStarEquationOfState() const { return OPT_VALUE("neutron-star-equation-of-state", m_NeutronStarEquationOfState.type, true); } + NS_ACCRETION_IN_CE NeutronStarAccretionInCE() const { return OPT_VALUE("neutron-star-accretion-in-ce", m_NSAccretionInCE.type, true); } std::string Notes(const size_t p_Idx) const { return OPT_VALUE("notes", m_Notes[p_Idx], true); } std::vector Notes() const { return OPT_VALUE("notes", m_Notes, true); } diff --git a/src/changelog.h b/src/changelog.h index 0e6f58346..147b95c3d 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1225,8 +1225,16 @@ // - Change TPAGB::IsSupernova() so that stars with base of AGB core masses below MCBUR1 remain on the TPAGB until they make WDs; remove ResolveTypeIIaSN() functionality. // - Add --evolve-main-sequence-mergers option which allows for main sequence merger products to continue evolution // - Update HG::CalculateRadialExtentConvectiveEnvelope() to use a combination of Hurley & Picker to avoid excessively high convective envelope densities +// 02.51.00 YS - Jul 03, 2024 - Update to neutron star accretion treatments: +// - Fixes to MSP formation/NS in mass transfer treatments: +// 1). Created a new function NS::PulsarAccretion() to calculate the pulsar evolution in stable mass transfer. +// 2). In UpdateMagneticFieldAndSpin(), splitting stable mass transfer into smaller steps so that no negative spin period is present. +// 3). Adding a new programming option "NS-ACCRETION-IN-CE" for different treatment of how neutron star would behave when in CE. + + +const std::string VERSION_STRING = "02.51.00"; + -const std::string VERSION_STRING = "02.50.00"; # endif // __changelog_h__ diff --git a/src/compasConfigDefault.yaml b/src/compasConfigDefault.yaml index 2264d93cc..c17a4272d 100644 --- a/src/compasConfigDefault.yaml +++ b/src/compasConfigDefault.yaml @@ -260,6 +260,7 @@ stringChoices: # --fryer-supernova-engine: 'DELAYED' # Default: 'DELAYED' # Options: ['DELAYED','RAPID'] # --kick-magnitude-distribution: 'MULLERMANDEL' # Default: 'MULLERMANDEL' # Options: ['MULLERMANDEL','MULLER2016MAXWELLIAN','MULLER2016','BRAYELDRIDGE','MAXWELLIAN','FLAT','FIXED','ZERO'] # --kick-direction: 'ISOTROPIC' # Default: 'ISOTROPIC' # Options: ['POLES','WEDGE','POWERLAW','PERPENDICULAR','INPLANE','ISOTROPIC'] +# --neutron-star-accretion-in-ce: 'ZERO' # Default: 'ZERO' # Options: ['ZERO', 'DISK', 'SURFACE'] # --neutron-star-equation-of-state: 'SSE' # Default: 'SSE' # Options: ['ARP3','SSE'] # --neutrino-mass-loss-BH-formation: 'FIXED_MASS' # Default: 'FIXED_MASS' # Options: ['FIXED_MASS','FIXED_FRACTION'] # --pulsar-birth-magnetic-field-distribution: 'ZERO' # Default: 'ZERO' # Options: ['LOGNORMAL','UNIFORM','FLATINLOG','FIXED','ZERO'] diff --git a/src/constants.h b/src/constants.h index d19a1c2ce..f87b584e6 100755 --- a/src/constants.h +++ b/src/constants.h @@ -685,6 +685,7 @@ enum class ERROR: int { UNKNOWN_MT_THERMALLY_LIMITED_VARIATION, // unknown mass transfer thermally limited variation UNKNOWN_NEUTRINO_MASS_LOSS_PRESCRIPTION, // unknown neutrino mass loss prescription UNKNOWN_NS_EOS, // unknown NS equation-of-state + UNKNOWN_NS_ACCRETION_IN_CE, // unknown NS accretion-in-ce UNKNOWN_PPI_PRESCRIPTION, // unknown pulsational pair instability prescription UNKNOWN_PROGRAM_OPTION, // unknown program option UNKNOWN_PROPERTY_TYPE, // unknown property type @@ -841,6 +842,7 @@ const COMPASUnorderedMap> ERROR_CATA { ERROR::UNKNOWN_MASS_LOSS_PRESCRIPTION, { ERROR_SCOPE::ALWAYS, "Unknown mass loss prescription" }}, { ERROR::UNKNOWN_NEUTRINO_MASS_LOSS_PRESCRIPTION, { ERROR_SCOPE::ALWAYS, "Unknown neutrino mass loss prescription" }}, { ERROR::UNKNOWN_NS_EOS, { ERROR_SCOPE::ALWAYS, "Unknown NS equation-of-state" }}, + { ERROR::UNKNOWN_NS_ACCRETION_IN_CE, { ERROR_SCOPE::ALWAYS, "Unknown NS accretion-in-ce" }}, { ERROR::UNKNOWN_PPI_PRESCRIPTION, { ERROR_SCOPE::ALWAYS, "Unknown pulsational pair instability prescription" }}, { ERROR::UNKNOWN_PROGRAM_OPTION, { ERROR_SCOPE::ALWAYS, "Unknown program option - property details not found" }}, { ERROR::UNKNOWN_PROPERTY_TYPE, { ERROR_SCOPE::ALWAYS, "Unknown property type - property details not found" }}, @@ -1282,6 +1284,14 @@ const COMPASUnorderedMap NS_EOSLabel = { { NS_EOS::ARP3, "ARP3" } }; +// Neutron Star accretion in CE +enum class NS_ACCRETION_IN_CE: int { ZERO, DISK, SURFACE }; +const COMPASUnorderedMap NS_ACCRETION_IN_CELabel = { + { NS_ACCRETION_IN_CE::ZERO, "ZERO" }, + { NS_ACCRETION_IN_CE::DISK, "DISK" }, + { NS_ACCRETION_IN_CE::SURFACE, "SURFACE" } +}; + // Orbital Period Distributions enum class ORBITAL_PERIOD_DISTRIBUTION: int { FLATINLOG }; @@ -2587,6 +2597,7 @@ enum class PROGRAM_OPTION: int { NOTES, NS_EOS, + NS_ACCRETION_IN_CE, ORBITAL_PERIOD, ORBITAL_PERIOD_DISTRIBUTION, diff --git a/src/typedefs.h b/src/typedefs.h index e52b93a7f..50299dbe7 100755 --- a/src/typedefs.h +++ b/src/typedefs.h @@ -312,5 +312,6 @@ typedef struct StellarCEDetails { // Common Envelope detail typedef std::vector< double > state_type; typedef boost::numeric::odeint::runge_kutta_cash_karp54< state_type > error_stepper_type; typedef boost::numeric::odeint::controlled_runge_kutta< error_stepper_type > controlled_stepper_type; +typedef boost::numeric::odeint::runge_kutta4< state_type > stepper; #endif // __typedefs_h__ diff --git a/src/yaml.h b/src/yaml.h index 52d854a5b..8b1793c43 100644 --- a/src/yaml.h +++ b/src/yaml.h @@ -329,6 +329,7 @@ namespace yaml { " --kick-magnitude-distribution", " --kick-direction", " --neutron-star-equation-of-state", + " --neutron-star-accretion-in-ce", " --neutrino-mass-loss-BH-formation", " --pulsar-birth-magnetic-field-distribution", " --pulsar-birth-spin-period-distribution",