diff --git a/compas_python_utils/preprocessing/compasConfigDefault.yaml b/compas_python_utils/preprocessing/compasConfigDefault.yaml
index 1f7a25e78..c9d691cb7 100644
--- 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/online-docs/pages/User guide/COMPAS output/standard-logfiles-default-record-specifications-BSE-detailed.rst b/online-docs/pages/User guide/COMPAS output/standard-logfiles-default-record-specifications-BSE-detailed.rst
index 66bca4bc5..b7e377b65 100644
--- a/online-docs/pages/User guide/COMPAS output/standard-logfiles-default-record-specifications-BSE-detailed.rst	
+++ b/online-docs/pages/User guide/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/online-docs/pages/User guide/COMPAS output/standard-logfiles-default-record-specifications-BSE-pulsars.rst b/online-docs/pages/User guide/COMPAS output/standard-logfiles-default-record-specifications-BSE-pulsars.rst
index c346acdd8..be92f1978 100644
--- a/online-docs/pages/User guide/COMPAS output/standard-logfiles-default-record-specifications-BSE-pulsars.rst	
+++ b/online-docs/pages/User guide/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/online-docs/pages/User guide/COMPAS output/standard-logfiles-default-record-specifications-SSE-pulsars.rst b/online-docs/pages/User guide/COMPAS output/standard-logfiles-default-record-specifications-SSE-pulsars.rst
index e26fc9054..66a1e16bc 100644
--- a/online-docs/pages/User guide/COMPAS output/standard-logfiles-default-record-specifications-SSE-pulsars.rst	
+++ b/online-docs/pages/User guide/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
index 8f3e31980..135894862 100644
--- 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
index dce5ceeed..12a20da5f 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	
@@ -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
index 5b457cd13..957658433 100644
--- 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
index b5e63d970..181d2e361 100644
--- 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
index f60af3eee..e3118c338 100755
--- 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
index ca4bd941c..4ea0cc34a 100644
--- 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
 
diff --git a/src/BinaryConstituentStar.cpp b/src/BinaryConstituentStar.cpp
index 8a2eca2de..672bd5c7d 100644
--- a/src/BinaryConstituentStar.cpp
+++ b/src/BinaryConstituentStar.cpp
@@ -398,6 +398,7 @@ void BinaryConstituentStar::SetRocheLobeFlags(const bool p_CommonEnvelope, const
 double BinaryConstituentStar::StarToRocheLobeRadiusRatio(const double p_SemiMajorAxis, const double p_Eccentricity) {
     // binary is unbound or semi-major axis is infinite (evolving single star as binary), so not in RLOF
     // JR: note, this will (probably) fail if option --fp-error-mode is not OFF (the calculation that resulted in p_SemiMajorAxis = inf will (probably) result in a trap)
+    // Maybe use !std::isfinite(p_SemiMajorAxis) to ensure p_SemiMajorAxis is not NaN or inf
     if ((utils::Compare(p_SemiMajorAxis, 0.0) <= 0) || (utils::Compare(p_Eccentricity, 1.0) > 0) || isinf(p_SemiMajorAxis)) return 0.0;
     
     double rocheLobeRadius = BaseBinaryStar::CalculateRocheLobeRadius_Static(Mass(), m_Companion->Mass());
diff --git a/src/BinaryConstituentStar.h b/src/BinaryConstituentStar.h
index 592bbd86c..ee1d830a3 100644
--- a/src/BinaryConstituentStar.h
+++ b/src/BinaryConstituentStar.h
@@ -230,7 +230,8 @@ class BinaryConstituentStar: virtual public Star {
                                                const double p_Epsilon)                  { Star::UpdateMagneticFieldAndSpin(p_CommonEnvelope, 
                                                                                                                            ExperiencedRecycledNS(), 
                                                                                                                            p_Stepsize, 
-                                                                                                                           m_MassTransferDiff * MSOL_TO_KG, p_Epsilon); }
+                                                                                                                           m_MassTransferDiff * MSOL_TO_G,
+                                                                                                                           p_Epsilon); }
 
     void            SetMassLossDiff(const double p_MassLossDiff)                        { m_MassLossDiff = p_MassLossDiff; }                        // JR: todo: better way?  Sanity check?
     void            SetObjectId(const OBJECT_ID p_ObjectId)                             { m_ObjectId = p_ObjectId; }
diff --git a/src/ErrorCatalog.h b/src/ErrorCatalog.h
index fcf66e6e4..a087c8021 100644
--- a/src/ErrorCatalog.h
+++ b/src/ErrorCatalog.h
@@ -128,6 +128,7 @@ enum class ERROR: int {
     TOO_MANY_OMEGA_ITERATIONS,                                      // too many iterations in OMEGA root finder
     TOO_MANY_OMEGA_TRIES,                                           // too many tries in OMEGA root finder
     TOO_MANY_PULSAR_SPIN_ITERATIONS,                                // too many iterations calculating the pulsar birth spin period
+    TOO_MANY_PULSAR_MAG_ITERATIONS,                                 // too many iterations calculating the pulsar birth magnetic field
     TOO_MANY_RETRIES,                                               // generic too many retries
     TOO_MANY_RLOF_ITERATIONS,                                       // too many iterations in RLOF root finder
     TOO_MANY_RLOF_TRIES,                                            // too many tries in RLOF root finder
@@ -171,6 +172,7 @@ enum class ERROR: int {
     UNKNOWN_MT_REJUVENATION_PRESCRIPTION,                           // unknown mass transfer rejuvenation prescription
     UNKNOWN_MT_THERMALLY_LIMITED_VARIATION,                         // unknown mass transfer thermally limited variation
     UNKNOWN_NEUTRINO_MASS_LOSS_PRESCRIPTION,                        // unknown neutrino mass loss prescription
+    UNKNOWN_NS_ACCRETION_IN_CE,                                     // unknown NS accretion in common envelope
     UNKNOWN_NS_EOS,                                                 // unknown NS equation-of-state
     UNKNOWN_OB_MASS_LOSS_PRESCRIPTION,                              // unknown OB mass loss prescription
     UNKNOWN_PPI_PRESCRIPTION,                                       // unknown pulsational pair instability prescription
@@ -344,6 +346,7 @@ const COMPASUnorderedMap<ERROR, std::tuple<ERROR_SCOPE, std::string>> ERROR_CATA
     { ERROR::UNKNOWN_MT_THERMALLY_LIMITED_VARIATION,                { ERROR_SCOPE::ALWAYS,              "Unknown mass transfer thermally limited variation" }},
     { 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_ACCRETION_IN_CE,                            { ERROR_SCOPE::ALWAYS,              "Unknown NS accretion in common envelope" }},
     { ERROR::UNKNOWN_NS_EOS,                                        { ERROR_SCOPE::ALWAYS,              "Unknown NS equation-of-state" }},
     { ERROR::UNKNOWN_OB_MASS_LOSS_PRESCRIPTION,                     { ERROR_SCOPE::ALWAYS,              "Unknown OB mass loss prescription" }},
     { ERROR::UNKNOWN_PPI_PRESCRIPTION,                              { ERROR_SCOPE::ALWAYS,              "Unknown pulsational pair instability prescription" }},
diff --git a/src/LogTypedefs.h b/src/LogTypedefs.h
index 3f4d3028a..df4bf73f6 100644
--- a/src/LogTypedefs.h
+++ b/src/LogTypedefs.h
@@ -860,6 +860,7 @@ enum class PROGRAM_OPTION: int {
 
     NOTES,
 
+    NS_ACCRETION_IN_CE,
     NS_EOS,
 
     ORBITAL_PERIOD,
@@ -1081,6 +1082,7 @@ const COMPASUnorderedMap<PROGRAM_OPTION, std::string> PROGRAM_OPTION_LABEL = {
 
     { PROGRAM_OPTION::NOTES,                                            "NOTES" },
 
+    { PROGRAM_OPTION::NS_ACCRETION_IN_CE,                               "NS_ACCRETION_IN_CE" },
     { PROGRAM_OPTION::NS_EOS,                                           "NS_EOS" },
 
     { PROGRAM_OPTION::ORBITAL_PERIOD,                                   "ORBITAL_PERIOD" },
@@ -1299,12 +1301,12 @@ const std::map<ANY_STAR_PROPERTY, PROPERTY_DETAILS> ANY_STAR_PROPERTY_DETAIL = {
     { ANY_STAR_PROPERTY::OMEGA_ZAMS,                                        { TYPENAME::DOUBLE,           "Omega@ZAMS",                      "Hz",               24, 15}},
     { ANY_STAR_PROPERTY::ORBITAL_ENERGY_POST_SUPERNOVA,                     { TYPENAME::DOUBLE,           "Orbital_Energy>SN",               "Msol^2AU^-1",      24, 15}},
     { ANY_STAR_PROPERTY::ORBITAL_ENERGY_PRE_SUPERNOVA,                      { TYPENAME::DOUBLE,           "Orbital_Energy<SN",               "Msol^2AU^-1",      24, 15}},
-    { ANY_STAR_PROPERTY::PULSAR_MAGNETIC_FIELD,                             { TYPENAME::DOUBLE,           "Pulsar_Mag_Field",                "Tesla",            24, 15}},
-    { ANY_STAR_PROPERTY::PULSAR_SPIN_DOWN_RATE,                             { TYPENAME::DOUBLE,           "Pulsar_Spin_Down",                "rad/s^2",          24, 15}},
+    { ANY_STAR_PROPERTY::PULSAR_MAGNETIC_FIELD,                             { TYPENAME::DOUBLE,           "Pulsar_Mag_Field",                "Gauss",            24, 15}},
+    { ANY_STAR_PROPERTY::PULSAR_SPIN_DOWN_RATE,                             { TYPENAME::DOUBLE,           "Pulsar_Spin_Down",                "s/s",              24, 15}},
     { ANY_STAR_PROPERTY::PULSAR_BIRTH_PERIOD,                               { TYPENAME::DOUBLE,           "Pulsar_Birth_Period",             "s",                24, 15}},
     { ANY_STAR_PROPERTY::PULSAR_BIRTH_SPIN_DOWN_RATE,                       { TYPENAME::DOUBLE,           "Pulsar_Birth_Spin_Down",          "s/s",              24, 15}},
     { ANY_STAR_PROPERTY::PULSAR_SPIN_FREQUENCY,                             { TYPENAME::DOUBLE,           "Pulsar_Spin_Freq",                "rad/s",            24, 15}},
-    { ANY_STAR_PROPERTY::PULSAR_SPIN_PERIOD,                                { TYPENAME::DOUBLE,           "Pulsar_Spin_Period",              "ms",               24, 15}},
+    { ANY_STAR_PROPERTY::PULSAR_SPIN_PERIOD,                                { TYPENAME::DOUBLE,           "Pulsar_Spin_Period",              "s",                24, 15}},
     { ANY_STAR_PROPERTY::RADIAL_EXPANSION_TIMESCALE,                        { TYPENAME::DOUBLE,           "Tau_Radial",                      "Myr",              24, 15}},
     { ANY_STAR_PROPERTY::RADIAL_EXPANSION_TIMESCALE_POST_COMMON_ENVELOPE,   { TYPENAME::DOUBLE,           "Tau_Radial>CE",                   "Myr",              24, 15}},
     { ANY_STAR_PROPERTY::RADIAL_EXPANSION_TIMESCALE_PRE_COMMON_ENVELOPE,    { TYPENAME::DOUBLE,           "Tau_Radial<CE",                   "Myr",              24, 15}},
@@ -1623,6 +1625,7 @@ const std::map<PROGRAM_OPTION, PROPERTY_DETAILS> PROGRAM_OPTION_DETAIL = {
     { PROGRAM_OPTION::NEUTRINO_MASS_LOSS_ASSUMPTION_BH,                         { TYPENAME::INT,        "PO_Neutrino_Mass_Loss_Assmptn",             "-",          4, 1 }},
     { PROGRAM_OPTION::NEUTRINO_MASS_LOSS_VALUE_BH,                              { TYPENAME::DOUBLE,     "PO_Neutrino_Mass_Loss_Value",               "-",         24, 15}},
 
+    { PROGRAM_OPTION::NS_ACCRETION_IN_CE,                                       { TYPENAME::INT,        "PO_NS_ACCRETION_IN_CE",                     "-",          4, 1 }},
     { PROGRAM_OPTION::NS_EOS,                                                   { TYPENAME::INT,        "PO_NS_EOS",                                 "-",          4, 1 }},
 
     { PROGRAM_OPTION::ORBITAL_PERIOD,                                           { TYPENAME::DOUBLE,     "PO_Orbital_Period",                         "days",      24, 15}},
@@ -1969,8 +1972,8 @@ const ANY_PROPERTY_VECTOR BSE_DETAILED_OUTPUT_REC = {
     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,
@@ -2018,8 +2021,8 @@ const ANY_PROPERTY_VECTOR BSE_PULSAR_EVOLUTION_REC = {
     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,
@@ -2204,7 +2207,7 @@ const ANY_PROPERTY_VECTOR SSE_PULSAR_EVOLUTION_REC = {
     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/src/NS.cpp b/src/NS.cpp
index 5e0bf5782..b199db7f0 100755
--- a/src/NS.cpp
+++ b/src/NS.cpp
@@ -136,11 +136,12 @@ DBL_DBL_DBL NS::CalculateCoreCollapseSNParams_Static(const double p_Mass) {
 
 /*
  * Calculate the spin period of a Pulsar at birth according to selected distribution (by commandline option)
- *
+ * Users should note that when choosing the NOSPIN option, 
+ * pulsar spin frequency is set to 0 and spin period is infinity. 
  *
  * double CalculateBirthSpinPeriod()
  *
- * @return                                      Birth spin period of Pulsar in ms
+ * @return                                      Birth spin period of Pulsar in s
  */
 double NS::CalculateBirthSpinPeriod() {
 
@@ -148,10 +149,6 @@ double NS::CalculateBirthSpinPeriod() {
 
     switch (OPTIONS->PulsarBirthSpinPeriodDistribution()) {                                                     // which distribution?
 
-        case PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION::ZERO:                                                       // ZERO
-            pSpin = 0.0;
-            break;
-
         case PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION::UNIFORM: {                                                  // UNIFORM distribution between minimum and maximum value as in Oslowski et al 2011 https://arxiv.org/abs/0903.3538 (default Pmin = and Pmax = )
                                                                                                                 // and also Kiel et al 2008 https://arxiv.org/abs/0805.0059 (default Pmin = 10 ms and Pmax 100 ms, section 3.4)
             double maximum = OPTIONS->PulsarBirthSpinPeriodDistributionMax();
@@ -162,12 +159,12 @@ double NS::CalculateBirthSpinPeriod() {
 
         case PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION::NORMAL: {                                                   // NORMAL distribution from Faucher-Giguere and Kaspi 2006 https://arxiv.org/abs/astro-ph/0512585
 
-            double mean  = 300.0;
-            double sigma = 150.0;
+            double mean  = OPTIONS->PulsarBirthSpinPeriodDistributionMean();
+            double sigma = OPTIONS->PulsarBirthSpinPeriodDistributionSigma();
 
             // this should terminate naturally, but just in case we add a guard
             std::size_t iterations = 0;
-            do { pSpin = RAND->RandomGaussian(sigma) + mean;} while (iterations++ < PULSAR_SPIN_ITERATIONS && utils::Compare(pSpin, 0.0) < 0);
+            do { pSpin = RAND->RandomGaussian(sigma) + mean;} while (iterations++ < PULSAR_SPIN_ITERATIONS && utils::Compare(pSpin, 0.0) <= 0);
             if (iterations >= PULSAR_SPIN_ITERATIONS) THROW_ERROR(ERROR::TOO_MANY_PULSAR_SPIN_ITERATIONS);
 
             } break;
@@ -183,7 +180,7 @@ double NS::CalculateBirthSpinPeriod() {
             THROW_ERROR(ERROR::UNKNOWN_PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION);                                  // throw error
     }
 
-    return pSpin;
+    return pSpin * SECONDS_IN_MS;
 }
 
 
@@ -202,15 +199,11 @@ double NS::CalculateBirthMagneticField() {
 
     switch (OPTIONS->PulsarBirthMagneticFieldDistribution()) {                                                  // which distribution?
 
-        case PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION::ZERO:                                                    // ZERO
-            log10B = 0.0;
-            break;
-
         case PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION::FLATINLOG: {                                             // FLAT IN LOG distribution from Oslowski et al 2011 https://arxiv.org/abs/0903.3538 (log10B0min = , log10B0max = )
 
             double maximum = OPTIONS->PulsarBirthMagneticFieldDistributionMax();
             double minimum = OPTIONS->PulsarBirthMagneticFieldDistributionMin();
-
+            
             log10B = minimum + (RAND->Random() * (maximum - minimum));
 
             } break;
@@ -225,10 +218,15 @@ double NS::CalculateBirthMagneticField() {
 
         case PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION::LOGNORMAL: {                                             // LOG NORMAL distribution from Faucher-Giguere and Kaspi 2006 https://arxiv.org/abs/astro-ph/0512585
 
-            double mean  = 12.65;
-            double sigma = 0.55;
+            double mean  = OPTIONS->PulsarBirthMagneticFieldDistributionMean();
+            double sigma = OPTIONS->PulsarBirthMagneticFieldDistributionSigma();
 
             log10B = RAND->RandomGaussian(sigma) + mean;
+
+            // add a guard to make sure magnetic field is always larger than the value set by --pulsar-minimum-magnetic-field
+            std::size_t iterations = 0;
+            do { log10B = RAND->RandomGaussian(sigma) + mean;} while (iterations++ < PULSAR_MAG_ITERATIONS && utils::Compare(log10B, log10(NS::NS_MAG_FIELD_LOWER_LIMIT)) <= 0);
+            if (iterations >= PULSAR_MAG_ITERATIONS) THROW_ERROR(ERROR::TOO_MANY_PULSAR_MAG_ITERATIONS);
             } break;
 
         default:                                                                                                // unknown prescription
@@ -253,65 +251,60 @@ double NS::CalculateBirthMagneticField() {
  * Uses m_Mass and m_Radius to calculate moment of inertia.
  *
  * Raithel et al. 2016, eq 8 in  https://arxiv.org/abs/1603.06594
- * https://tap.arizona.edu/sites/tap.arizona.edu/files/Raithel_2015_manuscript.pdf
  *
  *
- * double CalculateMomentOfInertiaCGS()
+ * double CalculateMomentOfInertiaCGS_Static(const double p_Mass, const double p_Radius)
  *
- * @return                                      Moment of inertia in g cm^2
+ * @param   [IN]    p_Mass                      Mass of the Neutron Star (g)
+ * @param   [IN]    p_Radius                    Radius of the Neutron Star in (cm)
+ * @return                                      Moment of inertia (g cm^2)
  */
-double NS::CalculateMomentOfInertiaCGS() const {
-
-    // pow() is slow - use multiplication
-
-    double r_km   = m_Radius * RSOL_TO_KM;
-    double m_r    = m_Mass / r_km;
-    double m_r_4  = m_r * m_r * m_r * m_r;
-
-    double r_cm   = m_Radius * RSOL_TO_CM;
-    double r_cm_2 = r_cm * r_cm;
-
-    return 0.237 * m_Mass * MSOL_TO_G * r_cm_2 * (1.0 + (4.2 * m_r) + 90.0 * m_r_4);
+double NS::CalculateMomentOfInertiaCGS_Static(const double p_Mass, const double p_Radius) {
+    double m_r = (p_Mass / MSOL_TO_G) / (p_Radius / KM_TO_CM);
+    return 0.237 * p_Mass * p_Radius * p_Radius * (1.0 + (4.2 * m_r) + 90.0 * m_r * m_r * m_r * m_r);
 }
 
 
 /*
  * Calculate the spin down rate for isolated Neutron Stars in cgs
  *
- * See Equation 2 in https://arxiv.org/pdf/1912.02415.pdf
- *
+ * See Equation 5 in https://arxiv.org/abs/2406.11428
+ * Note that magnetic and rotational axes are orthogonal, leading to sin^2(alpha) = 1 in this equation. 
+ * A model with evolving alpha will be implemented in a future version. 
+ * 
  * Calculates spindown with P and Pdot, then converts to OmegaDot for recording in the log file.
  * Evolution of the inclination between pulsar magnetic and rotational axes will be considered in a future version. 
  *
  * double CalculateSpinDownRate(const double p_Omega, const double p_MomentOfInteria, const double p_MagField, const double p_Radius)
  *
- * @param   [IN]    p_Omega                     Pulsar spin frequency. 
- * @param   [IN]    p_MomentOfInteria           Moment of Interia of the Neutron Star in kg m^2
- * @param   [IN]    p_MagField                  Magnetic field in Tesla
- * @param   [IN]    p_Radius                    Radius of the Neutron Star in kilometres
- * @return                                      Spin down rate (spin frequency derivative) of an isolated Neutron Star in s^(-2)
+ * @param   [IN]    p_Period                    Pulsar spin period (s). 
+ * @param   [IN]    p_MomentOfInteria           Moment of Inertia of the Neutron Star (g cm^2)
+ * @param   [IN]    p_MagField                  Magnetic field (Gauss )
+ * @param   [IN]    p_Radius                    Radius of the Neutron Star (kilometres)
+ * @return                                      Spin down rate (spin period derivative) of an isolated Neutron Star (s^-2)
  */
-double NS::CalculateSpinDownRate(const double p_Omega, const double p_MomentOfInteria, const double p_MagField, const double p_Radius) const {
+double NS::CalculateSpinDownRate(const double p_Period, const double p_MomentOfInteria, const double p_MagField, const double p_Radius) const {
 
    // pow() is slow - use multiplication
 
-   double period            = _2_PI / p_Omega;                                                                  // convert frequency to period
    double cgsRadius         = p_Radius * KM_TO_CM;                                                              // radius in cm
    double radius_6          = cgsRadius * cgsRadius * cgsRadius * cgsRadius * cgsRadius * cgsRadius;
-   double cgsMagField       = p_MagField * TESLA_TO_GAUSS;                                                      // B field in G
-   double magField_2        = cgsMagField * cgsMagField;
+   
+   double magField_2        = p_MagField * p_MagField;
    constexpr double _8_PI_2 = 8.0 * PI_2;
    constexpr double _3_C_3  = 3.0E6 * C * C * C;                                                                // 3.0 * (C * 100.0) * (C * 100.0) * (C * 100.0)
    double pDotTop           = _8_PI_2 * radius_6 * magField_2;
-   double pDotBottom        = _3_C_3 * p_MomentOfInteria * period;
+   double pDotBottom        = _3_C_3 * p_MomentOfInteria * p_Period;
    double pDot              = pDotTop / pDotBottom;                                                             // period derivative 
    
-   return(-pDot * p_Omega / period);                                                                            // convert period derivative to frequency derivative, which is what is recorded in the output
+   return pDot;                                                                           
 }
 
 
 /*
- * Calculates and sets pulsar parameters at birth of pulsar
+ * Calculate and set pulsar parameters at birth of pulsar.
+ * Users should note that when choosing the NOSPIN option, 
+ * pulsar spin frequency is set to 0 and spin period is infinity. 
  *
  * Modifies the following class member variables:
  *
@@ -328,29 +321,32 @@ double NS::CalculateSpinDownRate(const double p_Omega, const double p_MomentOfIn
  */
 void NS::CalculateAndSetPulsarParameters() {
 
-    m_PulsarDetails.magneticField     = PPOW(10.0, CalculateBirthMagneticField()) * GAUSS_TO_TESLA;         // magnetic field in Gauss -> convert to Tesla
-    m_PulsarDetails.spinPeriod        = CalculateBirthSpinPeriod();                                         // 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
+    m_PulsarDetails.magneticField     = PPOW(10.0, CalculateBirthMagneticField());                                                  // magnetic field in Gauss 
+    m_PulsarDetails.spinPeriod        = CalculateBirthSpinPeriod();                                                                 // spin period in s
+    m_MomentOfInertia_CGS             = CalculateMomentOfInertiaCGS();                                                              // MoI in CGS g cm^2
 	
-    // Note we convert neutronStarMomentOfInertia from CGS to SI here
-    m_PulsarDetails.spinDownRate      = CalculateSpinDownRate(m_PulsarDetails.spinFrequency, m_MomentOfInertia_CGS, m_PulsarDetails.magneticField, m_Radius * RSOL_TO_KM);  
+    m_PulsarDetails.spinFrequency     = _2_PI / m_PulsarDetails.spinPeriod;                              
+    m_PulsarDetails.birthPeriod       = m_PulsarDetails.spinPeriod;                                         
+
+    m_PulsarDetails.spinDownRate      = CalculateSpinDownRate(m_PulsarDetails.spinPeriod, m_MomentOfInertia_CGS, m_PulsarDetails.magneticField, m_Radius * RSOL_TO_KM);  
     m_PulsarDetails.birthSpinDownRate = m_PulsarDetails.spinDownRate; 
-    m_AngularMomentum_CGS             = m_MomentOfInertia_CGS * m_PulsarDetails.spinFrequency;              // in CGS g cm^2 s^-1
+    m_AngularMomentum_CGS             = m_MomentOfInertia_CGS * m_PulsarDetails.spinFrequency;                              // in CGS g cm^2 s^-1
 }
 
 
 /*
- * Update the magnetic field and spins of neutron stars when it's deemed to be an isolated pulsar. 
- *
- * This function is called in multiple situations in the NS::UpdateMagneticFieldAndSpin() function
+ * Update the magnetic field and spins of isolated pulsar
+ * Users should note that when pulsar is not spinning, 
+ * this function exits without changing anything. 
+ * Note that we assume the rotational and magnetic axis 
+ * are orthogonal and are not evolved in the current model.
+ * A model with evolving alpha will be implemented in a future version. 
  *
  * Modifies the following class member variables:
  *
  *    m_AngularMomentum_CGS
  *    m_PulsarDetails.spinFrequency
+ *    m_PulsarDetails.spinPeriod
  *    m_PulsarDetails.magneticField
  *    m_PulsarDetails.spinDownRate
  *
@@ -360,55 +356,96 @@ void NS::CalculateAndSetPulsarParameters() {
  * @param   [IN]    p_Stepsize                  Timestep size for integration (in seconds)
  */
 void NS::SpinDownIsolatedPulsar(const double p_Stepsize) {
+
+    double radius_IN_CM      = m_Radius * RSOL_TO_KM * KM_TO_CM;
+    double radius_3          = radius_IN_CM * radius_IN_CM * radius_IN_CM;
+    double radius_6          = radius_3 * radius_3;
+    constexpr double _8_PI_2 = 8.0 * PI_2;
+    constexpr double _3_C_3  = 3.0E6 * C * C * C;                                                                                           // 3.0 * (C * 100.0) * (C * 100.0) * (C * 100.0)
     
-    double NSradius_IN_CM         = m_Radius * RSOL_TO_KM * KM_TO_CM;
-    double NSradius_3             = NSradius_IN_CM * NSradius_IN_CM * NSradius_IN_CM;
-    double NSradius_6             = NSradius_3 * NSradius_3;
-    constexpr double _8_PI_2      = 8.0 * PI_2;
-    constexpr double _3_C_3       = 3.0E6 * C * C * C;                                                                      // 3.0 * (C * 100.0) * (C * 100.0) * (C * 100.0)
-    
-    double initialMagField        = m_PulsarDetails.magneticField;                                                          // (in T)
-    double initialMagField_G      = initialMagField * TESLA_TO_GAUSS;
-    double initialSpinPeriod      = _2_PI / m_PulsarDetails.spinFrequency;
-    double magFieldLowerLimit     = PPOW(10.0, OPTIONS->PulsarLog10MinimumMagneticField()) * GAUSS_TO_TESLA;    
-    double magFieldLowerLimit_G   = magFieldLowerLimit * TESLA_TO_GAUSS;                                   
-    double tau                    = OPTIONS->PulsarMagneticFieldDecayTimescale() * MYR_TO_YEAR * SECONDS_IN_YEAR;                                 
-
-    // calculate isolated decay of the magnetic field for a neutron star
+    double initialMagField   = m_PulsarDetails.magneticField;     
+    double initialSpinPeriod = m_PulsarDetails.spinPeriod;
+ 
+    // calculate the decay of magnetic field for an isolated neutron star
     // see Equation 6 in  arXiv:0903.3538v2       
-    m_PulsarDetails.magneticField = magFieldLowerLimit + (initialMagField - magFieldLowerLimit) * exp(-p_Stepsize / tau);   // update pulsar magnetic field in SI. 
-    
+    if (utils::Compare(initialMagField, NS::NS_MAG_FIELD_LOWER_LIMIT) < 0) {
+        // if magnetic field is already lower than the lower limit, 
+        // set it to the value at the beginning of the timestep.
+        m_PulsarDetails.magneticField = initialMagField;
+    }
+    else {
+        m_PulsarDetails.magneticField = NS::NS_MAG_FIELD_LOWER_LIMIT + (initialMagField - NS::NS_MAG_FIELD_LOWER_LIMIT) * std::exp(-p_Stepsize / NS::NS_DECAY_TIME_SCALE); // update pulsar magnetic field in cgs. 
+    }
     // calculate the spin down rate for isolated neutron stars
-    // see Equation 6 in arxiv:1912.02415
+    // see Equation 3 in arxiv:2406.11428
+    // Note that magnetic and rotational axes are orthogonal, leading to sin^2(alpha) = 1 in this equation. 
     // 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 constant2              = (_8_PI_2 * radius_6) / (_3_C_3 * m_MomentOfInertia_CGS);
+    double term1                  = NS::NS_MAG_FIELD_LOWER_LIMIT * NS::NS_MAG_FIELD_LOWER_LIMIT * p_Stepsize;
+    double term2                  = NS::NS_DECAY_TIME_SCALE * NS::NS_MAG_FIELD_LOWER_LIMIT * ( m_PulsarDetails.magneticField  - initialMagField);
+    double term3                  = (NS::NS_DECAY_TIME_SCALE / 2.0) * ((m_PulsarDetails.magneticField * m_PulsarDetails.magneticField) - (initialMagField * initialMagField));
     double Psquared               = 2.0 * constant2 * (term1 - term2 - term3) + (initialSpinPeriod * initialSpinPeriod);
     
-    double P_f                    = std::sqrt(Psquared);
-    m_PulsarDetails.spinFrequency = _2_PI / P_f;                                                                            // pulsar spin frequency
+    m_PulsarDetails.spinPeriod    = std::sqrt(Psquared);
+    m_PulsarDetails.spinFrequency = _2_PI / m_PulsarDetails.spinPeriod;                                                                     // 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.spinPeriod, m_MomentOfInertia_CGS, m_PulsarDetails.magneticField, m_Radius * RSOL_TO_KM) ; 
 
-    m_AngularMomentum_CGS         = m_PulsarDetails.spinFrequency * m_MomentOfInertia_CGS;                                  // angular momentum of star in CGS
+    m_AngularMomentum_CGS         = m_PulsarDetails.spinFrequency * m_MomentOfInertia_CGS;                                                  // angular momentum of star in CGS
+}
+
+
+/*
+ * Calculate the change in angular momentum wrt mass (dJ/dM) of a neutron star when it accretes mass through RLOF
+ * Note: this function uses CGS units and requires parameters in CGS units where applicable
+ * 
+ * See sec. 2.2.1 in arxiv:1912.02415
+ * 
+ *
+ * double DeltaJByAccretion_Static(const double p_Mass, const double p_Radius_6, const double p_MagField, const double p_SpinFrequency, const double p_mDot, const double p_Epsilon)
+ * 
+ * @param   [IN]    p_Mass                      Initial mass of the NS (g)
+ * @param   [IN]    p_Radius_6                  (Radius of the NS (cm))^6 (for performance - so it isn't calculated at every integration step)
+ * @param   [IN]    p_MagField                  NS magnetic field strength at the beginning of accretion (Gauss)
+ * @param   [IN]    p_SpinFrequency             Angular frequency for the NS at the beginning of accretion (rad/s)
+ * @param   [IN]    p_mDot                      Mass transfer rate (g s^-1)
+ * @param   [IN]    p_Epsilon                   Efficiency factor allowing for uncertainties of coupling magnetic field and matter
+ * @return                                      Change in angular momentum wrt mass (dJ/dM) of NS due to accretion
+ */
+double NS::DeltaJByAccretion_Static(const double p_Mass, const double p_Radius_6, const double p_MagField, const double p_SpinFrequency, const double p_mDot, const double p_Epsilon)  {
+
+    // calculate the Alfven radius for an accreting neutron star
+    // see eq 10 in arxiv:1912.02415 
+    double p              = p_Radius_6 * p_Radius_6 / (p_mDot * p_mDot * p_Mass);
+    double q              = PPOW(p, 1.0 / 7.0);
+    double magneticRadius = ALFVEN_CONST * q * PPOW(p_MagField, 4.0 / 7.0) / 2.0;               // Alfven radius / 2.0 (cm)
+    
+    // calculate the difference in the keplerian angular velocity at the magnetic radius and surface angular velocity of the NS
+    // see eq 2 in 1994MNRAS.269..455J / eq 9 in arxiv:1912.02415 
+    double omegaK = std::sqrt(G_CGS * p_Mass / magneticRadius) / magneticRadius;                // rad/s
+    double vDiff  = omegaK - p_SpinFrequency;                                                   // rad/s
+
+    return p_Epsilon * vDiff * magneticRadius * magneticRadius;                                 // eq 12 in arXiv:0805.0059 / eq 8 in arxiv:1912.02415
 }
 
 
 /*
  * Update the magnetic field and spins of neutron stars in the following situations:
+ * Note 1 : this function uses CGS units and requires parameters in CGS units where applicable
+ * Note 2 : Users should note that when pulsar is not spinning (frequency or magnetic field == 0), 
+ *          pulsar spin period is set to infinity and all other parameters to 0.
  * 
- * 1).  JR: What situations?
+ * 1).  Isolated or post mass-transfer spin-down of neutron star
+ * 2).  Neutron star in interacting binary system experiencing mass-transfer induced spin change for 
+ *      2.1).  Roche Lobe overflow, and 
+ *      2.2).  Common Envelope 
  * 
  * Modifies the following class member variables:
  *
  *    m_AngularMomentum_CGS
+ *    m_MomentOfInertia_CGS            
  *    m_PulsarDetails.spinFrequency
+ *    m_PulsarDetails.spinPeriod
  *    m_PulsarDetails.magneticField
  *    m_PulsarDetails.spinDownRate
  *
@@ -416,73 +453,113 @@ void NS::SpinDownIsolatedPulsar(const double p_Stepsize) {
  * void UpdateMagneticFieldAndSpin(const bool   p_CommonEnvelope, 
  *                                 const bool   p_RecycledNS, 
  *                                 const double p_Stepsize, 
- *                                 const double p_MassGainPerTimeStep, 
+ *                                 const double p_MassGain, 
  *                                 const double p_Epsilon)
  *
- * @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_CommonEnvelope            Boolean flag indicating whether there there is a common envelope
+ * @param   [IN]    p_RecycledNS                Boolean flag indicating whether this star is/was a recycled NS
+ * @param   [IN]    p_Stepsize                  Timestep size for integration (seconds)
+ * @param   [IN]    p_MassGain                  Required accretor mass gain for timestep (g)
  * @param   [IN]    p_Epsilon                   Uncertainty due to mass loss
- * @return                                      Tuple containing the Maximum Mass Acceptance Rate and the Accretion Efficiency Parameter
  */
-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);
+void NS::UpdateMagneticFieldAndSpin(const bool p_CommonEnvelope, const bool p_RecycledNS, double p_Stepsize, double p_MassGain, const double p_Epsilon) {
+
+    if ((!p_RecycledNS && !p_CommonEnvelope) || (!p_RecycledNS && utils::Compare(p_MassGain, 0.0) == 0 )) {                                 // 'classical' isolated pulsars
+        SpinDownIsolatedPulsar(p_Stepsize);                                                                                                 // spin down
     }
-    else if (utils::Compare(m_PulsarDetails.spinFrequency, _2_PI * 1000.0) < 0 && 
-             (p_RecycledNS || p_CommonEnvelope) && utils::Compare(p_MassGainPerTimeStep, 0.0) > 0) {
+    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 massG           = m_Mass * MSOL_TO_G;                                                                                        // mass in g
+        double radiusCM        = m_Radius * RSOL_TO_CM;                                                                                     // radius in cm
+        m_MomentOfInertia_CGS  = CalculateMomentOfInertiaCGS(); 
+        double jAcc            = std::sqrt(G_CGS * massG * radiusCM) * p_MassGain;    
+        m_AngularMomentum_CGS += jAcc;                                                                                                      // angular momentum of the accreted material as it falls onto the surface of the NS
+        if (utils::Compare(m_PulsarDetails.magneticField, NS::NS_MAG_FIELD_LOWER_LIMIT) < 0) {
+            // if magnetic field is already lower than the lower limit, 
+            // set it to the lower limit value.
+            m_PulsarDetails.magneticField = NS::NS_MAG_FIELD_LOWER_LIMIT;
+        }
+        else {
+            m_PulsarDetails.magneticField = (m_PulsarDetails.magneticField - NS::NS_MAG_FIELD_LOWER_LIMIT) * std::exp(-p_MassGain / G_TO_KG / NS::NS_DECAY_MASS_SCALE) + NS::NS_MAG_FIELD_LOWER_LIMIT; // eq. 12 in arxiv:1912.02415 
+        }
         
-        // 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 r_m_6                  =  r_m * r_m * r_m * r_m * r_m * r_m;
+        double previousSpinFrequency  = m_PulsarDetails.spinFrequency;
+        m_PulsarDetails.spinFrequency = m_AngularMomentum_CGS / m_MomentOfInertia_CGS;
+        m_PulsarDetails.spinPeriod    = _2_PI / m_PulsarDetails.spinFrequency;
+        double fDot                   = (m_PulsarDetails.spinFrequency - previousSpinFrequency) / p_Stepsize;
         
-        double MoI_SI                 = m_MomentOfInertia_CGS * CGS_SI;
-        double angularMomentum_SI     = m_AngularMomentum_CGS * CGS_SI;
+        m_PulsarDetails.spinDownRate  = -fDot * m_PulsarDetails.spinPeriod * m_PulsarDetails.spinPeriod / _2_PI;
+    } 
+    else if (utils::Compare(p_MassGain, 0.0) > 0                                                                          &&
+             utils::Compare(m_PulsarDetails.spinPeriod, 0.001) > 0                                                        && 
+             (!p_CommonEnvelope || (p_CommonEnvelope && OPTIONS->NeutronStarAccretionInCE() == NS_ACCRETION_IN_CE::DISK))) {                // pulsar recycling through accretion
+
+        // recycling happens for pulsars with spin period larger than 1 ms and in a binary system with mass transfer
+        // the pulsar being recycled is either in a common envelope, or should have started the recycling process in previous time steps
        
-        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 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);
-        
-        // calculate the difference in the keplerian angular velocity and surface angular velocity of the neutron star in m
-        // see Equation 2 in 1994MNRAS.269..455J       
-        double keplerianVelocityAtAlfvenRadius        = std::sqrt(2.0 * G * mass_kg / alfvenRadius); 
-        double keplerianAngularVelocityAtAlfvenRadius = 4.0 * M_PI * keplerianVelocityAtAlfvenRadius / alfvenRadius;
-        double velocityDifference                     = keplerianAngularVelocityAtAlfvenRadius - m_PulsarDetails.spinFrequency;
-
-        // calculate the change in angular momentum due to accretion
-        // see Equation 12 in arXiv:0805.0059/ Equation 8 in arxiv:1912.02415 
-        double Jdot        = p_Epsilon * velocityDifference * alfvenRadius * alfvenRadius * mDot;
-        angularMomentum_SI = angularMomentum_SI + Jdot * p_Stepsize;
-        
-        if (utils::Compare(angularMomentum_SI / MoI_SI, 0.0) > 0) {
-            m_PulsarDetails.magneticField = newPulsarMagneticField;
-            m_PulsarDetails.spinFrequency = angularMomentum_SI / MoI_SI;
-            m_PulsarDetails.spinDownRate  = Jdot / MoI_SI;
-            m_AngularMomentum_CGS         = angularMomentum_SI / CGS_SI;
-        } 
+        // solve for the angular momentum of the NS after accretion
+        // the accretor will gain p_MassGain g over p_Stepsize seconds
+
+        double initialAngularMomentum_CGS = m_AngularMomentum_CGS;                                                                          // initial angular momentum
+        double initialMagField            = m_PulsarDetails.magneticField;                                                                  // initial magnetic field
+        double initialMass                = m_Mass * MSOL_TO_G;                                                                             // initial mass of NS in g
+
+        double radius                     = m_Radius * RSOL_TO_CM;                                                                          // radius of NS in cm
+        double radius_6                   = radius * radius * radius * radius * radius * radius;                                            // for performance - do it once
+        double mDot                       = p_MassGain / p_Stepsize;                                                                        // required mass transfer rate (g s^-1)
+        double massFinal                  = initialMass + p_MassGain;                                                                       // required final mass of NS (after accretion) in g
+
+
+        // calculate initial mass slice size for integration
+        double jAcc      = DeltaJByAccretion_Static(initialMass, radius_6, m_PulsarDetails.magneticField, m_PulsarDetails.spinFrequency, mDot, p_Epsilon);
+        double massSlice = std::fabs(m_AngularMomentum_CGS / 1000.0 / jAcc);                                                                // abs(Jx10^-3 / dJdM)
+
+        // use the boost ODE solver for speed and accuracy
+
+        // define the ODE
+        // initial state
+        state_type x(1);
+        x[0] = initialAngularMomentum_CGS;                                                                                                  // angular momentum
+
+        // ODE
+        struct ode {
+            double p_Mass, p_Radius, p_Radius_6, p_MagField, p_Mdot, p_Epsilon;
+            ode(double mass, double radius, double radius_6, double magField, double mdot, double epsilon) :
+                p_Mass(mass), p_Radius(radius), p_Radius_6(radius_6), p_MagField(magField), p_Mdot(mdot), p_Epsilon(epsilon) { }
+
+            // x is the current state of the ODE (x[0] = angular momentum J)
+            // dxdm is the change of state wrt mass (dxdm[0] = dJdm)
+            // p_MassDelta is the cumulative change in mass of the NS
+            void operator () (const state_type& x, state_type& dxdm, double p_MassDelta ) const {
+                double m = p_Mass + p_MassDelta;
+                double B = (p_MagField - NS::NS_MAG_FIELD_LOWER_LIMIT) * std::exp(-p_MassDelta / NS::NS_DECAY_MASS_SCALE) + NS::NS_MAG_FIELD_LOWER_LIMIT;
+                double f = x[0] / CalculateMomentOfInertiaCGS_Static(m, p_Radius);
+                dxdm[0]  = DeltaJByAccretion_Static(m, p_Radius_6, B, f, p_Mdot, p_Epsilon);                                    
+            }
+        };
+
+        // integrate
+        controlled_stepper_type stepper;
+        (void)integrate_adaptive(stepper, ode{ initialMass, radius, radius_6, initialMagField, mDot, p_Epsilon }, x, 0.0, p_MassGain, massSlice);
+                
+        // final values
+        m_AngularMomentum_CGS = x[0];
+        m_MomentOfInertia_CGS = CalculateMomentOfInertiaCGS_Static(massFinal, radius);
+        if (utils::Compare(initialMagField, NS::NS_MAG_FIELD_LOWER_LIMIT) < 0) {
+            // if magnetic field is already lower than the lower limit, 
+            // set it to the lower limit value.
+            m_PulsarDetails.magneticField = NS::NS_MAG_FIELD_LOWER_LIMIT;
+        }
         else {
-            SpinDownIsolatedPulsar(p_Stepsize);
-        }        
-    }
-    else  {
-        //In all other conditions, treat the pulsar as isolated. 
-        SpinDownIsolatedPulsar(p_Stepsize);
+            m_PulsarDetails.magneticField = (initialMagField - NS::NS_MAG_FIELD_LOWER_LIMIT) * std::exp(-p_MassGain / NS::NS_DECAY_MASS_SCALE) + NS::NS_MAG_FIELD_LOWER_LIMIT;
+        }
+        m_PulsarDetails.spinFrequency = m_AngularMomentum_CGS / m_MomentOfInertia_CGS;
+        m_PulsarDetails.spinPeriod    = _2_PI / m_PulsarDetails.spinFrequency;
+        double fDot                   = (m_AngularMomentum_CGS - initialAngularMomentum_CGS) / m_MomentOfInertia_CGS / p_Stepsize;          // eq. 11 in arxiv:1912.02415 
+        m_PulsarDetails.spinDownRate  = -fDot * m_PulsarDetails.spinPeriod * m_PulsarDetails.spinPeriod / _2_PI;
+    }      
+    else {                                                                                                                                  // otherwise...    
+        SpinDownIsolatedPulsar(p_Stepsize);                                                                                                 // ...treat the pulsar as isolated - spin down
     }
 }
 
@@ -510,9 +587,6 @@ double NS::ResolveCommonEnvelopeAccretion(const double p_FinalMass,
                                           const double p_CompanionMass,
                                           const double p_CompanionRadius,
                                           const double p_CompanionEnvelope) {
-
-    double deltaMass = CalculateMassAccretedForCO(Mass(), p_CompanionMass, p_CompanionRadius, p_CompanionEnvelope);
-
-    return deltaMass;
+    return CalculateMassAccretedForCO(Mass(), p_CompanionMass, p_CompanionRadius, p_CompanionEnvelope);
 }
 
diff --git a/src/NS.h b/src/NS.h
index 2a6f14091..a232a3ea8 100755
--- a/src/NS.h
+++ b/src/NS.h
@@ -17,10 +17,27 @@ class NS: virtual public BaseStar, public Remnants {
 
 public:
 
-    NS() { m_StellarType = STELLAR_TYPE::NEUTRON_STAR; };
+    NS() {
+        m_StellarType = STELLAR_TYPE::NEUTRON_STAR;                                                                                                             // Set stellar type
+
+        // set NS values based on options provided (so we don't need to do this every timestep)
+        // JR: these are good candidates for the "globals/constants" singleton...
+        ALFVEN_CONST             = PPOW(2.0 * PI_2 / G_CGS, 1.0 / 7.0);
+        NS_MAG_FIELD_LOWER_LIMIT = PPOW(10.0, OPTIONS->PulsarLog10MinimumMagneticField());
+        NS_DECAY_MASS_SCALE      = OPTIONS->PulsarMagneticFieldDecayMassscale() * MSOL_TO_G;
+        NS_DECAY_TIME_SCALE      = OPTIONS->PulsarMagneticFieldDecayTimescale() * MYR_TO_YEAR * SECONDS_IN_YEAR; 
+    };
     
     NS(const BaseStar &p_BaseStar, const bool p_Initialise = true) : BaseStar(p_BaseStar), Remnants(p_BaseStar) {
         m_StellarType = STELLAR_TYPE::NEUTRON_STAR;                                                                                                             // Set stellar type
+
+        // set NS values based on options provided (so we don't need to do this every timestep)
+        // JR: these are good candidates for the "globals/constants" singleton...
+        ALFVEN_CONST             = PPOW(2.0 * PI_2 / G_CGS, 1.0 / 7.0);
+        NS_MAG_FIELD_LOWER_LIMIT = PPOW(10.0, OPTIONS->PulsarLog10MinimumMagneticField());
+        NS_DECAY_MASS_SCALE      = OPTIONS->PulsarMagneticFieldDecayMassscale() * MSOL_TO_G;
+        NS_DECAY_TIME_SCALE      = OPTIONS->PulsarMagneticFieldDecayTimescale() * MYR_TO_YEAR * SECONDS_IN_YEAR; 
+
         if (p_Initialise) Initialise();                                                                                                                         // Initialise if required
     }
 
@@ -37,18 +54,33 @@ class NS: virtual public BaseStar, public Remnants {
     }
 
 
+    // member variables
+
+    // static variables that only need to be calculated once
+    // be aware that these variables are global (because they're static), and are shared amongst all NS instances
+    // JR: these are good candidates for the "globals/constants" singleton...
+
+    static inline double ALFVEN_CONST { 0.0 };                                                                                                                  // Constant for calculating Alfven radius - CGS units (mu0 = 1.0)
+    static inline double NS_MAG_FIELD_LOWER_LIMIT { 0.0 };
+    static inline double NS_DECAY_MASS_SCALE { 0.0 };
+    static inline double NS_DECAY_TIME_SCALE { 0.0 }; 
+
+
     // member functions - alphabetically
     static  DBL_DBL_DBL CalculateCoreCollapseSNParams_Static(const double p_Mass);
 
+    static double       DeltaJByAccretion_Static(const double p_Mass, const double p_Radius_6, const double p_MagField, const double p_SpinFrequency, const double p_mDot, const double p_Epsilon);
+
     MT_CASE             DetermineMassTransferTypeAsDonor() const                { return MT_CASE::NONE; }                                                       // Always NONE
 
 
 protected:
     
     void Initialise() {
+
         CalculateTimescales();                                                                                                                                  // Initialise timescales
         
-        //Set internal properties to zero to avoid meaningless values
+        // set internal properties to zero to avoid meaningless values
         m_Age        = 0.0;
         m_COCoreMass = 0.0;
         m_HeCoreMass = 0.0;
@@ -60,6 +92,8 @@ class NS: virtual public BaseStar, public Remnants {
         CalculateAndSetPulsarParameters();
     }
 
+    // member variables
+
     double m_AngularMomentum_CGS;                                                                                                                               // Current angular momentum in CGS - only required in NS class
     double m_MomentOfInertia_CGS;                                                                                                                               // MoI in CGS - only required in NS class
 
@@ -75,14 +109,15 @@ class NS: virtual public BaseStar, public Remnants {
 
             double          CalculateMassLossRate()                             { return 0.0; }                                                                 // Ensure that NSs don't lose mass in winds
     
-            double          CalculateMomentOfInertiaCGS() const;                                                                                                // MoI in CGS
+    static  double          CalculateMomentOfInertiaCGS_Static(const double p_Mass, const double p_Radius);                                                     // MoI in CGS            
+            double          CalculateMomentOfInertiaCGS() const                 { return CalculateMomentOfInertiaCGS_Static(m_Mass * MSOL_TO_G, m_Radius * RSOL_TO_CM); } // MOI in CGS - use member variables
             double          CalculateMomentOfInertia() const                    { return CalculateMomentOfInertiaCGS() / MSOL_TO_G / RSOL_TO_CM / RSOL_TO_CM; } // MoI (default is solar units)
 
     static  double          CalculateRadiusOnPhaseInKM_Static(const double p_Mass);                                                                             // Radius on phase in km
     static  double          CalculateRadiusOnPhase_Static(const double p_Mass)  { return CalculateRadiusOnPhaseInKM_Static(p_Mass) * KM_TO_RSOL; }              // Radius on phase in Rsol
             double          CalculateRadiusOnPhase() const                      { return CalculateRadiusOnPhase_Static(m_Mass); }                               // Radius on phase in Rsol
     
-            double          CalculateRadiusOnPhase(double p_Mass, double p_Luminosity) const                        { return CalculateRadiusOnPhase(); }        // not a meaningful calculation for NS, ignore arguments
+            double          CalculateRadiusOnPhase(double p_Mass, double p_Luminosity) const { return CalculateRadiusOnPhase(); }                               // not a meaningful calculation for NS, ignore arguments
 
             double          CalculateSpinDownRate(const double p_Omega, const double p_MomentOfInteria, const double p_MagField, const double p_Radius) const;
 
diff --git a/src/Options.cpp b/src/Options.cpp
index 7d8180854..12e783e06 100644
--- a/src/Options.cpp
+++ b/src/Options.cpp
@@ -540,23 +540,32 @@ void Options::OptionValues::Initialise() {
     m_MetallicityDistributionMax                                    = MAXIMUM_METALLICITY;
 
 
+    // Neutron star accretion scenario in common envelope
+    m_NeutronStarAccretionInCE.type                                 = NS_ACCRETION_IN_CE::ZERO;
+    m_NeutronStarAccretionInCE.typeString                           = NS_ACCRETION_IN_CE_LABEL.at(m_NeutronStarAccretionInCE.type);
+
+
     // Neutron star equation of state
     m_NeutronStarEquationOfState.type                               = NS_EOS::SSE;
-    m_NeutronStarEquationOfState.typeString                         = NS_EOSLabel.at(m_NeutronStarEquationOfState.type);
+    m_NeutronStarEquationOfState.typeString                         = NS_EOS_LABEL.at(m_NeutronStarEquationOfState.type);
 
 
     // Pulsar birth magnetic field distribution
-    m_PulsarBirthMagneticFieldDistribution.type                     = PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION::ZERO;
+    m_PulsarBirthMagneticFieldDistribution.type                     = PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION::LOGNORMAL;
     m_PulsarBirthMagneticFieldDistribution.typeString               = PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION_LABEL.at(m_PulsarBirthMagneticFieldDistribution.type);
     m_PulsarBirthMagneticFieldDistributionMin                       = 11.0;
     m_PulsarBirthMagneticFieldDistributionMax                       = 13.0;
+    m_PulsarBirthMagneticFieldDistributionMean                      = 12.65;
+    m_PulsarBirthMagneticFieldDistributionSigma                     = 0.55;
 
 
     // Pulsar birth spin period distribution string
-    m_PulsarBirthSpinPeriodDistribution.type                        = PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION::ZERO;
+    m_PulsarBirthSpinPeriodDistribution.type                        = PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION::NORMAL;
     m_PulsarBirthSpinPeriodDistribution.typeString                  = PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION_LABEL.at(m_PulsarBirthSpinPeriodDistribution.type);
     m_PulsarBirthSpinPeriodDistributionMin                          = 10.0;
     m_PulsarBirthSpinPeriodDistributionMax                          = 100.0;
+    m_PulsarBirthSpinPeriodDistributionMean                         = 75.0;
+    m_PulsarBirthSpinPeriodDistributionSigma                        = 25.0;
 
     m_PulsarMagneticFieldDecayTimescale                             = 1000.0;
     m_PulsarMagneticFieldDecayMassscale                             = 0.025;
@@ -1506,6 +1515,15 @@ bool Options::AddOptions(OptionValues *p_Options, po::options_description *p_Opt
             "pulsar-birth-magnetic-field-distribution-min",                
             po::value<double>(&p_Options->m_PulsarBirthMagneticFieldDistributionMin)->default_value(p_Options->m_PulsarBirthMagneticFieldDistributionMin),                                        
             ("Minimum pulsar birth magnetic field, in log10(Gauss) (default = " + std::to_string(p_Options->m_PulsarBirthMagneticFieldDistributionMin) + ")").c_str()
+        )(
+            "pulsar-birth-magnetic-field-distribution-mean",                
+            po::value<double>(&p_Options->m_PulsarBirthMagneticFieldDistributionMean)->default_value(p_Options->m_PulsarBirthMagneticFieldDistributionMean),                                        
+            ("Mean of normal or lognormal distribution for birth magnetic field (log10 B/G) (default = " + std::to_string(p_Options->m_PulsarBirthMagneticFieldDistributionMean) + ")").c_str()
+        )
+        (
+            "pulsar-birth-magnetic-field-distribution-sigma",                
+            po::value<double>(&p_Options->m_PulsarBirthMagneticFieldDistributionSigma)->default_value(p_Options->m_PulsarBirthMagneticFieldDistributionSigma),                                        
+            ("Standard deviation of normal or lognormal distribution for birth magnetic field (log10 B/G) (default = " + std::to_string(p_Options->m_PulsarBirthMagneticFieldDistributionSigma) + ")").c_str()
         )
         (
             "pulsar-birth-spin-period-distribution-max",                   
@@ -1517,6 +1535,16 @@ bool Options::AddOptions(OptionValues *p_Options, po::options_description *p_Opt
             po::value<double>(&p_Options->m_PulsarBirthSpinPeriodDistributionMin)->default_value(p_Options->m_PulsarBirthSpinPeriodDistributionMin),                                              
             ("Minimum pulsar birth spin period, in ms (default = " + std::to_string(p_Options->m_PulsarBirthSpinPeriodDistributionMin) + ")").c_str()
         )
+        (
+            "pulsar-birth-spin-period-distribution-mean",                   
+            po::value<double>(&p_Options->m_PulsarBirthSpinPeriodDistributionMean)->default_value(p_Options->m_PulsarBirthSpinPeriodDistributionMean),                                              
+            ("Mean of normal or lognormal distribution for birth spin period (ms) (default = " + std::to_string(p_Options->m_PulsarBirthSpinPeriodDistributionMax) + ")").c_str()
+        )
+        (
+            "pulsar-birth-spin-period-distribution-sigma",                   
+            po::value<double>(&p_Options->m_PulsarBirthSpinPeriodDistributionSigma)->default_value(p_Options->m_PulsarBirthSpinPeriodDistributionSigma),                                              
+            ("Standard deviation of normal or lognormal distribution for birth spin period (ms) (default = " + std::to_string(p_Options->m_PulsarBirthSpinPeriodDistributionSigma) + ")").c_str()
+        )
         (
             "pulsar-magnetic-field-decay-massscale",                       
             po::value<double>(&p_Options->m_PulsarMagneticFieldDecayMassscale)->default_value(p_Options->m_PulsarMagneticFieldDecayMassscale),                                                    
@@ -1843,6 +1871,11 @@ bool Options::AddOptions(OptionValues *p_Options, po::options_description *p_Opt
             po::value<std::string>(&p_Options->m_NeutrinoMassLossAssumptionBH.typeString)->default_value(p_Options->m_NeutrinoMassLossAssumptionBH.typeString),                                                  
             ("Assumption about neutrino mass loss during BH formation (" + AllowedOptionValuesFormatted("neutrino-mass-loss-BH-formation") + ", default = '" + p_Options->m_NeutrinoMassLossAssumptionBH.typeString + "')").c_str()
         )
+        (
+            "neutron-star-accretion-in-ce",                              
+            po::value<std::string>(&p_Options->m_NeutronStarAccretionInCE.typeString)->default_value(p_Options->m_NeutronStarAccretionInCE.typeString),                                                      
+            ("Neutron star accretion in common envelope to use (" + AllowedOptionValuesFormatted("neutron-star-accretion-in-ce") + ", default = '" + p_Options->m_NeutronStarAccretionInCE.typeString + "')").c_str()
+        )
         (
             "neutron-star-equation-of-state",                              
             po::value<std::string>(&p_Options->m_NeutronStarEquationOfState.typeString)->default_value(p_Options->m_NeutronStarEquationOfState.typeString),                                                      
@@ -2293,8 +2326,13 @@ std::string Options::OptionValues::CheckAndSetOptions() {
             COMPLAIN_IF(!found, "Unknown Neutrino Mass Loss Assumption");
         }
 
+        if (!DEFAULTED("neutron-star-accretion-in-ce")) {                                                                         // neutron star accretion in common envelope
+            std::tie(found, m_NeutronStarAccretionInCE.type) = utils::GetMapKey(m_NeutronStarAccretionInCE.typeString, NS_ACCRETION_IN_CE_LABEL, m_NeutronStarAccretionInCE.type);
+            COMPLAIN_IF(!found, "Unknown Neutron Star Accretion in Common Envelope");
+        }
+
         if (!DEFAULTED("neutron-star-equation-of-state")) {                                                                         // neutron star equation of state
-            std::tie(found, m_NeutronStarEquationOfState.type) = utils::GetMapKey(m_NeutronStarEquationOfState.typeString, NS_EOSLabel, m_NeutronStarEquationOfState.type);
+            std::tie(found, m_NeutronStarEquationOfState.type) = utils::GetMapKey(m_NeutronStarEquationOfState.typeString, NS_EOS_LABEL, m_NeutronStarEquationOfState.type);
             COMPLAIN_IF(!found, "Unknown Neutron Star Equation of State");
         }
 
@@ -2439,7 +2477,13 @@ std::string Options::OptionValues::CheckAndSetOptions() {
         COMPLAIN_IF(m_OverallWindMassLossMultiplier < 0.0, "Overall wind mass loss multiplier (--overall-wind-mass-loss-multiplier) < 0.0");
 
         COMPLAIN_IF(!DEFAULTED("pulsar-magnetic-field-decay-timescale") && m_PulsarMagneticFieldDecayTimescale <= 0.0, "Pulsar magnetic field decay timescale (--pulsar-magnetic-field-decay-timescale) <= 0");
-        COMPLAIN_IF(!DEFAULTED("pulsar-magnetic-field-decay-massscale") && m_PulsarMagneticFieldDecayMassscale <= 0.0, "Pulsar Magnetic field decay massscale (--pulsar-magnetic-field-decay-massscale) <= 0");
+        COMPLAIN_IF(!DEFAULTED("pulsar-magnetic-field-decay-massscale") && m_PulsarMagneticFieldDecayMassscale <= 0.0, "Pulsar magnetic field decay massscale (--pulsar-magnetic-field-decay-massscale) <= 0");
+
+        COMPLAIN_IF(m_PulsarBirthMagneticFieldDistributionMax <= m_PulsarBirthMagneticFieldDistributionMin, "Pulsar birth magnetic field max (--pulsar-birth-magnetic-field-distribution-max) <= min (--pulsar-birth-magnetic-field-distribution-max)");
+        COMPLAIN_IF(m_PulsarBirthMagneticFieldDistributionMin <= m_PulsarLog10MinimumMagneticField, "Pulsar birth magnetic field min (--pulsar-birth-magnetic-field-distribution-min) <= lower limit (--pulsar-minimum-magnetic-field)");
+        
+        COMPLAIN_IF(m_PulsarBirthSpinPeriodDistributionMax <= m_PulsarBirthSpinPeriodDistributionMin, "Pulsar birth spin period max (--pulsar-birth-spin-period-distribution-max) <= min (--pulsar-birth-spin-period-distribution-max)");
+        COMPLAIN_IF(m_PulsarBirthMagneticFieldDistributionMin <= 0.0, "Pulsar birth magnetic field min (--pulsar-birth-spin-period-distribution-min) <= 0");
 
         COMPLAIN_IF(m_RadialChangeFraction <= 0.0, "Radial change fraction per timestep (--radial-change-fraction) <= 0");
         
@@ -2605,7 +2649,8 @@ std::vector<std::string> Options::AllowedOptionValues(const std::string p_Option
         case _("metallicity-distribution")                          : POPULATE_RET(METALLICITY_DISTRIBUTION_LABEL);                 break;
         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_CE_LABEL);                       break;
+        case _("neutron-star-equation-of-state")                    : POPULATE_RET(NS_EOS_LABEL);                                   break;
         case _("OB-mass-loss-prescription")                         : POPULATE_RET(OB_MASS_LOSS_PRESCRIPTION_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;
@@ -4721,6 +4766,7 @@ COMPAS_VARIABLE Options::OptionValue(const T_ANY_PROPERTY p_Property) const {
 
         case PROGRAM_OPTION::NOTES                                          : value = Notes();                                                              break;
 
+        case PROGRAM_OPTION::NS_ACCRETION_IN_CE                             : value = static_cast<int>(NeutronStarAccretionInCE());                         break;
         case PROGRAM_OPTION::NS_EOS                                         : value = static_cast<int>(NeutronStarEquationOfState());                       break;
 
         case PROGRAM_OPTION::ORBITAL_PERIOD                                 : value = OrbitalPeriod();                                                      break;
diff --git a/src/Options.h b/src/Options.h
index e843e5201..8c611ff4b 100755
--- a/src/Options.h
+++ b/src/Options.h
@@ -211,33 +211,14 @@ class Options {
     //       is only shown once per COMPAS run).
 
     std::vector<std::tuple<std::string, std::string, bool>> deprecatedOptionStrings = {
-        { "black-hole-kicks",                            "black-hole-kicks-mode",                           false },
-        { "chemically-homogeneous-evolution",            "chemically-homogeneous-evolution-mode",           false },
-        { "kick-direction",                              "kick-direction-distribution",                     false },
-        { "luminous-blue-variable-prescription",         "LBV-mass-loss-prescription",                      false },
-        { "mass-transfer",                               "use-mass-transfer",                               false },
-        { "mass-transfer-thermal-limit-accretor",        "mass-transfer-thermal-limit-accretor-multiplier", false },
-        { "OB-mass-loss",                                "OB-mass-loss-prescription",                       false },
-        { "retain-core-mass-during-caseA-mass-transfer", "",                                                false },
-        { "RSG-mass-loss",                               "RSG-mass-loss-prescription",                      false },
-        { "VMS-mass-loss",                               "VMS-mass-loss-prescription",                      false },
-        { "WR-mass-loss",                                "WR-mass-loss-prescription",                       false }
+        { "retain-core-mass-during-caseA-mass-transfer", "", false }
     };
 
     std::vector<std::tuple<std::string, std::string, std::string, bool>> deprecatedOptionValues = {
-        { "critical-mass-ratio-prescription",    "GE20", "GE", false },
-        { "critical-mass-ratio-prescription",    "GE20_IC", "GE_IC", false },
-        { "LBV-mass-loss-prescription",          "NONE", "ZERO", false },
-        { "luminous-blue-variable-prescription", "NONE", "ZERO", false },
+        { "critical-mass-ratio-prescription",          "GE20", "GE", false },
+        { "critical-mass-ratio-prescription",          "GE20_IC", "GE_IC", false },
         { "pulsational-pair-instability-prescription", "COMPAS", "WOOSLEY", false},
-        { "OB-mass-loss",                        "NONE", "ZERO", false },
-        { "OB-mass-loss-prescription",           "NONE", "ZERO", false },
-        { "RSG-mass-loss",                       "NONE", "ZERO", false },
-        { "RSG-mass-loss-prescription",          "NONE", "ZERO", false },
-        { "VMS-mass-loss",                       "NONE", "ZERO", false },
-        { "VMS-mass-loss-prescription",          "NONE", "ZERO", false },
-        { "WR-mass-loss",                        "NONE", "ZERO", false },
-        { "WR-mass-loss-prescription",           "NONE", "ZERO", false }
+	{ "pulsar-birth-spin-period-distribution",     "ZERO", "NOSPIN", false }
     };
 
     // the following vector is used to replace deprecated options in the logfile-definitions file
@@ -504,6 +485,8 @@ class Options {
         "maximum-mass-donor-nandez-ivanova",
         "minimum-secondary-mass",
 
+        "neutron-star-accretion-in-ce",
+
         "orbital-period",
         "orbital-period-distribution",
         "orbital-period-max",
@@ -1101,6 +1084,8 @@ class Options {
 
             double                                              m_mCBUR1;                                                       // Minimum core mass at base of the AGB to avoid fully degenerate CO core formation
 
+            // Neutron star accretion in common envelope
+            ENUM_OPT<NS_ACCRETION_IN_CE>                        m_NeutronStarAccretionInCE;                                     // NS accretion in common envelope
 
             // Neutron star equation of state
             ENUM_OPT<NS_EOS>                                    m_NeutronStarEquationOfState;                                   // NS EOS
@@ -1110,11 +1095,15 @@ class Options {
             ENUM_OPT<PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION>  m_PulsarBirthMagneticFieldDistribution;                         // Birth magnetic field distribution for pulsars
             double                                              m_PulsarBirthMagneticFieldDistributionMin;                      // Minimum birth magnetic field (log10 B/G)
             double                                              m_PulsarBirthMagneticFieldDistributionMax;                      // Maximum birth magnetic field (log10 B/G)
+            double                                              m_PulsarBirthMagneticFieldDistributionMean;                     // Mean of normal or lognormal distribution for birth magnetic field (log10 B/G)
+            double                                              m_PulsarBirthMagneticFieldDistributionSigma;                    // Standard deviation of normal or lognormal distribution for birth magnetic field (log10 B/G)
 
             // Pulsar birth spin period distribution string
             ENUM_OPT<PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION>     m_PulsarBirthSpinPeriodDistribution;                            // Birth spin period distribution for pulsars
             double                                              m_PulsarBirthSpinPeriodDistributionMin;                         // Minimum birth spin period (ms)
             double                                              m_PulsarBirthSpinPeriodDistributionMax;                         // Maximum birth spin period (ms)
+            double                                              m_PulsarBirthSpinPeriodDistributionMean;                        // Mean of normal or lognormal distribution for birth spin period (ms)
+            double                                              m_PulsarBirthSpinPeriodDistributionSigma;                       // Standard deviation of normal or lognormal distribution for birth spin period (ms)
 
             double                                              m_PulsarMagneticFieldDecayTimescale;                            // Timescale on which magnetic field decays (Myr)
             double                                              m_PulsarMagneticFieldDecayMassscale;                            // Mass scale on which magnetic field decays during accretion (solar masses)
@@ -1572,6 +1561,7 @@ class Options {
     NEUTRINO_MASS_LOSS_PRESCRIPTION             NeutrinoMassLossAssumptionBH() const                                    { return OPT_VALUE("neutrino-mass-loss-BH-formation", m_NeutrinoMassLossAssumptionBH.type, true); }
     double                                      NeutrinoMassLossValueBH() const                                         { return OPT_VALUE("neutrino-mass-loss-BH-formation-value", m_NeutrinoMassLossValueBH, true); }
 
+    NS_ACCRETION_IN_CE                          NeutronStarAccretionInCE() const                                        { return OPT_VALUE("neutron-star-accretion-in-ce", m_NeutronStarAccretionInCE.type, true); }
     NS_EOS                                      NeutronStarEquationOfState() const                                      { return OPT_VALUE("neutron-star-equation-of-state", m_NeutronStarEquationOfState.type, true); }
 
     std::string                                 Notes(const size_t p_Idx) const                                         { return OPT_VALUE("notes", m_Notes[p_Idx], true); }
@@ -1602,10 +1592,14 @@ class Options {
     PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION    PulsarBirthMagneticFieldDistribution() const                            { return OPT_VALUE("pulsar-birth-magnetic-field-distribution", m_PulsarBirthMagneticFieldDistribution.type, true); }
     double                                      PulsarBirthMagneticFieldDistributionMax() const                         { return OPT_VALUE("pulsar-birth-magnetic-field-distribution-max", m_PulsarBirthMagneticFieldDistributionMax, true); }
     double                                      PulsarBirthMagneticFieldDistributionMin() const                         { return OPT_VALUE("pulsar-birth-magnetic-field-distribution-min", m_PulsarBirthMagneticFieldDistributionMin, true); }
+    double                                      PulsarBirthMagneticFieldDistributionMean() const                        { return OPT_VALUE("pulsar-birth-magnetic-field-distribution-mean", m_PulsarBirthMagneticFieldDistributionMean, true); }
+    double                                      PulsarBirthMagneticFieldDistributionSigma() const                       { return OPT_VALUE("pulsar-birth-magnetic-field-distribution-sigma", m_PulsarBirthMagneticFieldDistributionSigma, true); }
 
     PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION       PulsarBirthSpinPeriodDistribution() const                               { return OPT_VALUE("pulsar-birth-spin-period-distribution", m_PulsarBirthSpinPeriodDistribution.type, true); }
     double                                      PulsarBirthSpinPeriodDistributionMax() const                            { return OPT_VALUE("pulsar-birth-spin-period-distribution-max", m_PulsarBirthSpinPeriodDistributionMax, true); }
     double                                      PulsarBirthSpinPeriodDistributionMin() const                            { return OPT_VALUE("pulsar-birth-spin-period-distribution-min", m_PulsarBirthSpinPeriodDistributionMin, true); }
+    double                                      PulsarBirthSpinPeriodDistributionMean() const                           { return OPT_VALUE("pulsar-birth-spin-period-distribution-mean", m_PulsarBirthSpinPeriodDistributionMean, true); }
+    double                                      PulsarBirthSpinPeriodDistributionSigma() const                          { return OPT_VALUE("pulsar-birth-spin-period-distribution-sigma", m_PulsarBirthSpinPeriodDistributionSigma, true); }
 
     double                                      PulsarLog10MinimumMagneticField() const                                 { return OPT_VALUE("pulsar-minimum-magnetic-field", m_PulsarLog10MinimumMagneticField, true); }
 
diff --git a/src/changelog.h b/src/changelog.h
index fcac9f75b..a91166565 100644
--- a/src/changelog.h
+++ b/src/changelog.h
@@ -1475,7 +1475,26 @@
 //                                      - Changed order of calls to stellar evolution and wind mass loss in SSE to match BSE
 // 03.14.01   IM - Mar 9, 2025     - Defect Repair:
 //                                      - Added a check to prevent a divide-by-zero error from the previous PR (resolves issue #1345)
+// 03.15.00   YS/JR - Mar 03, 2025 - Defect repairs, Enhancement:
+//                                      - Fixed the issue that during mass transfer, the spin-up of a neutron star sometimes created a negative spin period
+//                                      - Updated NS::UpdateMagneticFieldAndSpin() for spin-up/recycling: added Boost integration of angular momentum of neutron star during mass transfer 
+//                                      - Fix for issue #1002
+//                                      - Fix for issue #1257
+//                                      - Updated references to pulsar calculations. 
+//                                      - Added safeguards to make sure the inputs of birth spin period and magnetic field inputs are valid. If not, raise error messages and stop run. 
+//                                      - Consider neutron star not spinning when spin period is infinity, spin frequency is 0 or magnetic field is 0, and all subsequent pulsar parameters are set to 0.
+//                                      - Changes in program options:
+//                                        1). Added program option "--neutron-star-accretion-in-ce" to account for how a neutron star accretes mass during a common envelope event
+//                                        2). Default pulsar birth spin period distribution is set to NORMAL instead of ZERO; ZERO is now deprecated, and non-spinning pulsars are no longer allowed when evolving pulsars. 
+//                                        3). Added program options "--pulsar-birth-spin-period-distribution-mean" (default 75ms) and "--"pulsar-birth-spin-period-distribution-sigma" (default 25ms) to determine the birth distribution of pulsar period when it's normal or lognormal. 
+//                                        4). Default pulsar birth magnetic field distribution is set to LOGNORMAL instead of ZERO; ZERO is now deprecated, and pulsars with zero magnetic field are no longer allowed when evolving pulsars. 
+//                                        5). New command line options "--pulsar-birth-magnetic-field-distribution-mean" (default 12.65)  and "--"pulsar-birth-magnetic-field-distribution-sigma" (default 0.55) to determine the birth distribution of pulsar magnetic field when it's normal or lognormal. 
+//                                      - Changes to SSE/BSE_Pulsar_Evolution file:
+//                                        1). Pulsar magnetic field strength is now recorded in Gauss instead of Tesla 
+//                                        2). Spin of pulsar is now by default recorded with period (s) instead of frequency (Hz). Spin frequency is still tracked and can be added as an output in the logfiles.
+//                                        3). Spin-down of pulsar (m_PulsarDetails.spinDownRate) is now tracking period derivative (p-dot, s/s) instead of frequency derivative (omega-dot, rad/s^2)
+//                                      - Fixed incorrect declarations of BaseStar::CalculateLambdaLoveridgeEnergyFormalism()
 
-const std::string VERSION_STRING = "03.14.01";
+const std::string VERSION_STRING = "03.15.00";
 
 # endif // __changelog_h__
diff --git a/src/constants.h b/src/constants.h
index 6d5411aff..9b30cb857 100755
--- a/src/constants.h
+++ b/src/constants.h
@@ -278,6 +278,7 @@ constexpr double ADAPTIVE_MASS0_SEARCH_FACTOR_FRAC      = 1.0;
 constexpr int    MULLERMANDEL_REMNANT_MASS_MAX_ITERATIONS = 1000;                                                   // Maximum number of iterations to find remnant mass in GiantBranch::CalculateRemnantMassByMullerMandel()
 
 constexpr int    PULSAR_SPIN_ITERATIONS                 = 100;                                                      // Maximum number of iterations to find pulsar birth spin period in NS::CalculatePulsarBirthSpinPeriod()
+constexpr int    PULSAR_MAG_ITERATIONS                  = 100;                                                      // Maximum number of iterations to find pulsar birth magnetic field in NS::CalculateBirthMagneticField()
 
 constexpr int    SEMI_MAJOR_AXIS_SAMPLES                = 100;                                                      // Maximum number of samples when sampling period/semi-major axis in utils::SampleSemiMajorAxis()
 
diff --git a/src/typedefs.h b/src/typedefs.h
index c98b6fec0..987254d44 100755
--- a/src/typedefs.h
+++ b/src/typedefs.h
@@ -708,9 +708,17 @@ const COMPASUnorderedMap<NEUTRINO_MASS_LOSS_PRESCRIPTION, std::string> NEUTRINO_
     { NEUTRINO_MASS_LOSS_PRESCRIPTION::FIXED_MASS,     "FIXED_MASS" }
 };
 
+// neutron star accretion scenario under common envelope
+enum class NS_ACCRETION_IN_CE: int { ZERO, SURFACE, DISK };
+const COMPASUnorderedMap<NS_ACCRETION_IN_CE, std::string> NS_ACCRETION_IN_CE_LABEL = {
+    { NS_ACCRETION_IN_CE::ZERO,  "ZERO" },
+    { NS_ACCRETION_IN_CE::SURFACE, "SURFACE" },
+    { NS_ACCRETION_IN_CE::DISK, "DISK" },
+};
+
 // neutron star equations of state
 enum class NS_EOS: int { SSE, ARP3 };
-const COMPASUnorderedMap<NS_EOS, std::string> NS_EOSLabel = {
+const COMPASUnorderedMap<NS_EOS, std::string> NS_EOS_LABEL = {
     { NS_EOS::SSE,  "SSE" },
     { NS_EOS::ARP3, "ARP3" }
 };
@@ -775,18 +783,16 @@ const COMPASUnorderedMap<PPI_PRESCRIPTION, std::string> PPI_PRESCRIPTION_LABEL =
 enum class PROGRAM_STATUS: int { SUCCESS, CONTINUE, STOPPED, ERROR_IN_COMMAND_LINE, LOGGING_FAILED, ERROR_UNHANDLED_EXCEPTION };
 
 // pulsar birth magnetic field distributions
-enum class PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION: int { ZERO, FLATINLOG, UNIFORM, LOGNORMAL };
+enum class PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION: int { FLATINLOG, UNIFORM, LOGNORMAL };
 const COMPASUnorderedMap<PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION, std::string> PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION_LABEL = {
-    { PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION::ZERO,      "ZERO" },
     { PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION::FLATINLOG, "FLATINLOG" },
     { PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION::UNIFORM,   "UNIFORM" },
     { PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION::LOGNORMAL, "LOGNORMAL" }
 };
 
 // pulsar birth spin period distributions
-enum class PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION: int { ZERO, UNIFORM, NORMAL };
+enum class PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION: int { UNIFORM, NORMAL };
 const COMPASUnorderedMap<PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION, std::string> PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION_LABEL = {
-    { PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION::ZERO,    "ZERO" },
     { PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION::UNIFORM, "UNIFORM" },
     { PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION::NORMAL,  "NORMAL" }
 };
diff --git a/src/yaml.h b/src/yaml.h
index ffacfd476..d8af41fb2 100644
--- a/src/yaml.h
+++ b/src/yaml.h
@@ -275,8 +275,12 @@ namespace yaml {
             "    ### PULSAR PARAMETERS",
             "    --pulsar-birth-magnetic-field-distribution-min                  # [log10(B/G)]",
             "    --pulsar-birth-magnetic-field-distribution-max                  # [log10(B/G)]",
+            "    --pulsar-birth-magnetic-field-distribution-mean                 # [log10(B/G)]",
+            "    --pulsar-birth-magnetic-field-distribution-sigma                # [log10(B/G)]",
             "    --pulsar-birth-spin-period-distribution-min                     # [ms]",
             "    --pulsar-birth-spin-period-distribution-max                     # [ms]",
+            "    --pulsar-birth-spin-period-distribution-mean                    # [ms]",
+            "    --pulsar-birth-spin-period-distribution-sigma                   # [ms]",
             "    --pulsar-magnetic-field-decay-timescale                         # [Myr]",
             "    --pulsar-magnetic-field-decay-massscale                         # [Msol]",
             "    --pulsar-minimum-magnetic-field                                 # [log10(B/G)]",
@@ -337,6 +341,7 @@ namespace yaml {
             "    --fryer-supernova-engine",
             "    --kick-magnitude-distribution",
             "    --kick-direction-distribution",
+            "    --neutron-star-accretion-in-ce",
             "    --neutron-star-equation-of-state",
             "    --neutrino-mass-loss-BH-formation",
             "    --pulsar-birth-magnetic-field-distribution",