Flexible initial stellar type

online-docs/pages/User guide/COMPAS output/standard-logfiles-default-record-specifications-BSE-detailed.rst

@@ -10,12 +10,12 @@ Default record definition for the BSE Detailed Output log file::
         BINARY_PROPERTY::UNBOUND,
         BINARY_PROPERTY::SEMI_MAJOR_AXIS_RSOL,
         BINARY_PROPERTY::ECCENTRICITY,
-        STAR_1_PROPERTY::MZAMS,
-        STAR_2_PROPERTY::MZAMS,
-        STAR_1_PROPERTY::MASS_0,
-        STAR_2_PROPERTY::MASS_0,
         STAR_1_PROPERTY::MASS,
         STAR_2_PROPERTY::MASS,
+        STAR_1_PROPERTY::MASS_0,
+        STAR_2_PROPERTY::MASS_0,
+        STAR_1_PROPERTY::MZAMS,
+        STAR_2_PROPERTY::MZAMS,
         STAR_1_PROPERTY::ENV_MASS,
         STAR_2_PROPERTY::ENV_MASS,
         STAR_1_PROPERTY::CORE_MASS,

online-docs/pages/User guide/COMPAS output/standard-logfiles-default-record-specifications-BSE-sysparms.rst

@@ -36,4 +36,3 @@ Default record definition for the BSE System Parameters log file::
         BINARY_PROPERTY::ERROR,
         PROGRAM_OPTION::NOTES
     };
-

online-docs/pages/User guide/COMPAS output/standard-logfiles-record-specification-stellar.rst

@@ -872,6 +872,22 @@ Following is an alphabetical list of stellar properties available for inclusion
 `Note that this property has the same header string as BINARY_PROPERTY::ID & BINARY_PROPERTY::RLOF_CURRENT_ID. It is expected that one or 
 the other is printed in any file, but not both. If both are printed then the file will contain two columns with the same header string.`
 
+.. flat-table::
+   :widths: 25 75 1 1
+   :header-rows: 0
+   :class: aligned-text
+
+   * - :cspan:`2` **INITIAL_RADIUS**
+     -
+   * - Data type:
+     - DOUBLE
+   * - COMPAS variable:
+     - BaseStar::m_InitialRadius
+   * - Description:
+     - Initial radius of the star at the beginning of the simulation
+   * - Header Strings:
+     - Initial_Radius(1), Initial_Radius(2)
+
 .. flat-table::
    :widths: 25 75 1 1
    :header-rows: 0

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

@@ -535,6 +535,24 @@ Default = 5.0
 Single power law power to generate primary mass using ``POWERLAW`` IMF. |br|
 Default = 0.0
 
+**--initial-stellar-type** |br|
+Initial stellar type for a single star when evolving in SSE mode. |br|
+Options: { MS, HeMS, HeWD, COWD, ONeWD, NS, BH } |br|
+Corresponding to Main Sequence, Helium Main Sequence, Helium White Dwarf, Carbon-Oxygen White Dwarf, Oxygen-Neon White Dwarf, Neutron Star, and Black Hole.
+Default = MS
+
+**--initial-stellar-type-1** |br|
+Initial stellar type for the primary star when evolving in BSE mode |br|
+Options: { MS, HeMS, HeWD, COWD, ONeWD, NS, BH } |br|
+Corresponding to Main Sequence, Helium Main Sequence, Helium White Dwarf, Carbon-Oxygen White Dwarf, Oxygen-Neon White Dwarf, Neutron Star, and Black Hole.
+Default = MS
+
+**--initial-stellar-type-2** |br|
+Initial stellar type for the secondary star when evolving in BSE mode |br|
+Options: { MS, HeMS, HeWD, COWD, ONeWD, NS, BH } |br|
+Corresponding to Main Sequence, Helium Main Sequence, Helium White Dwarf, Carbon-Oxygen White Dwarf, Oxygen-Neon White Dwarf, Neutron Star, and Black Hole.
+Default = MS
+
 .. _options-props-J:
 
 .. _options-props-K:
@@ -1464,6 +1482,8 @@ Go to :ref:`the top of this page <options-props-top>` for the full alphabetical
 
 --initial-mass-function, --initial-mass, --initial-mass-1, --initial-mass-2, --initial-mass-min, --initial-mass-max, --initial-mass-power
 
+--initial-stellar-type, --initial-stellar-type-1, --initial-stellar-type-2
+
 --mass-ratio-distribution, --mass-ratio, --mass-ratio-min, --mass-ratio-max, --minimum-secondary-mass
 
 --eccentricity-distribution, --eccentricity, --eccentricity-min, --eccentricity-max

src/BH.h

@@ -55,9 +55,20 @@ class BH: virtual public BaseStar, public Remnants {
         m_Mass0      = 0.0;
 
         EvolveOnPhase(0.0);
-    }
-
+   }
+
+    void FastForward() {                                                                                                                                                        // Set stellar attributes for stars initialized to this stellar type
+
+        m_Radius                                   = CalculateRadiusOnPhase();
+        m_Luminosity                               = CalculateLuminosityOnPhase();
 
+        m_InitialLuminosity                        = m_Luminosity;
+        m_InitialRadius                            = m_Radius;
+        m_InitialStellarType                       = m_StellarType;
+        m_StellarTypePrev                          = m_StellarType;
+    }
+
+
     // member functions - alphabetically
     double  CalculateConvergedMassStepZetaNuclear() const                           { return 0.0; }
     double  CalculateEddingtonCriticalRate() const                                  { return 2.6E-8 * m_Mass * MYR_TO_YEAR; }                           // E.g., Marchant+, 2017, Eq. 3, assuming accretion efficiency of 10%

src/BaseBinaryStar.cpp

@@ -110,6 +110,9 @@ BaseBinaryStar::BaseBinaryStar(const unsigned long int p_Seed, const long int p_
 
             mass2 = mass1 * q;                                                                                                          // calculate mass2 using mass ratio                                                                     
         }
+
+        STELLAR_TYPE initialStellarType1 = OPTIONS->InitialStellarType1();                                                              // Get primary stellar type, possibly user specified
+        STELLAR_TYPE initialStellarType2 = OPTIONS->InitialStellarType2();                                                              // Get secondary stellar type, possibly user specified
 
         double metallicity = !OPTIONS->OptionDefaulted("metallicity")                                                                   // user specified metallicity?
                                 ? OPTIONS->Metallicity()                                                                                // yes, use it
@@ -158,12 +161,12 @@ BaseBinaryStar::BaseBinaryStar(const unsigned long int p_Seed, const long int p_
         // binary star contains two instances of star to hold masses, radii and luminosities.
         // star 1 initially more massive
         m_Star1 = !OPTIONS->OptionDefaulted("rotational-frequency-1")                                                                   // user specified primary rotational frequency?
-                    ? 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())
-                    : new BinaryConstituentStar(m_RandomSeed, mass1, metallicity, kickParameters1);                                     // no - let it be calculated
+                    ? new BinaryConstituentStar(m_RandomSeed, mass1, initialStellarType1, metallicity, kickParameters1, OPTIONS->RotationalFrequency1() * SECONDS_IN_YEAR) // yes - use it (convert from Hz to cycles per year - see BaseStar::CalculateZAMSAngularFrequency())
+                    : new BinaryConstituentStar(m_RandomSeed, mass1, initialStellarType1, metallicity, kickParameters1);                                                   // no - let it be calculated
 
         m_Star2 = !OPTIONS->OptionDefaulted("rotational-frequency-2")                                                                   // user specified secondary rotational frequency?
-                    ? 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())
-                    : new BinaryConstituentStar(m_RandomSeed, mass2, metallicity, kickParameters2);                                     // no - let it be calculated
+                    ? new BinaryConstituentStar(m_RandomSeed, mass2, initialStellarType2, metallicity, kickParameters2, OPTIONS->RotationalFrequency2() * SECONDS_IN_YEAR) // yes - use it (convert from Hz to cycles per year - see BaseStar::CalculateZAMSAngularFrequency())
+                    : new BinaryConstituentStar(m_RandomSeed, mass2, initialStellarType2, metallicity, kickParameters2);                                                   // no - let it be calculated
 
         double starToRocheLobeRadiusRatio1 = (m_Star1->Radius() * RSOL_TO_AU) / (m_SemiMajorAxis * (1.0 - m_Eccentricity) * CalculateRocheLobeRadius_Static(mass1, mass2));
         double starToRocheLobeRadiusRatio2 = (m_Star2->Radius() * RSOL_TO_AU) / (m_SemiMajorAxis * (1.0 - m_Eccentricity) * CalculateRocheLobeRadius_Static(mass2, mass1));
@@ -188,13 +191,13 @@ BaseBinaryStar::BaseBinaryStar(const unsigned long int p_Seed, const long int p_
             // create new stars with equal masses - all other ZAMS values recalculated
             delete m_Star1;
             m_Star1 = !OPTIONS->OptionDefaulted("rotational-frequency-1")                                                               // user specified primary rotational frequency?
-                        ? 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())
-                        : new BinaryConstituentStar(m_RandomSeed, mass1, metallicity, kickParameters1);                                 // no - let it be calculated
+                        ? new BinaryConstituentStar(m_RandomSeed, mass1, initialStellarType1, metallicity, kickParameters1, OPTIONS->RotationalFrequency1() * SECONDS_IN_YEAR) // yes - use it (convert from Hz to cycles per year - see BaseStar::CalculateZAMSAngularFrequency())
+                        : new BinaryConstituentStar(m_RandomSeed, mass1, initialStellarType1, metallicity, kickParameters1);                                 // no - let it be calculated
 
             delete m_Star2;
             m_Star2 = !OPTIONS->OptionDefaulted("rotational-frequency-2")                                                               // user specified secondary rotational frequency?
-                        ? 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())
-                        : new BinaryConstituentStar(m_RandomSeed, mass2, metallicity, kickParameters2);                                 // no - let it be calculated
+                        ? new BinaryConstituentStar(m_RandomSeed, mass2, initialStellarType2, metallicity, kickParameters2, OPTIONS->RotationalFrequency2() * SECONDS_IN_YEAR) // yes - use it (convert from Hz to cycles per year - see BaseStar::CalculateZAMSAngularFrequency())
+                        : new BinaryConstituentStar(m_RandomSeed, mass2, initialStellarType2, metallicity, kickParameters2);                                 // no - let it be calculated
 
             starToRocheLobeRadiusRatio1 = (m_Star1->Radius() * RSOL_TO_AU) / (m_SemiMajorAxis * CalculateRocheLobeRadius_Static(mass1, mass2)); //eccentricity already zero
             starToRocheLobeRadiusRatio2 = (m_Star2->Radius() * RSOL_TO_AU) / (m_SemiMajorAxis * CalculateRocheLobeRadius_Static(mass2, mass1));
@@ -301,33 +304,42 @@ void BaseBinaryStar::SetRemainingValues() {
 
         // check for CHE
         //
+        // because we've changed the rotational frequency of the constituent stars we
+        // have to reset the stellar type - at this stage, based on their rotational
+        // frequency at birth, they may have already been assigned one of MS_LTE_07,
+        // MS_GT_07, or CHEMICALLY_HOMOGENEOUS (ignore if assigned a more evolved type)
+        //
         // here we need to change from MS_* -> CH, or from CH->MS* based on the binary orbital frequency, assuming that the stars are tidally locked
         // set the spin to the orbital frequency, unless the user has specified a spin frequency
 
         // star 1
-        if (utils::Compare(omega, m_Star1->OmegaCHE()) >= 0) {                                                                                          // star 1 CH?
-            if (m_Star1->StellarType() != STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) {
-                (void)m_Star1->SwitchTo(STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS, true);                                                                    // yes, switch if not already Chemically Homogeneous
-                if (OPTIONS->OptionDefaulted("rotational-frequency-1")) m_Star1->SetOmega(omega); }                                                     // set spin to orbital frequency unless user specified
-        }
-        else if (m_Star1->MZAMS() <= 0.7) {                                                                                                             // no - MS - initial mass determines actual type  (don't use utils::Compare() here)
-            if (m_Star1->StellarType() != STELLAR_TYPE::MS_LTE_07) (void)m_Star1->SwitchTo(STELLAR_TYPE::MS_LTE_07, true);                              // MS <= 0.7 Msol - switch if necessary
-        }
-        else {
-            if (m_Star1->StellarType() != STELLAR_TYPE::MS_GT_07) (void)m_Star1->SwitchTo(STELLAR_TYPE::MS_GT_07, true);                                // MS > 0.7 Msol - switch if necessary
+        if (m_Star1->IsOneOf({ STELLAR_TYPE::MS_LTE_07, STELLAR_TYPE::MS_GT_07, STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS })) {
+            if (utils::Compare(omega, m_Star1->OmegaCHE()) >= 0) {                                                                                          // star 1 CH?
+                if (m_Star1->StellarType() != STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) {
+                    (void)m_Star1->SwitchTo(STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS, true);                                                                    // yes, switch if not already Chemically Homogeneous
+                    if (OPTIONS->OptionDefaulted("rotational-frequency-1")) m_Star1->SetOmega(omega); }                                                     // set spin to orbital frequency unless user specified
+            }
+            else if (m_Star1->MZAMS() <= 0.7) {                                                                                                             // no - MS - initial mass determines actual type  (don't use utils::Compare() here)
+                if (m_Star1->StellarType() != STELLAR_TYPE::MS_LTE_07) (void)m_Star1->SwitchTo(STELLAR_TYPE::MS_LTE_07, true);                              // MS <= 0.7 Msol - switch if necessary
+            }
+            else {
+                if (m_Star1->StellarType() != STELLAR_TYPE::MS_GT_07) (void)m_Star1->SwitchTo(STELLAR_TYPE::MS_GT_07, true);                                // MS > 0.7 Msol - switch if necessary
+            }
         }
 
         // star 2
-        if (utils::Compare(omega, m_Star2->OmegaCHE()) >= 0) {                                                                                          // star 2 CH?
-            if (m_Star2->StellarType() != STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) {
-                (void)m_Star2->SwitchTo(STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS, true);                                                                    // yes, switch if not already Chemically Homogeneous
-                if (OPTIONS->OptionDefaulted("rotational-frequency-2")) m_Star2->SetOmega(omega); }                                                     // set spin to orbital frequency unless user specified
-        }
-        else if (m_Star2->MZAMS() <= 0.7) {                                                                                                             // no - MS - initial mass determines actual type  (don't use utils::Compare() here)
-            if (m_Star2->StellarType() != STELLAR_TYPE::MS_LTE_07) (void)m_Star2->SwitchTo(STELLAR_TYPE::MS_LTE_07, true);                              // MS <= 0.0 Msol - switch if necessary
-        }
-        else {
-            if (m_Star2->StellarType() != STELLAR_TYPE::MS_GT_07) (void)m_Star2->SwitchTo(STELLAR_TYPE::MS_GT_07, true);                                // MS > 0.7 Msol - switch if necessary
+        if (m_Star2->IsOneOf({ STELLAR_TYPE::MS_LTE_07, STELLAR_TYPE::MS_GT_07, STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS })) {
+            if (utils::Compare(omega, m_Star2->OmegaCHE()) >= 0) {                                                                                          // star 2 CH?
+                if (m_Star2->StellarType() != STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) {
+                    (void)m_Star2->SwitchTo(STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS, true);                                                                    // yes, switch if not already Chemically Homogeneous
+                    if (OPTIONS->OptionDefaulted("rotational-frequency-2")) m_Star2->SetOmega(omega); }                                                     // set spin to orbital frequency unless user specified
+            }
+            else if (m_Star2->MZAMS() <= 0.7) {                                                                                                             // no - MS - initial mass determines actual type  (don't use utils::Compare() here)
+                if (m_Star2->StellarType() != STELLAR_TYPE::MS_LTE_07) (void)m_Star2->SwitchTo(STELLAR_TYPE::MS_LTE_07, true);                              // MS <= 0.0 Msol - switch if necessary
+            }
+            else {
+                if (m_Star2->StellarType() != STELLAR_TYPE::MS_GT_07) (void)m_Star2->SwitchTo(STELLAR_TYPE::MS_GT_07, true);                                // MS > 0.7 Msol - switch if necessary
+            }
         }
 
     }