TeamCOMPAS · ilyamandel · Mar 13, 2025 · Dec 18, 2024 · Dec 18, 2024 · Dec 18, 2024
diff --git a/compas_python_utils/preprocessing/compasConfigDefault.yaml b/compas_python_utils/preprocessing/compasConfigDefault.yaml
@@ -1,5 +1,5 @@
 ##~!!~## COMPAS option values
-##~!!~## File Created Wed Feb 19 11:27:27 2025 by COMPAS v03.13.02
+##~!!~## File Created Wed Mar 12 14:41:45 2025 by COMPAS v03.15.00
 ##~!!~##
 ##~!!~## The default COMPAS YAML file (``compasConfigDefault.yaml``), as distributed, has
 ##~!!~## all COMPAS option entries commented so that the COMPAS default value for the
@@ -82,11 +82,11 @@ numericalChoices:
 #    --hdf5-chunk-size: 100000                                             # Default: 100000
 #    --hdf5-buffer-size: 1                                                 # Default: 1
 #    --log-level: 0                                                        # Default: 0
-#    --mass-change-fraction: 0.000000                                      # Default: 0.000000                                                                                                                                                                         # approximate desired fractional changes in stellar mass per timestep
+#    --mass-change-fraction: 0.001000                                      # Default: 0.001000                                                                                                                                                                         # approximate desired fractional changes in stellar mass per timestep
 #    --maximum-evolution-time: 13700.00                                    # Default: 13700.00                                                                                                                                                                         # maximum physical time a system can be evolved [Myr]
 #    --maximum-number-timestep-iterations: 99999                           # Default: 99999
 #    --number-of-systems: 10                                               # Default: 10                                                                                                                                                                               # number of systems per batch
-#    --radial-change-fraction: 0.000000                                    # Default: 0.000000                                                                                                                                                                         # approximate desired fractional changes in stellar radius per timestep
+#    --radial-change-fraction: 0.100000                                    # Default: 0.100000                                                                                                                                                                         # approximate desired fractional changes in stellar radius per timestep
 #    --timestep-multiplier: 1.000000                                       # Default: 1.000000                                                                                                                                                                         # optional multiplier relative to default time step duration
 
     ### STELLAR PROPERTIES
@@ -209,8 +209,12 @@ numericalChoices:
     ### PULSAR PARAMETERS
 #    --pulsar-birth-magnetic-field-distribution-min: 11.000000             # Default: 11.000000                                                                                                                                                                        # [log10(B/G)]
 #    --pulsar-birth-magnetic-field-distribution-max: 13.000000             # Default: 13.000000                                                                                                                                                                        # [log10(B/G)]
+#    --pulsar-birth-magnetic-field-distribution-mean: 12.650000            # Default: 12.650000                                                                                                                                                                        # [log10(B/G)]
+#    --pulsar-birth-magnetic-field-distribution-sigma: 0.550000            # Default: 0.550000                                                                                                                                                                         # [log10(B/G)]
 #    --pulsar-birth-spin-period-distribution-min: 10.000000                # Default: 10.000000                                                                                                                                                                        # [ms]
 #    --pulsar-birth-spin-period-distribution-max: 100.00                   # Default: 100.00                                                                                                                                                                           # [ms]
+#    --pulsar-birth-spin-period-distribution-mean: 75.000000               # Default: 75.000000                                                                                                                                                                        # [ms]
+#    --pulsar-birth-spin-period-distribution-sigma: 25.000000              # Default: 25.000000                                                                                                                                                                        # [ms]
 #    --pulsar-magnetic-field-decay-timescale: 1000.00                      # Default: 1000.00                                                                                                                                                                          # [Myr]
 #    --pulsar-magnetic-field-decay-massscale: 0.025000                     # Default: 0.025000                                                                                                                                                                         # [Msol]
 #    --pulsar-minimum-magnetic-field: 8.000000                             # Default: 8.000000                                                                                                                                                                         # [log10(B/G)]
@@ -271,13 +275,14 @@ 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-distribution: 'ISOTROPIC'                            # Default: 'ISOTROPIC'                   # Options: ['POLES','WEDGE','POWERLAW','PERPENDICULAR','INPLANE','ISOTROPIC']
+#    --neutron-star-accretion-in-ce: 'ZERO'                                # Default: 'ZERO'                        # Options: ['DISK','SURFACE','ZERO']
 #    --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','ZERO']
-#    --pulsar-birth-spin-period-distribution: 'ZERO'                       # Default: 'ZERO'                        # Options: ['NORMAL','UNIFORM','ZERO']
+#    --pulsar-birth-magnetic-field-distribution: 'LOGNORMAL'               # Default: 'LOGNORMAL'                   # Options: ['LOGNORMAL','UNIFORM','FLATINLOG']
+#    --pulsar-birth-spin-period-distribution: 'NORMAL'                     # Default: 'NORMAL'                      # Options: ['NORMAL','UNIFORM']
 #    --remnant-mass-prescription: 'MULLERMANDEL'                           # Default: 'MULLERMANDEL'                # Options: ['MALTSEV2024','SCHNEIDER2020ALT','SCHNEIDER2020','MULLERMANDEL','MULLER2016','FRYER2022','FRYER2012','BELCZYNSKI2002','HURLEY2000']
 
-    ### LOGFILES AND OUTPUTS
+    ### LOGFILES AND OUTPUTS 
 #    --logfile-type: 'HDF5'                                                # Default: 'HDF5'                        # Options: ['TXT','TSV','CSV','HDF5']
 #    --logfile-name-prefix: ''                                             # Default: ''
 #    --logfile-definitions: ''                                             # Default: ''
@@ -292,7 +297,7 @@ stringChoices:
 #    --output-path: '.'                                                    # Default: '.'
 
 
-listChoices:
+listChoices: 
 
 #    --log-classes: { }                                                    # Default: { }
 #    --debug-classes: { }                                                  # Default: { }
diff --git a/.../COMPAS output/standard-logfiles-default-record-specifications-BSE-detailed.rst b/.../COMPAS output/standard-logfiles-default-record-specifications-BSE-detailed.rst
@@ -71,8 +71,8 @@ Default record definition for the BSE Detailed Output log file::
         BINARY_PROPERTY::MASS_TRANSFER_TRACKER_HISTORY,
         STAR_1_PROPERTY::PULSAR_MAGNETIC_FIELD,
         STAR_2_PROPERTY::PULSAR_MAGNETIC_FIELD,
-        STAR_1_PROPERTY::PULSAR_SPIN_FREQUENCY,
-        STAR_2_PROPERTY::PULSAR_SPIN_FREQUENCY,
+        STAR_1_PROPERTY::PULSAR_SPIN_PERIOD,
+        STAR_2_PROPERTY::PULSAR_SPIN_PERIOD,
         STAR_1_PROPERTY::PULSAR_SPIN_DOWN_RATE,
         STAR_2_PROPERTY::PULSAR_SPIN_DOWN_RATE,
         STAR_1_PROPERTY::PULSAR_BIRTH_PERIOD,

diff --git a/...e/COMPAS output/standard-logfiles-default-record-specifications-BSE-pulsars.rst b/...e/COMPAS output/standard-logfiles-default-record-specifications-BSE-pulsars.rst
@@ -13,8 +13,8 @@ Default record definition for the BSE Pulsar Evolution log file::
         BINARY_PROPERTY::MASS_TRANSFER_TRACKER_HISTORY,
         STAR_1_PROPERTY::PULSAR_MAGNETIC_FIELD,
         STAR_2_PROPERTY::PULSAR_MAGNETIC_FIELD,
-        STAR_1_PROPERTY::PULSAR_SPIN_FREQUENCY,
-        STAR_2_PROPERTY::PULSAR_SPIN_FREQUENCY,
+        STAR_1_PROPERTY::PULSAR_SPIN_PERIOD,
+        STAR_2_PROPERTY::PULSAR_SPIN_PERIOD,
         STAR_1_PROPERTY::PULSAR_SPIN_DOWN_RATE,
         STAR_2_PROPERTY::PULSAR_SPIN_DOWN_RATE,
         STAR_1_PROPERTY::PULSAR_BIRTH_PERIOD,

diff --git a/...e/COMPAS output/standard-logfiles-default-record-specifications-SSE-pulsars.rst b/...e/COMPAS output/standard-logfiles-default-record-specifications-SSE-pulsars.rst
@@ -8,7 +8,7 @@ Default record definition for the SSE Pulsar Evolution log file::
         STAR_PROPERTY::MASS,
         STAR_PROPERTY::STELLAR_TYPE,
         STAR_PROPERTY::PULSAR_MAGNETIC_FIELD,
-        STAR_PROPERTY::PULSAR_SPIN_FREQUENCY,
+        STAR_PROPERTY::PULSAR_SPIN_PERIOD,
         STAR_PROPERTY::PULSAR_SPIN_DOWN_RATE,
         STAR_PROPERTY::PULSAR_BIRTH_PERIOD,
         STAR_PROPERTY::PULSAR_BIRTH_SPIN_DOWN_RATE,

diff --git a/online-docs/pages/User guide/Handling errors/error-table.rst b/online-docs/pages/User guide/Handling errors/error-table.rst
@@ -183,6 +183,8 @@ Following is a list of COMPAS error numbers, corresponding symbolic name, and me
    Reached maximum number of tries when looking for omega when circularising and synchronising for tides
 #. TOO_MANY_PULSAR_SPIN_ITERATIONS |br|
    Reached maximum number of iterations calculating the pulsar birth spin period
+#. TOO_MANY_PULSAR_MAG_ITERATIONS |br|
+   Reached maximum number of iterations calculating the pulsar birth magnetic field
 #. TOO_MANY_REMNANT_MASS_ITERATIONS |br|
    Reached maximum number of iterations when calcuating remnant mass (MULLERMANDEL)
 #. TOO_MANY_RETRIES |br|

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
@@ -976,6 +976,14 @@ 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|
+Assumption about neutron star accretion in CE. |br|
+ZERO indicates no accretion onto NS in CE. |br|
+DISK indicates a RLOF like disk accretion onto NS at Alfven radius. |br|
+SURFACE indicates mass is directly accreted onto the surface of the NS. |br|
+Options: { ZERO, DISK, SURFACE } |br|
+Default = ZERO
+
 **--neutron-star-equation-of-state** |br|
 Neutron star equation of state. |br|
 Options: { SSE, ARP3 } |br|
@@ -1089,8 +1097,8 @@ Default = FALSE
 
 **--pulsar-birth-magnetic-field-distribution** |br|
 Pulsar birth magnetic field distribution. |br|
-Options: { ZERO, FLATINLOG, UNIFORM, LOGNORMAL } |br|
-Default = ZERO
+Options: { FLATINLOG, UNIFORM, LOGNORMAL } |br|
+Default = LOGNORMAL
 
 **--pulsar-birth-magnetic-field-distribution-max** |br|
 Maximum (:math:`log_{10}`) pulsar birth magnetic field. |br|
@@ -1100,10 +1108,18 @@ Default = 13.0
 Minimum (:math:`log_{10}`) pulsar birth magnetic field. |br|
 Default = 11.0
 
+**--pulsar-birth-magnetic-field-distribution-mean** |br|
+Mean of lognormal (:math:`log_{10}`) pulsar birth magnetic field. |br|
+Default = 12.65
+
+**--pulsar-birth-magnetic-field-distribution-sigma** |br|
+Sigma of lognormal (:math:`log_{10}`) pulsar birth magnetic field. |br|
+Default = 0.55
+
 **--pulsar-birth-spin-period-distribution** |br|
 Pulsar birth spin period distribution. |br|
-Options: { ZERO, UNIFORM, NORMAL } |br|
-Default = ZERO
+Options: { UNIFORM, NORMAL } |br|
+Default = NORMAL
 
 **--pulsar-birth-spin-period-distribution-max** |br|
 Maximum pulsar birth spin period (ms). |br|
@@ -1113,6 +1129,14 @@ Default = 100.0
 Minimum pulsar birth spin period (ms). |br|
 Default = 10.0
 
+**--pulsar-birth-spin-period-distribution-mean** |br|
+Mean of normal pulsar birth spin period (ms) distribution. |br|
+Default = 75.0
+
+**--pulsar-birth-spin-period-distribution-sigma** |br|
+Sigma of normal pulsar birth spin period (ms) distribution. |br|
+Default = 25.0
+
 **--pulsar-magnetic-field-decay-massscale** |br|
 Mass scale on which magnetic field decays during accretion (:math:`M_\odot`). |br|
 Default = 0.025

diff --git a/online-docs/pages/whats-new.rst b/online-docs/pages/whats-new.rst
@@ -3,6 +3,27 @@ What's new
 
 Following is a brief list of important updates to the COMPAS code.  A complete record of changes can be found in the file ``changelog.h``.
 
+**03.15.00 Mar 5, 2025**
+
+Changes to the treatment of Neutron Star evolution.
+
+New command line options:
+* `--neutron-star-accretion-in-ce` to determine how a NS accretes mass in a common envelope
+* `--pulsar-birth-magnetic-field-distribution-mean` and `--pulsar-birth-magnetic-field-distribution-sigma` to determine the birth distribution of the pulsar magnetic field (only relevant when the `--pulsar-birth-magnetic-field-distribution` option value is `NORMAL` or `LOGNORMAL`) 
+* `--pulsar-birth-spin-period-distribution-mean` and `--pulsar-birth-spin-period-distribution-sigma` to determine the birth distribution of the pulsar period (only relevant when the `--pulsar-birth-spin-period-distribution` option value is `NORMAL` or `LOGNORMAL`)
+
+Changed command line option values and defaults:
+* `--pulsar-birth-spin-period-distribution` option value `ZERO` is now deprecated.
+* `--pulsar-birth-magnetic-field-distribution` option value `ZERO` is now deprecated.
+* `--pulsar-birth-spin-period-distribution` default option value is now `NORMAL` (was `ZERO`)
+* `--pulsar-birth-magnetic-field-distribution` default option value is now `LOGNORMAL` (was `ZERO`)
+
+Changes to the NS-related values in log files:
+* The pulsar magnetic field strength is now recorded in `Gauss` (was `Tesla`) 
+* The pulsar spin down rate now tracks the pulsar spin period derivative (was spin frequency derivative).
+* The SSE/BSE_Pulsar_Evolution file default record now includes the pulsar spin period (s) instead of spin frequency. Spin frequency is still tracked and can be added using the `logfile-definitions` option
+* The period of non-spinning neutron stars is now reported as infinity instead of zero
+
 **03.14.00 Mar 3, 2025**
 
 * Updates to improve convergence without sacrificing computational speed, including updates to default mass and radial change fractions 

diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp
@@ -1878,11 +1878,14 @@ double BaseBinaryStar::CalculateMassTransferOrbit(const double                 p
             double p_MassDonor0, p_MassAccretor0, p_FractionAccreted, p_IsAccretorDegenerate;
             ode(double massDonor0, double massAccretor0, double fractionAccreted, bool isAccretorDegenerate) : p_MassDonor0(massDonor0), p_MassAccretor0(massAccretor0), p_FractionAccreted(fractionAccreted), p_IsAccretorDegenerate(isAccretorDegenerate) {}
 
-            void operator()(state_type const& x, state_type& dxdt, double p_MassChange ) const {
+            // x is the current state of the ODE (x[0] = semi-major axis a)
+            // dxdm is the change of state wrt mass (dxdm[0] = dadm)
+            // p_MassChange is the cumulative change in mass of the star(s)
+            void operator()(state_type const& x, state_type& dxdm, double p_MassChange ) const {
                 double massD = p_MassDonor0 + p_MassChange;
                 double massA = p_MassAccretor0 - p_MassChange * p_FractionAccreted;
                 double jLoss = CalculateGammaAngularMomentumLoss_Static(massD, massA, p_IsAccretorDegenerate);
-                dxdt[0]      = (-2.0 / massD) * (1.0 - (p_FractionAccreted * (massD / massA)) - ((1.0 - p_FractionAccreted) * (jLoss + 0.5) * (massD / (massA + massD)))) * x[0];
+                dxdm[0]      = (-2.0 / massD) * (1.0 - (p_FractionAccreted * (massD / massA)) - ((1.0 - p_FractionAccreted) * (jLoss + 0.5) * (massD / (massA + massD)))) * x[0];
             }
         };
 
@@ -2521,6 +2524,7 @@ void BaseBinaryStar::ResolveMassChanges() {
 
     // update binary separation, but only if semimajor axis not already infinite and binary does not contain a massless remnant
     // JR: note, this will (probably) fail if option --fp-error-mode is not OFF (the calculation that resulted in m_SemiMajorAxis = inf will (probably) result in a trap)
+    // Maybe use std::isfinite(p_SemiMajorAxis) to ensure p_SemiMajorAxis is not NaN or inf
     if (!isinf(m_SemiMajorAxis) && !HasOneOf({STELLAR_TYPE::MASSLESS_REMNANT}))
         m_SemiMajorAxis = m_SemiMajorAxisPrev + m_aMassLossDiff + m_aMassTransferDiff;
 

diff --git a/src/BaseStar.cpp b/src/BaseStar.cpp
@@ -1111,13 +1111,13 @@ double BaseStar::CalculateLogBindingEnergyLoveridge(bool p_IsMassLoss) const {
  * Binding energy from detailed models (Loveridge et al. 2011) is given in [E]=erg, so use cgs
  *
  *
- * double CalculateLambdaLoveridgeEnergyFormalism(const double p_EnvMass, const double p_IsMassLoss)
+ * double CalculateLambdaLoveridgeEnergyFormalism(const double p_EnvMass, const bool p_IsMassLoss)
  *
  * @param   [IN]    p_EnvMass                   Envelope mass (Msol)
  * @param   [IN]    p_IsMassLoss                Boolean indicating whether mass-loss correction should be applied
  * @return                                      Common envelope lambda parameter
  */
-double BaseStar::CalculateLambdaLoveridgeEnergyFormalism(const double p_EnvMass, const double p_IsMassLoss) const {
+double BaseStar::CalculateLambdaLoveridgeEnergyFormalism(const double p_EnvMass, const bool p_IsMassLoss) const {
 
     double bindingEnergy = PPOW(10.0, CalculateLogBindingEnergyLoveridge(p_IsMassLoss));
     return utils::Compare(bindingEnergy, 0.0) > 0 ? (G_CGS * m_Mass * MSOL_TO_G * p_EnvMass * MSOL_TO_G) / (m_Radius * RSOL_TO_AU * AU_TO_CM * bindingEnergy) : 1.0E-20;

diff --git a/src/BaseStar.h b/src/BaseStar.h
@@ -564,7 +564,7 @@ class BaseStar {
 
     virtual double              CalculateLambdaDewi() const                                                             { return 1.0; }                                                             // Default for stellar types with no LamdaDewi definitions - 1.0 is benign
             double              CalculateLambdaKruckow(const double p_Radius, const double p_Alpha) const;
-            double              CalculateLambdaLoveridgeEnergyFormalism(const double p_EnvMass, const double p_IsMassLoss = false) const;
+            double              CalculateLambdaLoveridgeEnergyFormalism(const double p_EnvMass, const bool p_IsMassLoss = false) const;
     virtual double              CalculateLambdaNanjingStarTrack(const double p_Mass, const double p_Metallicity) const  { return 1.0; }                                                             // Default for stellar types with no LamdaNanjing definitions - 1.0 is benign
     virtual double              CalculateLambdaNanjingEnhanced(const int p_MassIndex, const STELLAR_POPULATION p_StellarPop) const { return 1.0; }                                                  // Default for stellar types with no LamdaNanjing definitions - 1.0 is benign