Merge pull request #520 from TeamCOMPAS/lambda_enhancements

Lambda enhancements

Author: ilyamandel

online-docs/pages/User guide/Program options/program-options-list-defaults.rst

@@ -137,6 +137,22 @@ Default = 0.1
 Multiplicative constant to be applied to the common envelope lambda parameter for any prescription. |br|
 Default = 1.0
 
+**--common-envelope-lambda-nanjing-enhanced** |br|
+Continuous extrapolation beyond maximum radius range in Nanjing lambda's as implemented in StarTrack. Only used when ``--common-envelope-lambda-prescription = LAMBDA_NANJING``. |br|
+Default = FALSE
+
+**--common-envelope-lambda-nanjing-interpolate-in-mass** |br|
+Interpolate Nanjing lambda parameters across different mass models. Only used when ``--common-envelope-lambda-prescription = LAMBDA_NANJING``. |br|
+Default = FALSE
+
+**--common-envelope-lambda-nanjing-interpolate-in-metallicity** |br|
+Interpolate Nanjing lambda parameters across population I and population II metallicity models. Only used when ``--common-envelope-lambda-prescription = LAMBDA_NANJING``. |br|
+Default = FALSE
+
+**--common-envelope-lambda-nanjing-use_rejuvenated-mass** |br|
+Use rejuvenated or effective ZAMS mass instead of true birth mass when computing Nanjing lambda parameters. Only used when ``--common-envelope-lambda-prescription = LAMBDA_NANJING``. |br|
+Default = FALSE
+
 **--common-envelope-lambda-prescription** |br|
 CE lambda (envelope binding energy) prescription. |br|
 Options: { LAMBDA_FIXED, LAMBDA_LOVERIDGE, LAMBDA_NANJING, LAMBDA_KRUCKOW, LAMBDA_DEWI } |br|
@@ -955,7 +971,7 @@ Go to :ref:`the top of this page <options-props-top>` for the full alphabetical
 
 --circulariseBinaryDuringMassTransfer, --angular-momentum-conservation-during-circularisation
 
---envelope-state-prescription, --common-envelope-alpha, --common-envelope-alpha-thermal, --common-envelope-lambda-prescription, --common-envelope-lambda, --common-envelope-slope-kruckow, --common-envelope-lambda-multiplier, --common-envelope-allow-main-sequence-survive, --common-envelope-allow-radiative-envelope-survive*, --common-envelope-allow-immediate-RLOF-post-CE-survive, --common-envelope-mass-accretion-prescription, --common-envelope-mass-accretion-constant, --common-envelope-mass-accretion-min, --common-envelope-mass-accretion-max, --common-envelope-recombination-energy-density, --maximum-mass-donor-nandez-ivanova, --revised-energy-formalism-nandez-ivanova
+--envelope-state-prescription, --common-envelope-alpha, --common-envelope-alpha-thermal, --common-envelope-lambda-prescription, --common-envelope-lambda, --common-envelope-slope-kruckow, --common-envelope-lambda-multiplier, --common-envelope-lambda-nanjing-enhanced, --common-envelope-lambda-nanjing-interpolate-in-mass, --common-envelope-lambda-nanjing-interpolate-in-metallicity, --common-envelope-lambda-nanjing-use_rejuvenated-mass, --common-envelope-allow-main-sequence-survive, --common-envelope-allow-radiative-envelope-survive*, --common-envelope-allow-immediate-RLOF-post-CE-survive, --common-envelope-mass-accretion-prescription, --common-envelope-mass-accretion-constant, --common-envelope-mass-accretion-min, --common-envelope-mass-accretion-max, --common-envelope-recombination-energy-density, --maximum-mass-donor-nandez-ivanova, --revised-energy-formalism-nandez-ivanova
 
 :ref:`Back to Top <options-props-top>`
 

src/BaseBinaryStar.cpp

@@ -223,7 +223,7 @@ BaseStar::BaseStar(const unsigned long int p_RandomSeed,
 
 
 /*
- * Determine the value of the requested property of the constitusnt star (parameter p_Property)
+ * Determine the value of the requested property of the constituent star (parameter p_Property)
  *
  * The property is a boost variant variable, and is one of the following types:
  *
@@ -1084,6 +1084,161 @@ double BaseStar::CalculateLambdaLoveridgeEnergyFormalism(const double p_EnvMass,
 }
 
 
+/* 
+ * Wrapper function to return Nanjing lambda based on options
+ * 
+ * 
+ * double BaseStar::CalculateLambdaNanjing()
+ * 
+ * @return                                      Common envelope lambda parameter
+ */ 
+double BaseStar::CalculateLambdaNanjing() const {
+
+    double mass = m_MZAMS;
+    double lambda = 0.0;
+    if (OPTIONS->CommonEnvelopeLambdaNanjingUseRejuvenatedMass()) {mass = m_Mass0;}                              // Use rejuvenated mass to calculate lambda instead of true birth mass
+    
+    if (OPTIONS->CommonEnvelopeLambdaNanjingEnhanced()) {                                                        // If using enhanced Nanjing lambda's
+        if (OPTIONS->CommonEnvelopeLambdaNanjingInterpolateInMass()) {
+            if (OPTIONS->CommonEnvelopeLambdaNanjingInterpolateInMetallicity()) {
+                lambda = BaseStar::CalculateMassAndZInterpolatedLambdaNanjing(mass, m_Metallicity);
+            }
+            else {
+                int Zind = 0;
+                if (utils::Compare(m_Metallicity, LAMBDA_NANJING_ZLIMIT) < 0) {Zind = 0;} else {Zind = 1;}
+                lambda = BaseStar::CalculateMassInterpolatedLambdaNanjing(mass, Zind);
+            }
+        }
+        else {
+            int massInd = BaseStar::FindLambdaNanjingNearestMassIndex(mass);                                    // Do not interpolate in mass, so need to use nearest mass bin
+            if (OPTIONS->CommonEnvelopeLambdaNanjingInterpolateInMetallicity()) {
+                lambda = BaseStar::CalculateZInterpolatedLambdaNanjing(m_Metallicity, massInd);
+            }
+            else {
+                int Zind = 0;
+                if (utils::Compare(m_Metallicity, LAMBDA_NANJING_ZLIMIT) < 0) {Zind = 0;} else {Zind = 1;}
+                lambda = BaseStar::CalculateLambdaNanjingEnhanced(massInd, Zind);
+            }
+        }
+    }
+    else { lambda = CalculateLambdaNanjingStarTrack(mass, m_Metallicity); }
+    return lambda;
+}
+
+
+/* 
+ * Calculate mass- and metallicity-interpolated Nanjing lambda
+ * 
+ * 
+ * double BaseStar::CalculateMassAndZInterpolatedLambdaNanjing(const double p_Mass, const double p_Z)
+ * 
+ * @return                                      Common envelope lambda parameter
+ */ 
+double BaseStar::CalculateMassAndZInterpolatedLambdaNanjing(const double p_Mass, const double p_Z) const {
+
+    double lambda = 0.0;
+    if (utils::Compare(m_Metallicity, LAMBDA_NANJING_POPII_Z) < 0) {
+        lambda = BaseStar::CalculateMassInterpolatedLambdaNanjing(p_Mass, 0);                    // Use lambda for pop. II metallicity
+    }
+    else if (utils::Compare(m_Metallicity, LAMBDA_NANJING_POPI_Z) > 0) {
+        lambda = BaseStar::CalculateMassInterpolatedLambdaNanjing(p_Mass, 1);                    // Use lambda for pop. I metallicity
+    }
+    else {                                                                                       // Linear interpolation in logZ between pop. I and pop. II metallicities
+        const double logZ = log(m_Metallicity);
+        double lambdaLow = BaseStar::CalculateMassInterpolatedLambdaNanjing(p_Mass, 0);
+        double lambdaUp  = BaseStar::CalculateMassInterpolatedLambdaNanjing(p_Mass, 1);
+        lambda = lambdaLow + (logZ - LAMBDA_NANJING_POPII_LOGZ) / (LAMBDA_NANJING_POPI_LOGZ - LAMBDA_NANJING_POPII_LOGZ) * (lambdaUp - lambdaLow);
+    }
+    return lambda;
+}
+
+
+/* 
+ * Interpolate Nanjing lambda in mass for a given metallicity
+ * 
+ * 
+ * double BaseStar::CalculateMassInterpolatedLambdaNanjing(const int p_Zind)
+ * 
+ * @param   [IN]    p_Zind                      Index of metallicity (0 = pop II metallicity, 1 = pop I metallicity)
+ * @param   [IN]    p_Mass                      Mass / Msun to evaluate lambda with
+ * @return                                      Common envelope lambda parameter
+ */ 
+double BaseStar::CalculateMassInterpolatedLambdaNanjing(const double p_Mass, const int p_Zind) const {
+
+    double lambda = 0.0;
+    std::vector<int> ind = utils::binarySearch(NANJING_MASSES, p_Mass);
+    int low = ind[0];
+    int up = ind[1];
+    if ( (low < 0)  && (up >= 0) ) {                                                            // Mass below range calculated by Xu & Li (2010)
+        lambda = CalculateLambdaNanjingEnhanced(0, p_Zind);                                     // Use lambda for minimum mass
+    }
+    else if ( (low >= 0) && (up < 0) ) {                                                        // Mass above range calculated by Xu & Li (2010)
+        lambda = CalculateLambdaNanjingEnhanced(NANJING_MASSES.size() - 1, p_Zind);             // Use lambda for maximum mass
+    }
+    else if (low == up) {                                                                       // Mass is exactly equal to the mass of a model evolved by Xu & Li (2010)
+        lambda = CalculateLambdaNanjingEnhanced(low, p_Zind);
+    }
+    else {                                                                                      // Linear interpolation between upper and lower mass bins
+        double lambdaLow = CalculateLambdaNanjingEnhanced(low, p_Zind);
+        double lambdaUp  = CalculateLambdaNanjingEnhanced(up, p_Zind);
+        lambda = lambdaLow + (p_Mass - NANJING_MASSES[low]) / (NANJING_MASSES[up] - NANJING_MASSES[low]) * (lambdaUp - lambdaLow);
+    }
+    return lambda;
+}
+
+
+/* 
+ * Interpolate Nanjing lambda in metallicity for a given mass
+ * 
+ * 
+ * double BaseStar::CalculateZInterpolatedLambdaNanjing(const int p_Zind)
+ * 
+ * @param   [IN]    p_Z                         Metallicity
+ * @param   [IN]    p_MassInd                   Index specifying donor mass (see NANJING_MASSES in constants.h)
+ * @return                                      Common envelope lambda parameter
+ */ 
+double BaseStar::CalculateZInterpolatedLambdaNanjing(const double p_Z, const int p_MassInd) const {
+
+    double lambda = 0.0;
+    if (utils::Compare(m_Metallicity, LAMBDA_NANJING_POPII_Z) < 0) {
+        lambda = CalculateLambdaNanjingEnhanced(p_MassInd, 0);                       // Use lambda for pop. II metallicity
+    }
+    else if (utils::Compare(m_Metallicity, LAMBDA_NANJING_POPI_Z) > 0) {
+        lambda = CalculateLambdaNanjingEnhanced(p_MassInd, 1);                       // Use lambda for pop. I metallicity
+    }
+    else {                                                                           // Linear interpolation in logZ between pop. I and pop. II metallicities
+        const double logZ = log(m_Metallicity);
+        double lambdaLow = CalculateLambdaNanjingEnhanced(p_MassInd, 0);
+        double lambdaUp  = CalculateLambdaNanjingEnhanced(p_MassInd, 1);
+        lambda = lambdaLow + (logZ - LAMBDA_NANJING_POPII_LOGZ) / (LAMBDA_NANJING_POPI_LOGZ - LAMBDA_NANJING_POPII_LOGZ) * (lambdaUp - lambdaLow);
+    }
+    return lambda;
+}
+
+
+/* 
+ * Returns index in NANJING_MASSES corresponding to nearest mass model computed by Xu & Li (2010)
+ * 
+ * 
+ * double BaseStar::FindLambdaNanjingNearestMassIndex(const double p_Mass)
+ * 
+ * @return                                      Index in NANJING_MASSES
+ */ 
+double BaseStar::FindLambdaNanjingNearestMassIndex(const double p_Mass) const {
+
+    if (p_Mass < NANJING_MASSES_MIDPOINTS[0]) {                                                                  // M < 1.5 Msun, use lambda for the 1 Msun model
+        return 0;
+    }
+    else if (p_Mass >= NANJING_MASSES_MIDPOINTS.back()) {                                                        // M >= 75 Msun, use lambda for the 100 Msun model
+        return NANJING_MASSES.size() - 1;
+    }
+    else {                                                                                                       // Search for upper and lower mass bin edges
+        std::vector<int> ind = utils::binarySearch(NANJING_MASSES_MIDPOINTS, p_Mass);
+        return ind[1];
+    }
+}
+
+
 ///////////////////////////////////////////////////////////////////////////////////////
 //                                                                                   //
 //                                 ZETA CALCULATIONS                                 //

src/BaseBinaryStar.h

@@ -386,7 +386,8 @@ class BaseStar {
     virtual double              CalculateLambdaDewi() const                                                             { SHOW_WARN(ERROR::NO_LAMBDA_DEWI, "Default used: 1.0"); return 1.0; }      // Not supported: show error
             double              CalculateLambdaKruckow(const double p_Radius, const double p_Alpha) const;
             double              CalculateLambdaLoveridgeEnergyFormalism(const double p_EnvMass, const double p_IsMassLoss = false) const;
-    virtual double              CalculateLambdaNanjing() const                                                          { SHOW_WARN(ERROR::NO_LAMBDA_NANJING, "Default used: 1.0"); return 1.0; }   // Not supported: show error
+    virtual double              CalculateLambdaNanjingStarTrack(const double p_Mass, const double p_Metallicity) const           { SHOW_WARN(ERROR::NO_LAMBDA_NANJING, "Default used: 1.0"); return 1.0; }   // Not supported: show error
+    virtual double              CalculateLambdaNanjingEnhanced(const int p_MassInd, const int p_Zind) const             { SHOW_WARN(ERROR::NO_LAMBDA_NANJING, "Default used: 1.0"); return 1.0; }   // Not supported: show error
 
             void                CalculateLCoefficients(const double p_LogMetallicityXi, DBL_VECTOR &p_LCoefficients);
 
@@ -401,6 +402,10 @@ class BaseStar {
             double              CalculateLuminosityGivenCoreMass(const double p_CoreMass) const;
     virtual double              CalculateLuminosityOnPhase() const                                                      { return m_Luminosity; }                                                    // Default is NO-OP
 
+            double              CalculateMassAndZInterpolatedLambdaNanjing(const double p_Mass, const double p_Z) const;
+            double              CalculateMassInterpolatedLambdaNanjing(const double p_Mass, const int p_Zind) const;
+            double              CalculateZInterpolatedLambdaNanjing(const double p_Z, const int p_MassInd) const;
+
             void                CalculateMassCutoffs(const double p_Metallicity, const double p_LogMetallicityXi, DBL_VECTOR &p_MassCutoffs);
 
     static  double              CalculateMassLoss_Static(const double p_Mass, const double p_Mdot, const double p_Dt);
@@ -498,12 +503,14 @@ class BaseStar {
                                                     const double p_EjectaMass,
                                                     const double p_RemnantMass) const;
 
+            double              CalculateLambdaNanjing() const;
     virtual void                EvolveOneTimestepPreamble() { };                                                                                                                                    // Default is NO-OP
 
             STELLAR_TYPE        EvolveOnPhase();
 
     virtual STELLAR_TYPE        EvolveToNextPhase()                                                                     { return m_StellarType; }
 
+            double              FindLambdaNanjingNearestMassIndex(const double p_Mass) const;
 
     virtual bool                IsEndOfPhase() const                                                                    { return false; }
     virtual bool                IsSupernova() const                                                                     { return false; }