diff --git a/docs/build/doctrees/binary_c_parameters.doctree b/docs/build/doctrees/binary_c_parameters.doctree
index eec593e1d6815ee01beb9bc568682ab241abafc5..96cf0b54c6fdf1e6b95972a386be276e838aa73c 100644
Binary files a/docs/build/doctrees/binary_c_parameters.doctree and b/docs/build/doctrees/binary_c_parameters.doctree differ
diff --git a/docs/build/doctrees/custom_logging_functions.doctree b/docs/build/doctrees/custom_logging_functions.doctree
index 5bffadc340db71d391048a1f7b0ec70f04afac6f..49c9e8a7f23a14b54ebb7a54e9db0a46011dced7 100644
Binary files a/docs/build/doctrees/custom_logging_functions.doctree and b/docs/build/doctrees/custom_logging_functions.doctree differ
diff --git a/docs/build/doctrees/distribution_functions.doctree b/docs/build/doctrees/distribution_functions.doctree
index 525ada99590d8728eb613d2b7ef166a7424241c2..1da5ab487d8e3e6976aa07ca1ccce68e92453225 100644
Binary files a/docs/build/doctrees/distribution_functions.doctree and b/docs/build/doctrees/distribution_functions.doctree differ
diff --git a/docs/build/doctrees/environment.pickle b/docs/build/doctrees/environment.pickle
deleted file mode 100644
index a55afdd521a7e085438117dc19577df0381194c4..0000000000000000000000000000000000000000
Binary files a/docs/build/doctrees/environment.pickle and /dev/null differ
diff --git a/docs/build/doctrees/functions.doctree b/docs/build/doctrees/functions.doctree
index 8bc0c89c20ece60fa5d27b30ce029e2232d1958c..845353ed4f085d68aa65bef2d356bbea88818ed1 100644
Binary files a/docs/build/doctrees/functions.doctree and b/docs/build/doctrees/functions.doctree differ
diff --git a/docs/source/binary_c_parameters.rst b/docs/source/binary_c_parameters.rst
index 4fe498e19b7be769db2dba6961d758d454815b8e..d13f4f30c3e6f22acea55f0168310125a6d6b6a3 100644
--- a/docs/source/binary_c_parameters.rst
+++ b/docs/source/binary_c_parameters.rst
@@ -3,8 +3,8 @@ Binary\_c parameters
 The following chapter contains all the parameters that the current version of binary\_c can handle, along with their descriptions and other properties.
 
 
-This information was obtained by the following binary_c build: 
-	**binary_c git branch**: branch_david	**binary_c git revision**: 5542:20210311:f3401ead4	**Built on**: Mar 22 2021 12:07:51
+This information was obtained by the following binary_c build:
+	**binary_c git branch**: branch_david	**binary_c git revision**: 5891:20210625:0148a2a16	**Built on**: Jun 25 2021 23:33:49
 
 
 Section: stars
@@ -45,23 +45,25 @@ Section: stars
 | **Description**: Equatorial rotational speed of star 1 (km/s). If 0.0, the Hurley et al 2000/2002 prescription is used to set the main-sequence velocity, so for a truly non-rotating star, set to something small (e.g. 0.001). Requires MANUAL_VROT. See also vrot2.
 | **Parameter input type**: Float
 | **Default value**: 0
-| **Macros**: ['VROT_BSE = 0', 'VROT_NON_ROTATING = 1e-10', 'VROT_BREAKUP = -1', 'VROT_SYNC = -2', 'binary_c help for variable : vrot1 <Float>', 'Equatorial rotational speed of star 1 (km/s). If 0.0, the Hurley et al 2000/2002 prescription is used to set the main-sequence velocity, so for a truly non-rotating star, set to something small (e.g. 0.001). Requires MANUAL_VROT. See also vrot2.', 'Default : 0']
+| **Macros**: ['VROT_BSE = 0', 'VROT_BREAKUP = -1', 'VROT_SYNC = -2', 'VROT_NON_ROTATING = -3', 'binary_c help for variable : vrot1 <Float>', 'Equatorial rotational speed of star 1 (km/s). If 0.0, the Hurley et al 2000/2002 prescription is used to set the main-sequence velocity, so for a truly non-rotating star, set to something small (e.g. 0.001). Requires MANUAL_VROT. See also vrot2.', 'Default : 0']
 
 | **Parameter**: vrot2
 | **Description**: Equatorial rotational speed of star 2 (km/s). If 0.0, the Hurley et al 2000/2002 prescription is used to set the main-sequence velocity, so for a truly non-rotating star, set to something small (e.g. 0.001). Requires MANUAL_VROT. See also vrot1.
 | **Parameter input type**: Float
 | **Default value**: 0
-| **Macros**: ['VROT_BSE = 0', 'VROT_NON_ROTATING = 1e-10', 'VROT_BREAKUP = -1', 'VROT_SYNC = -2', 'binary_c help for variable : vrot2 <Float>', 'Equatorial rotational speed of star 2 (km/s). If 0.0, the Hurley et al 2000/2002 prescription is used to set the main-sequence velocity, so for a truly non-rotating star, set to something small (e.g. 0.001). Requires MANUAL_VROT. See also vrot1.', 'Default : 0']
+| **Macros**: ['VROT_BSE = 0', 'VROT_BREAKUP = -1', 'VROT_SYNC = -2', 'VROT_NON_ROTATING = -3', 'binary_c help for variable : vrot2 <Float>', 'Equatorial rotational speed of star 2 (km/s). If 0.0, the Hurley et al 2000/2002 prescription is used to set the main-sequence velocity, so for a truly non-rotating star, set to something small (e.g. 0.001). Requires MANUAL_VROT. See also vrot1.', 'Default : 0']
 
 | **Parameter**: vrot3
 | **Description**: The initial equatorial rotational velocity of star three (in km/s, internally this is star index 2). If 0.0, the Hurley et al 2000/2002 prescription is used to set the main-sequence velocity, so for a truly non-rotating star, set to something small (e.g. 0.001). Requires MANUAL_VROT. See also vrot1,2,4.
 | **Parameter input type**: Float
 | **Default value**: 0
+| **Macros**: ['VROT_BSE = 0', 'VROT_BREAKUP = -1', 'VROT_SYNC = -2', 'VROT_NON_ROTATING = -3']
 
 | **Parameter**: vrot4
 | **Description**: The initial equatorial rotational velocity of star four (in km/s, internally this is star index 3). If 0.0, the Hurley et al 2000/2002 prescription is used to set the main-sequence velocity, so for a truly non-rotating star, set to something small (e.g. 0.001). Requires MANUAL_VROT. See also vrot1,2,3.
 | **Parameter input type**: Float
 | **Default value**: 0
+| **Macros**: ['VROT_BSE = 0', 'VROT_BREAKUP = -1', 'VROT_SYNC = -2', 'VROT_NON_ROTATING = -3']
 
 | **Parameter**: inclination1
 | **Description**: The initial inclination of star one (in degrees).
@@ -127,23 +129,49 @@ Section: stars
 | **Description**: Set the stellar type of star 1 (internal index 0), usually MAIN_SEQUENCE (main sequence). Note that setting the stellar type only works for stars with both age=0 and core_mass=0, i.e. main sequence (hydrogen or helium), white dwarfs, black holes and neutrn stars.
 | **Parameter input type**: Integer
 | **Default value**: 0
-| **Macros**: ['LOW_MASS_MAIN_SEQUENCE = 0', 'LOW_MASS_MS = 0', 'MAIN_SEQUENCE = 1', 'MS = 1', 'HG = 2', 'HERTZSPRUNG_GAP = 2', 'GIANT_BRANCH = 3', 'FIRST_GIANT_BRANCH = 3', 'CHeB = 4', 'CORE_HELIUM_BURNING = 4', 'EAGB = 5', 'EARLY_ASYMPTOTIC_GIANT_BRANCH = 5', 'TPAGB = 6', 'THERMALLY_PULSING_ASYMPTOTIC_GIANT_BRANCH = 6', 'HeMS = 7', 'NAKED_MAIN_SEQUENCE_HELIUM_STAR = 7', 'HeHG = 8', 'NAKED_HELIUM_STAR_HERTZSPRUNG_GAP = 8', 'HeGB = 9', 'NAKED_HELIUM_STAR_GIANT_BRANCH = 9', 'HeWD = 10', 'HELIUM_WHITE_DWARF = 10', 'COWD = 11', 'CARBON_OXYGEN_WHITE_DWARF = 11', 'ONeWD = 12', 'OXYGEN_NEON_WHITE_DWARF = 12', 'NS = 13', 'NEUTRON_STAR = 13', 'BH = 14', 'BLACK_HOLE = 14', 'MASSLESS_REMNANT = 15']
+| **Macros**: ['LOW_MASS_MS = 0', 'MS = 1', 'HG = 2', 'GIANT_BRANCH = 3', 'CHeB = 4', 'EAGB = 5', 'TPAGB = 6', 'HeMS = 7', 'HeHG = 8', 'HeGB = 9', 'HeWD = 10', 'COWD = 11', 'ONeWD = 12', 'NS = 13', 'BH = 14', 'MASSLESS_REMNANT = 15', 'LOW_MASS_MAIN_SEQUENCE = 0', 'MAIN_SEQUENCE = 1', 'HERTZSPRUNG_GAP = 2', 'FIRST_GIANT_BRANCH = 3', 'CORE_HELIUM_BURNING = 4', 'EARLY_ASYMPTOTIC_GIANT_BRANCH = 5', 'THERMALLY_PULSING_ASYMPTOTIC_GIANT_BRANCH = 6', 'NAKED_MAIN_SEQUENCE_HELIUM_STAR = 7', 'NAKED_HELIUM_STAR_HERTZSPRUNG_GAP = 8', 'NAKED_HELIUM_STAR_GIANT_BRANCH = 9', 'HELIUM_WHITE_DWARF = 10', 'CARBON_OXYGEN_WHITE_DWARF = 11', 'OXYGEN_NEON_WHITE_DWARF = 12', 'NEUTRON_STAR = 13', 'BLACK_HOLE = 14', 'STAR_WITH_NO_MASS = 15']
 
 | **Parameter**: stellar_type_2
 | **Description**: Set the stellar type of star 2 (internal index 1), usually MAIN_SEQUENCE (main sequence). Note that setting the stellar type only works for stars with both age=0 and core_mass=0, i.e. main sequence (hydrogen or helium), white dwarfs, black holes and neutrn stars.
 | **Parameter input type**: Integer
 | **Default value**: 0
-| **Macros**: ['LOW_MASS_MAIN_SEQUENCE = 0', 'LOW_MASS_MS = 0', 'MAIN_SEQUENCE = 1', 'MS = 1', 'HG = 2', 'HERTZSPRUNG_GAP = 2', 'GIANT_BRANCH = 3', 'FIRST_GIANT_BRANCH = 3', 'CHeB = 4', 'CORE_HELIUM_BURNING = 4', 'EAGB = 5', 'EARLY_ASYMPTOTIC_GIANT_BRANCH = 5', 'TPAGB = 6', 'THERMALLY_PULSING_ASYMPTOTIC_GIANT_BRANCH = 6', 'HeMS = 7', 'NAKED_MAIN_SEQUENCE_HELIUM_STAR = 7', 'HeHG = 8', 'NAKED_HELIUM_STAR_HERTZSPRUNG_GAP = 8', 'HeGB = 9', 'NAKED_HELIUM_STAR_GIANT_BRANCH = 9', 'HeWD = 10', 'HELIUM_WHITE_DWARF = 10', 'COWD = 11', 'CARBON_OXYGEN_WHITE_DWARF = 11', 'ONeWD = 12', 'OXYGEN_NEON_WHITE_DWARF = 12', 'NS = 13', 'NEUTRON_STAR = 13', 'BH = 14', 'BLACK_HOLE = 14', 'MASSLESS_REMNANT = 15']
+| **Macros**: ['LOW_MASS_MS = 0', 'MS = 1', 'HG = 2', 'GIANT_BRANCH = 3', 'CHeB = 4', 'EAGB = 5', 'TPAGB = 6', 'HeMS = 7', 'HeHG = 8', 'HeGB = 9', 'HeWD = 10', 'COWD = 11', 'ONeWD = 12', 'NS = 13', 'BH = 14', 'MASSLESS_REMNANT = 15', 'LOW_MASS_MAIN_SEQUENCE = 0', 'MAIN_SEQUENCE = 1', 'HERTZSPRUNG_GAP = 2', 'FIRST_GIANT_BRANCH = 3', 'CORE_HELIUM_BURNING = 4', 'EARLY_ASYMPTOTIC_GIANT_BRANCH = 5', 'THERMALLY_PULSING_ASYMPTOTIC_GIANT_BRANCH = 6', 'NAKED_MAIN_SEQUENCE_HELIUM_STAR = 7', 'NAKED_HELIUM_STAR_HERTZSPRUNG_GAP = 8', 'NAKED_HELIUM_STAR_GIANT_BRANCH = 9', 'HELIUM_WHITE_DWARF = 10', 'CARBON_OXYGEN_WHITE_DWARF = 11', 'OXYGEN_NEON_WHITE_DWARF = 12', 'NEUTRON_STAR = 13', 'BLACK_HOLE = 14', 'STAR_WITH_NO_MASS = 15']
 
 | **Parameter**: stellar_type_3
 | **Description**: Set the stellar type of star 3 (internal index 2), usually MAIN_SEQUENCE (main sequence). Note that setting the stellar type only works for stars with both age=0 and core_mass=0, i.e. main sequence (hydrogen or helium), white dwarfs, black holes and neutrn stars.
 | **Parameter input type**: Integer
 | **Default value**: 0
+| **Macros**: ['LOW_MASS_MS = 0', 'MS = 1', 'HG = 2', 'GIANT_BRANCH = 3', 'CHeB = 4', 'EAGB = 5', 'TPAGB = 6', 'HeMS = 7', 'HeHG = 8', 'HeGB = 9', 'HeWD = 10', 'COWD = 11', 'ONeWD = 12', 'NS = 13', 'BH = 14', 'MASSLESS_REMNANT = 15', 'LOW_MASS_MAIN_SEQUENCE = 0', 'MAIN_SEQUENCE = 1', 'HERTZSPRUNG_GAP = 2', 'FIRST_GIANT_BRANCH = 3', 'CORE_HELIUM_BURNING = 4', 'EARLY_ASYMPTOTIC_GIANT_BRANCH = 5', 'THERMALLY_PULSING_ASYMPTOTIC_GIANT_BRANCH = 6', 'NAKED_MAIN_SEQUENCE_HELIUM_STAR = 7', 'NAKED_HELIUM_STAR_HERTZSPRUNG_GAP = 8', 'NAKED_HELIUM_STAR_GIANT_BRANCH = 9', 'HELIUM_WHITE_DWARF = 10', 'CARBON_OXYGEN_WHITE_DWARF = 11', 'OXYGEN_NEON_WHITE_DWARF = 12', 'NEUTRON_STAR = 13', 'BLACK_HOLE = 14', 'STAR_WITH_NO_MASS = 15']
 
 | **Parameter**: stellar_type_4
 | **Description**: Set the stellar type of star 4 (internal index 3), usually MAIN_SEQUENCE (main sequence). Note that setting the stellar type only works for stars with both age=0 and core_mass=0, i.e. main sequence (hydrogen or helium), white dwarfs, black holes and neutrn stars.
 | **Parameter input type**: Integer
 | **Default value**: 0
+| **Macros**: ['LOW_MASS_MS = 0', 'MS = 1', 'HG = 2', 'GIANT_BRANCH = 3', 'CHeB = 4', 'EAGB = 5', 'TPAGB = 6', 'HeMS = 7', 'HeHG = 8', 'HeGB = 9', 'HeWD = 10', 'COWD = 11', 'ONeWD = 12', 'NS = 13', 'BH = 14', 'MASSLESS_REMNANT = 15', 'LOW_MASS_MAIN_SEQUENCE = 0', 'MAIN_SEQUENCE = 1', 'HERTZSPRUNG_GAP = 2', 'FIRST_GIANT_BRANCH = 3', 'CORE_HELIUM_BURNING = 4', 'EARLY_ASYMPTOTIC_GIANT_BRANCH = 5', 'THERMALLY_PULSING_ASYMPTOTIC_GIANT_BRANCH = 6', 'NAKED_MAIN_SEQUENCE_HELIUM_STAR = 7', 'NAKED_HELIUM_STAR_HERTZSPRUNG_GAP = 8', 'NAKED_HELIUM_STAR_GIANT_BRANCH = 9', 'HELIUM_WHITE_DWARF = 10', 'CARBON_OXYGEN_WHITE_DWARF = 11', 'OXYGEN_NEON_WHITE_DWARF = 12', 'NEUTRON_STAR = 13', 'BLACK_HOLE = 14', 'STAR_WITH_NO_MASS = 15']
+
+| **Parameter**: max_stellar_type_1
+| **Description**: The maximum stellar type of star 1 (internal index 0). Evolution is stopped when the star reaches this stellar type. If this is negative, massless remnants are allowed, and the maximum stellar type is the absolute value. 
+| **Parameter input type**: Integer
+| **Default value**: 16
+| **Macros**: ['LOW_MASS_MS = 0', 'MS = 1', 'HG = 2', 'GIANT_BRANCH = 3', 'CHeB = 4', 'EAGB = 5', 'TPAGB = 6', 'HeMS = 7', 'HeHG = 8', 'HeGB = 9', 'HeWD = 10', 'COWD = 11', 'ONeWD = 12', 'NS = 13', 'BH = 14', 'MASSLESS_REMNANT = 15', 'LOW_MASS_MAIN_SEQUENCE = 0', 'MAIN_SEQUENCE = 1', 'HERTZSPRUNG_GAP = 2', 'FIRST_GIANT_BRANCH = 3', 'CORE_HELIUM_BURNING = 4', 'EARLY_ASYMPTOTIC_GIANT_BRANCH = 5', 'THERMALLY_PULSING_ASYMPTOTIC_GIANT_BRANCH = 6', 'NAKED_MAIN_SEQUENCE_HELIUM_STAR = 7', 'NAKED_HELIUM_STAR_HERTZSPRUNG_GAP = 8', 'NAKED_HELIUM_STAR_GIANT_BRANCH = 9', 'HELIUM_WHITE_DWARF = 10', 'CARBON_OXYGEN_WHITE_DWARF = 11', 'OXYGEN_NEON_WHITE_DWARF = 12', 'NEUTRON_STAR = 13', 'BLACK_HOLE = 14', 'STAR_WITH_NO_MASS = 15']
+
+| **Parameter**: max_stellar_type_2
+| **Description**: The maximum stellar type of star 2 (internal index 1). Evolution is stopped when the star reaches this stellar type. If this is negative, massless remnants are allowed, and the maximum stellar type is the absolute value.
+| **Parameter input type**: Integer
+| **Default value**: 16
+| **Macros**: ['LOW_MASS_MS = 0', 'MS = 1', 'HG = 2', 'GIANT_BRANCH = 3', 'CHeB = 4', 'EAGB = 5', 'TPAGB = 6', 'HeMS = 7', 'HeHG = 8', 'HeGB = 9', 'HeWD = 10', 'COWD = 11', 'ONeWD = 12', 'NS = 13', 'BH = 14', 'MASSLESS_REMNANT = 15', 'LOW_MASS_MAIN_SEQUENCE = 0', 'MAIN_SEQUENCE = 1', 'HERTZSPRUNG_GAP = 2', 'FIRST_GIANT_BRANCH = 3', 'CORE_HELIUM_BURNING = 4', 'EARLY_ASYMPTOTIC_GIANT_BRANCH = 5', 'THERMALLY_PULSING_ASYMPTOTIC_GIANT_BRANCH = 6', 'NAKED_MAIN_SEQUENCE_HELIUM_STAR = 7', 'NAKED_HELIUM_STAR_HERTZSPRUNG_GAP = 8', 'NAKED_HELIUM_STAR_GIANT_BRANCH = 9', 'HELIUM_WHITE_DWARF = 10', 'CARBON_OXYGEN_WHITE_DWARF = 11', 'OXYGEN_NEON_WHITE_DWARF = 12', 'NEUTRON_STAR = 13', 'BLACK_HOLE = 14', 'STAR_WITH_NO_MASS = 15']
+
+| **Parameter**: max_stellar_type_3
+| **Description**: The maximum stellar type of star 3 (internal index 2). Evolution is stopped when the star reaches this stellar type. If this is negative, massless remnants are allowed, and the maximum stellar type is the absolute value.
+| **Parameter input type**: Integer
+| **Default value**: 16
+| **Macros**: ['LOW_MASS_MS = 0', 'MS = 1', 'HG = 2', 'GIANT_BRANCH = 3', 'CHeB = 4', 'EAGB = 5', 'TPAGB = 6', 'HeMS = 7', 'HeHG = 8', 'HeGB = 9', 'HeWD = 10', 'COWD = 11', 'ONeWD = 12', 'NS = 13', 'BH = 14', 'MASSLESS_REMNANT = 15', 'LOW_MASS_MAIN_SEQUENCE = 0', 'MAIN_SEQUENCE = 1', 'HERTZSPRUNG_GAP = 2', 'FIRST_GIANT_BRANCH = 3', 'CORE_HELIUM_BURNING = 4', 'EARLY_ASYMPTOTIC_GIANT_BRANCH = 5', 'THERMALLY_PULSING_ASYMPTOTIC_GIANT_BRANCH = 6', 'NAKED_MAIN_SEQUENCE_HELIUM_STAR = 7', 'NAKED_HELIUM_STAR_HERTZSPRUNG_GAP = 8', 'NAKED_HELIUM_STAR_GIANT_BRANCH = 9', 'HELIUM_WHITE_DWARF = 10', 'CARBON_OXYGEN_WHITE_DWARF = 11', 'OXYGEN_NEON_WHITE_DWARF = 12', 'NEUTRON_STAR = 13', 'BLACK_HOLE = 14', 'STAR_WITH_NO_MASS = 15']
+
+| **Parameter**: max_stellar_type_4
+| **Description**: The maximum stellar type of star 4 (internal index 3). Evolution is stopped when the star reaches this stellar type. If this is negative, massless remnants are allowed, and the maximum stellar type is the absolute value.
+| **Parameter input type**: Integer
+| **Default value**: 16
+| **Macros**: ['LOW_MASS_MS = 0', 'MS = 1', 'HG = 2', 'GIANT_BRANCH = 3', 'CHeB = 4', 'EAGB = 5', 'TPAGB = 6', 'HeMS = 7', 'HeHG = 8', 'HeGB = 9', 'HeWD = 10', 'COWD = 11', 'ONeWD = 12', 'NS = 13', 'BH = 14', 'MASSLESS_REMNANT = 15', 'LOW_MASS_MAIN_SEQUENCE = 0', 'MAIN_SEQUENCE = 1', 'HERTZSPRUNG_GAP = 2', 'FIRST_GIANT_BRANCH = 3', 'CORE_HELIUM_BURNING = 4', 'EARLY_ASYMPTOTIC_GIANT_BRANCH = 5', 'THERMALLY_PULSING_ASYMPTOTIC_GIANT_BRANCH = 6', 'NAKED_MAIN_SEQUENCE_HELIUM_STAR = 7', 'NAKED_HELIUM_STAR_HERTZSPRUNG_GAP = 8', 'NAKED_HELIUM_STAR_GIANT_BRANCH = 9', 'HELIUM_WHITE_DWARF = 10', 'CARBON_OXYGEN_WHITE_DWARF = 11', 'OXYGEN_NEON_WHITE_DWARF = 12', 'NEUTRON_STAR = 13', 'BLACK_HOLE = 14', 'STAR_WITH_NO_MASS = 15']
 
 | **Parameter**: probability
 | **Description**: The probability is a weighting applied to the star based on, say, the initial mass function. When running a grid of stars to simulate *all* stars, the summed probability of all the stars should be 1.0.
@@ -182,11 +210,21 @@ Section: stars
 | **Parameter input type**: True|False
 | **Default value**: True
 
+| **Parameter**: disable_debug
+| **Description**: Disables debug output. Only has an effect when DEBUG is 1, which probably requires a rebuild. Default FALSE.
+| **Parameter input type**: True|False
+| **Default value**: False
+
 | **Parameter**: timestep_logging
 | **Description**: Turn on timestep logging (default is False).
 | **Parameter input type**: True|False
 | **Default value**: False
 
+| **Parameter**: rejects_in_log
+| **Description**: Show timestep rejections in the main log (default is False).
+| **Parameter input type**: True|False
+| **Default value**: False
+
 | **Parameter**: vandenHeuvel_logging
 | **Description**: Turn on van den Heuvel logging (default is False).
 | **Parameter input type**: True|False
@@ -292,7 +330,7 @@ Section: stars
 | **Description**: Wind prescription during the TPAGB. 0=Karakas 2002 (a modified Vassiliadis and Wood 1993), 1=Hurley et al 2000/2002 (Vassiliadis and Wood 1993), 2=Reimers, 3=Bloecker, 4=Van Loon,  5=Rob's C-wind (broken?), 6,7=Vassiliadis and Wood 1993 (Karakas,Hurley variants respectively) when C/O>1, 8=Mattsson, 9 = Goldman et al. (2017), 10 = Beasor et al. (2020).
 | **Parameter input type**: Integer
 | **Default value**: 0
-| **Macros**: ['TPAGB_WIND_BLOECKER = 3', 'TPAGB_WIND_VW93_KARAKAS = 0', 'TPAGB_WIND_VW93_ORIG = 1', 'TPAGB_WIND_REIMERS = 2', 'TPAGB_WIND_VAN_LOON = 4', 'TPAGB_WIND_ROB_CWIND = 5', 'TPAGB_WIND_VW93_KARAKAS_CARBON_STARS = 6', 'TPAGB_WIND_VW93_ORIG_CARBON_STARS = 7', 'TPAGB_WIND_MATTSSON = 8', 'TPAGB_WIND_GOLDMAN_ETAL_2017 = 9', 'TPAGB_WIND_BEASOR_ETAL_2020 = 10']
+| **Macros**: ['TPAGB_WIND_VW93_KARAKAS = 0', 'TPAGB_WIND_VW93_ORIG = 1', 'TPAGB_WIND_REIMERS = 2', 'TPAGB_WIND_BLOECKER = 3', 'TPAGB_WIND_VAN_LOON = 4', 'TPAGB_WIND_ROB_CWIND = 5', 'TPAGB_WIND_VW93_KARAKAS_CARBON_STARS = 6', 'TPAGB_WIND_VW93_ORIG_CARBON_STARS = 7', 'TPAGB_WIND_MATTSSON = 8', 'TPAGB_WIND_GOLDMAN_ETAL_2017 = 9', 'TPAGB_WIND_BEASOR_ETAL_2020 = 10']
 
 | **Parameter**: eagbwind
 | **Description**: Wind prescription during the EAGB. 0=BSE (Hurley+2002, based on VW93), 1 = Goldman et al. (2017), 2 = Beasor et al. (2020).
@@ -384,11 +422,18 @@ Section: stars
 | **Default value**: 1
 
 | **Parameter**: BH_prescription
-| **Description**: Black hole mass prescription: relates the mass of a newly formed black hole to its progenitor's (CO) core mass. 0=Hurley et al 2000/2002, 1=Belczynski (early 2000s).
+| **Description**: Black hole mass prescrition: relates the mass of a newly formed black hole to its progenitor's (CO) core mass. 0=Hurley et al 2000/2002, 1=Belczynski (early 2000s).
 | **Parameter input type**: Integer
 | **Default value**: 0
 | **Macros**: ['BH_HURLEY2002 = 0', 'BH_BELCZYNSKI = 1', 'BH_SPERA2015 = 2', 'BH_FRYER12_DELAYED = 3', 'BH_FRYER12_RAPID = 4']
 
+| **Parameter**: PPISN_prescription
+| **Description**: (Pulsational) Pair-Instability Supernova prescription: Relates initial helium core mass of star to whether the star undergoes PPISN or PISN. Requires PPISN flag to be True (see binary_c_parameters.h). 0=no ppisn, 1=Farmer et al 2019.
+| **Parameter input type**: Integer
+| **Default value**: 1
+| **Macros**: ['PPISN_NONE = 0', 'PPISN_FARMER19 = 1']
+| **Extra**: Ignore
+
 | **Parameter**: sn_kick_distribution_II
 | **Description**: Set the distribution of speeds applied to kick type II core collapse supernova systems. 0=fixed, 1=maxwellian (hurley/BSE), 2=custom function (see monte_carlo_kicks.c).
 | **Parameter input type**: Integer
@@ -628,7 +673,7 @@ Section: stars
 | **Description**: Set the speed (if >=0) of, or the algothim (if <0) used to calculate the kick on the companion, if it survives, in a subluminous hybrid He-COWD type Ia explosion. 0 = none, 1 = Liu+2015, 2 = Wheeler+ 1975.
 | **Parameter input type**: Float
 | **Default value**: 0
-| **Macros**: ['SN_IMPULSE_NONE = 0', 'SN_IMPULSE_LIU2015 = 1', 'SN_IMPULSE_WHEELER1975 = 2', 'SN_IMPULSE_WHEELER1975 = 2', 'SN_IMPULSE_WHEELER1975 = 2']
+| **Macros**: ['SN_IMPULSE_NONE = 0', 'SN_IMPULSE_LIU2015 = 1', 'SN_IMPULSE_WHEELER1975 = 2']
 
 | **Parameter**: wd_sigma
 | **Description**: Set the speed at which white dwarfs are kicked when they form, in km/s. Default is zero (i.e. no kick). Requires WD_KICKS.
@@ -639,7 +684,7 @@ Section: stars
 | **Description**: Set the direction of white dwarf kicks. 0 = random, 1 = up, 2 = forward, 3 = backward, 4 = inward, 5 = outward. Requires WD_KICKS.
 | **Parameter input type**: Integer
 | **Default value**: 0
-| **Macros**: ['KICK_RANDOM = 0', 'KICK_FORWARD = 2', 'KICK_BACKWARD = 3', 'KICK_STRAIGHT_UP = 1', 'KICK_INWARD = 4', 'KICK_OUTWARD = 5']
+| **Macros**: ['KICK_RANDOM = 0', 'KICK_STRAIGHT_UP = 1', 'KICK_FORWARD = 2', 'KICK_BACKWARD = 3', 'KICK_INWARD = 4', 'KICK_OUTWARD = 5']
 
 | **Parameter**: wd_kick_when
 | **Description**: Decide when to kick a white dwarf. 0=at birth, 1=at first RLOF, 2=at given pulse number (see wd_kick_pulse_number), 3 at every pulse Requires WD_KICKS.
@@ -690,22 +735,22 @@ Section: stars
 | **Parameter**: delta_mcmin
 | **Description**: A parameter to reduce the minimum core mass for third dredge up to occur on the TPAGB. As used by Izzard and Tout (2004) to increase the amount of dredge up, hence carbon, in Magellanic cloud stars.
 | **Parameter input type**: Float
-| **Default value**: NULL
+| **Default value**: 0
 
 | **Parameter**: lambda_min
 | **Description**: A parameter to increase the efficiency of third dredge up on the TPAGB. The efficiency is lambda * lambda_mult, and setting lambda_min>0 implies that, once Mc>Mcmin (see delta_mcmin) lambda=Max(lambda(fit to Karakas), lambda_min). As used by Izzard and Tout (2004) to increase the amount of dredge up, hence carbon, in Magellanic cloud stars. See also lambda_multiplier.
 | **Parameter input type**: Float
-| **Default value**: NULL
+| **Default value**: 0
 
 | **Parameter**: lambda_multiplier
 | **Description**: A parameter to increase the efficiency of third dredge up on the TPAGB. The efficiency is lambda * lambda_mult, and setting lambda_min>0 implies that, once Mc>Mcmin (see delta_mcmin) lambda=Max(lambda(fit to Karakas), lambda_min). As used by Izzard and Tout (2004) to increase the amount of dredge up, hence carbon, in Magellanic cloud stars.
 | **Parameter input type**: Float
-| **Default value**: NULL
+| **Default value**: 1
 
 | **Parameter**: minimum_envelope_mass_for_third_dredgeup
 | **Description**: The minimum envelope mass for third dredge up on the TPAGB. Early, solar metallicity models by Straniero et al suggested 0.5Msun is typical. However, circumstantial evidence (Izzard et al 2009) as well as newer models by Stancliffe and Karakas suggest that at low metallicity a value nearer zero is more appropriate.
 | **Parameter input type**: Float
-| **Default value**: NULL
+| **Default value**: 0.5
 
 | **Parameter**: mass_of_pmz
 | **Description**: The mass in the partial mixing zone of a TPAGB star, using the Karakas 2012 tables. Ask Carlo Abate for more details, or see the series of papers Abate et al 2012, 2013, 2014. Requires NUCSYN and USE_TABULAR_INTERSHELL_ABUNDANCES_KARAKAS_2012.
@@ -715,18 +760,24 @@ Section: stars
 | **Parameter**: c13_eff
 | **Description**: The "efficiency" of partial mixing in a TPAGB star intershell region, when using the s-process tables of Gallino, Busso, Lugaro et al. as provided by Maria Lugaro for the Izzard et al. 2009 paper. Requires NUCSYN and NUCSYN_S_PROCESS.
 | **Parameter input type**: Float
-| **Default value**: NULL
+| **Default value**: 1
 
 | **Parameter**: mc13_pocket_multiplier
 | **Description**: Multiplies the mass in the partial mixing zone of a TPAGB star, when using the s-process tables of Gallino, Busso, Lugaro et al. as provided by Maria Lugaro for the Izzard et al. 2009 paper. Requires NUCSYN and NUCSYN_S_PROCESS.
 | **Parameter input type**: Float
-| **Default value**: NULL
+| **Default value**: 1
+
+| **Parameter**: tides_convective_damping
+| **Description**: Tidal convective damping algorithm. 0=TIDES_HURLEY2002 Zahn 197x timescales + Hut, as in Hurley et al (2002), 1 = TIDES_ZAHN1989: Zahn 1989 lambdas + Hut.
+| **Parameter input type**: Integer
+| **Default value**: 0
+| **Macros**: ['TIDES_HURLEY2002 = 0', 'TIDES_ZAHN1989 = 1']
 
 | **Parameter**: E2_prescription
-| **Description**: Choose how to calculate the E2 structural parameter (used in tidal timescale calculations). 0=Hurley 1=Izzard (see Siess et al 2013)
+| **Description**: Choose how to calculate the E2 structural parameter (used in tidal timescale calculations). 0=Hurley 1=Izzard (see Siess et al 2013).
 | **Parameter input type**: Integer
 | **Default value**: 0
-| **Macros**: ['E2_HURLEY_2002 = 0', 'E2_IZZARD = 1']
+| **Macros**: ['E2_HURLEY_2002 = 0', 'E2_IZZARD = 1', 'E2_MINT = 2']
 
 | **Parameter**: dtfac
 | **Description**: A parameter to decrease the timestep ONLY during the TPAGB phase.
@@ -736,7 +787,7 @@ Section: stars
 | **Parameter**: hbbtfac
 | **Description**: A parameter to modulate the temperature at the base of the hot-bottom burning zone in TPAGB stars. (Works only if NUCSYN is defined)
 | **Parameter input type**: Float
-| **Default value**: NULL
+| **Default value**: 1
 
 | **Parameter**: wind_multiplier_%d
 | **Description**: Wind multiplier for the stellar type specified by %d. By default these are all 1.0.
@@ -800,6 +851,11 @@ Section: stars
 | **Default value**: 0
 | **Macros**: ['NOVA_RETENTION_ALGORITHM_CONSTANT = 0', 'NOVA_RETENTION_ALGORITHM_CLAEYS2014 = 1', 'NOVA_RETENTION_ALGORITHM_HILLMAN2015 = 2']
 
+| **Parameter**: MINT_metallicity
+| **Description**: This sets the metallicity for MINT. It is ignored if set to -1.0, the default, in which case the normal metallicity parameter is used.
+| **Parameter input type**: Float
+| **Default value**: NULL
+
 | **Parameter**: gaia_Teff_binwidth
 | **Description**: log10(Effective temperature) bin width used to make Gaia-like HRDs
 | **Parameter input type**: Float
@@ -817,33 +873,33 @@ Section: stars
 | **Macros**: ['GAIA_CONVERSION_UBVRI_UNIVARIATE_JORDI2010 = 0', 'GAIA_CONVERSION_UBVRI_BIVARIATE_JORDI2010 = 1', 'GAIA_CONVERSION_ugriz_UNIVARIATE_JORDI2010 = 2', 'GAIA_CONVERSION_ugriz_BIVARIATE_JORDI2010 = 3', 'GAIA_CONVERSION_UBVRI_UNIVARIATE_EVANS2018 = 4', 'GAIA_CONVERSION_ugriz_UNIVARIATE_EVANS2018 = 5', 'GAIA_CONVERSION_UBVRI_RIELLO2020 = 6', 'GAIA_CONVERSION_ugriz_RIELLO2020 = 7']
 
 | **Parameter**: rotationally_enhanced_mass_loss
-| **Description**: Set to 1 to enable rotationally enhanced mass loss rate algorithms: 0= none, 1=formula cf. Langer models (=ROTATIONALLY_ENHNACED_MASSLOSS_LANGER_FORMULA), 2=limit accretion rate before wind loss is applied, 3 = both 1 and 2. See also rotationally_enhanced_exponent
+| **Description**: Set to 1 to enable rotationally enhanced mass loss rate algorithms: 0= none, 1=formula cf. Langer models (=ROTATIONALLY_ENHANCED_MASSLOSS_LANGER_FORMULA), 2=limit accretion rate before wind loss is applied, 3 = both 1 and 2. See also rotationally_enhanced_exponent
 | **Parameter input type**: Integer
 | **Default value**: 0
-| **Macros**: ['ROTATIONALLY_ENHNACED_MASSLOSS_NONE = 0', 'ROTATIONALLY_ENHNACED_MASSLOSS_LANGER_FORMULA = 1', 'ROTATIONALLY_ENHNACED_MASSLOSS_ANGMOM = 2', 'ROTATIONALLY_ENHNACED_MASSLOSS_LANGER_FORMULA_AND_ANGMOM = 3']
+| **Macros**: ['ROTATIONALLY_ENHANCED_MASSLOSS_NONE = 0', 'ROTATIONALLY_ENHANCED_MASSLOSS_LANGER_FORMULA = 1', 'ROTATIONALLY_ENHANCED_MASSLOSS_ANGMOM = 2', 'ROTATIONALLY_ENHANCED_MASSLOSS_LANGER_FORMULA_AND_ANGMOM = 3']
 
 | **Parameter**: AGB_core_algorithm
 | **Description**: Algorithm to use for calculating AGB core masses. 0=Hurley et al. 2002 if no NUCSYN, Karakas 2002 if NUCSYN is defined; 1=Hurley et al. 2002 (overshooting models); 1=Karakas 2002 (non-overshooting models).
 | **Parameter input type**: Integer
-| **Default value**: 1
+| **Default value**: 2
 | **Macros**: ['AGB_CORE_ALGORITHM_DEFAULT = 0', 'AGB_CORE_ALGORITHM_HURLEY = 1', 'AGB_CORE_ALGORITHM_KARAKAS = 2']
 
 | **Parameter**: AGB_radius_algorithm
 | **Description**: Algorithm to use for calculating radii on the TPAGB.
 | **Parameter input type**: Integer
-| **Default value**: 1
+| **Default value**: 2
 | **Macros**: ['AGB_RADIUS_ALGORITHM_DEFAULT = 0', 'AGB_RADIUS_ALGORITHM_HURLEY = 1', 'AGB_RADIUS_ALGORITHM_KARAKAS = 2']
 
 | **Parameter**: AGB_luminosity_algorithm
 | **Description**: Algorithm to use for calculating luminosities on the TPAGB.
 | **Parameter input type**: Integer
-| **Default value**: 1
+| **Default value**: 2
 | **Macros**: ['AGB_LUMINOSITY_ALGORITHM_DEFAULT = 0', 'AGB_LUMINOSITY_ALGORITHM_HURLEY = 1', 'AGB_LUMINOSITY_ALGORITHM_KARAKAS = 2']
 
 | **Parameter**: AGB_3dup_algorithm
 | **Description**: Algorithm to use for calculating third dredge up efficiency on the TPAGB.
 | **Parameter input type**: Integer
-| **Default value**: 1
+| **Default value**: 2
 | **Macros**: ['AGB_THIRD_DREDGE_UP_ALGORITHM_DEFAULT = 0', 'AGB_THIRD_DREDGE_UP_ALGORITHM_HURLEY = 1', 'AGB_THIRD_DREDGE_UP_ALGORITHM_KARAKAS = 2', 'AGB_THIRD_DREDGE_UP_ALGORITHM_STANCLIFFE = 3']
 
 | **Parameter**: overspin_algorithm
@@ -877,19 +933,6 @@ Section: stars
 | **Parameter input type**: Boolean(scanf)
 | **Default value**: NULL
 
-| **Parameter**: PPISN_prescription
-| **Description**: (Pulsational) Pair-Instability Supernova prescription: Relates initial helium core mass of star to whether the star undergoes PPISN or PISN. Requires PPISN. 0=no ppisn, 1=farmer 2019, 2=Marchant 2018, 3=Woosley 2019
-| **Parameter input type**: Integer
-| **Default value**: 1
-| **Macros**: ['PPISN_DISABLED = 0', 'PPISN_FARMER19 = 1']
-| **Extra**: Ignore
-
-| **Parameter**: angmom_to_orbit_factor
-| **Description**: Parameter to control the fraction of the excess angular momentum thats returned to the orbit during a disk mass transfer episode
-| **Parameter input type**: Float
-| **Default value**: NULL
-| **Extra**: Ignore
-
 Section: binary
 ---------------
 
@@ -1054,6 +1097,37 @@ Section: binary
 | **Default value**: /tmp/
 | **Extra**: /tmp/
 
+| **Parameter**: post_ce_adaptive_menv
+| **Description**: If TRUE, and if post_ce_objects_have_envelopes is TRUE, then the envelope mass of a post-CE star is such that it sits just inside its Roche lobe. If FALSE then a fixed (thin) envelope mass is applied that depends on the stellar type (see macros POST_CE_ENVELOPE_DM_GB, POST_CE_ENVELOPE_DM_EAGB and POST_CE_ENVELOPE_DM_TPAGB).
+| **Parameter input type**: True|False
+| **Default value**: False
+
+| **Parameter**: post_ce_objects_have_envelopes
+| **Description**: If TRUE then post-common-envelope objects have thin envelopes. You need this if you are to have post-CE post-AGB stars. Note that this *may* be unstable, i.e. you may end up having many CEEs. The mass in the envelope is controlled by post_ce_adaptive_menv. TRUE by default.
+| **Parameter input type**: True|False
+| **Default value**: True
+
+| **Parameter**: PN_comenv_transition_time
+| **Description**: post-common envelope transition time in years (1e2).  This is the time taken to move from CEE ejection to Teff > 30e4 K. Hall et al. (2013) suggest ~100 years.
+| **Parameter input type**: Float
+| **Default value**: 100
+
+| **Parameter**: minimum_time_between_PNe
+| **Description**: The minimum time (Myr) between planetary nebula detections. This prevents multiple, fast common envelopes triggering two PNe (0.1).
+| **Parameter input type**: Float
+| **Default value**: 0.1
+
+| **Parameter**: PN_Hall_fading_time_algorithm
+| **Description**: In stars with low mass (<0.45Msun) cores, you can choose to set the PN fading time to either the minimum (PN_HALL_FADING_TIME_ALGORITHM_MINIMUM) or maximum (PN_HALL_FADING_TIME_ALGORITHM_MAXIMUM) as shown in Fig. 6 of Hall et al. (2013).
+| **Parameter input type**: Integer
+| **Default value**: 0
+| **Macros**: ['PN_HALL_FADING_TIME_ALGORITHM_MINIMUM = 0', 'PN_HALL_FADING_TIME_ALGORITHM_MAXIMUM = 1']
+
+| **Parameter**: PPN_envelope_mass
+| **Description**: Desired pre-planetary nebula (post-AGB) envelope mass.
+| **Parameter input type**: Float
+| **Default value**: 0.01
+
 | **Parameter**: cbdisc_eccentricity_pumping_method
 | **Description**: Select from various eccentricity-pumping methods when there is a circumbinary disc. Requires DISCS. 0 = off.
 | **Parameter input type**: Integer
@@ -1124,7 +1198,7 @@ Section: binary
 | **Parameter**: multiplicity
 | **Description**: Multiplicity: 1=single star, 2=binary, 3=triple, 4=quadruple.
 | **Parameter input type**: Integer
-| **Default value**: 2
+| **Default value**: 0
 
 | **Parameter**: accretion_limit_eddington_steady_multiplier
 | **Description**: Steady accretion is limited by the Eddington instability, with limiting rate given by the accretion_limit_eddington_steady_multiplier * the normal (spherically symmetric) Eddington rate. This is known in the trade as the Eddington factor, and anything greater than 1.0 potentially gives you super-Eddington accretion.
@@ -1201,19 +1275,19 @@ Section: binary
 | **Parameter**: type_Ia_MCh_supernova_algorithm
 | **Description**: Algorithm to be used when calculating type Ia yields from Chandrasekhar-mass exploders. 0 = DD7 (Iwamoto 1999), 1 = Seitenzahl 2013 3D hydro yields (you must also set Seitenzahl2013_model) 
 | **Parameter input type**: Integer
-| **Default value**: NULL
+| **Default value**: 0
 | **Macros**: ['TYPE_IA_MCH_SUPERNOVA_ALGORITHM_DD2 = 0', 'TYPE_IA_MCH_SUPERNOVA_ALGORITHM_SEITENZAHL2013 = 1', 'TYPE_IA_MCH_SUPERNOVA_ALGORITHM_SEITENZAHL2013_AUTOMATIC = 2']
 
 | **Parameter**: Seitenzahl2013_model
 | **Description**: Which of Seitenzahl et al. 2013's models to use? One of N1,N3,N5,N10,N20,N40,N100L,N100,N100H,N150,N200,N300C,N1600,N1600C,N100_Z0.5,N100_Z0.1,N100_Z0.01 (defaults to N100).
 | **Parameter input type**: String
-| **Default value**: NULL
+| **Default value**: N100
 | **Extra**: N1
 
 | **Parameter**: type_Ia_sub_MCh_supernova_algorithm
 | **Description**: Algorithm to be used when calculating type Ia yields from sub-Chandrasekhar-mass exploders. (Currently unused.)
 | **Parameter input type**: Integer
-| **Default value**: NULL
+| **Default value**: 0
 | **Macros**: ['TYPE_IA_SUB_MCH_SUPERNOVA_ALGORITHM_LIVNE_ARNETT_1995 = 0']
 
 | **Parameter**: max_HeWD_mass
@@ -1221,6 +1295,11 @@ Section: binary
 | **Parameter input type**: Float
 | **Default value**: 0.7
 
+| **Parameter**: merger_mass_loss_fraction
+| **Description**: Fraction of the total mass which is lost when stars merge.
+| **Parameter input type**: Float
+| **Default value**: 0
+
 | **Parameter**: merger_angular_momentum_factor
 | **Description**: When two stars merge the resulting single star retains a fraction of the total system angular momentum (or the critical spin angular momentum, if it is smaller) multiplied by this factor.
 | **Parameter input type**: Float
@@ -1494,7 +1573,6 @@ Section: binary
 | **Description**: Common envelope parameter. The binding energy of the common envelope is G*M*Menv/(lambda*R). Typically this is taken to be 0.5, but if set to -1 binary_c uses the Dewi and Tauris fits instead, -2 uses the formalism of Wang, Jia and Li (2016) and if -3 then a polytropic formalism is used (see also comenv_splitmass).
 | **Parameter input type**: Float(scanf)
 | **Default value**: NULL
-| **Macros**: ['LAMBDA_CE_DEWI_TAURIS = -1', 'LAMBDA_CE_WANG_2016 = -2', 'LAMBDA_CE_POLYTROPE = -3']
 
 | **Parameter**: comenv_splitmass
 | **Description**: When lambda_ce=-2, the envelope binding energy, lambda, is calculated using a polytropic formalism. The comenv_splitmass defines the point, in the units of the core mass, above which material is ejected.
@@ -1690,7 +1768,7 @@ Section: binary
 | **Parameter**: minimum_timestep
 | **Description**: The minimum timestep (Myr).
 | **Parameter input type**: Float
-| **Default value**: 1e-08
+| **Default value**: 1e-06
 
 | **Parameter**: timestep_solver_factor
 | **Description**: Factor applied in timestep_limits, e.g. to prevent X changing too fast by comparing to X/dX/dt, which is usually 1 but can be higher to lengthen timesteps when using an alternative solver.
@@ -1834,9 +1912,19 @@ Section: nucsyn
 | **Parameter**: NeNaMgAl
 | **Description**: Enables NeNaMgAl reaction network. Requires NUCSYN and NUCSYN_HBB.
 | **Parameter input type**: True|False
-| **Default value**: NULL
+| **Default value**: True
 | **Extra**: Ignore
 
+| **Parameter**: nucsyn_network%d
+| **Description**: Usage: --nucsyn_network%d <boolean>. Turn a nuclear network on or off.
+| **Parameter input type**: Boolean(scanf)
+| **Default value**: NULL
+
+| **Parameter**: nucsyn_network_error%d
+| **Description**: Usage: --nucsyn_network_error%d <f>. Threshold error in nuclear network solver for network %d.
+| **Parameter input type**: Float(scanf)
+| **Default value**: NULL
+
 | **Parameter**: nucreacmult%d
 | **Description**: Usage: --nucreacmult%d <f>. Multiply nuclear reaction given by the integer %d (integer) by f (float). 
 | **Parameter input type**: Float(scanf)
@@ -1845,13 +1933,20 @@ Section: nucsyn
 | **Parameter**: nucsyn_metallicity
 | **Description**: This sets the metallicity of the nucleosynthesis algorithms, i.e. the amount (by mass) of matter which is not hydrogen or helium. Usually you'd just set this with the metallicity parameter, but if you want the nucleosynthesis to be outside the range of the stellar evolution algorithm (e.g. Z=0 or Z=0.04) then you need to use nucsyn_metallicity. That said, it's also outside the range of some of the nucleosynthesis algorithms as well, so you have been warned!
 | **Parameter input type**: Float
-| **Default value**: NULL
+| **Default value**: -1
 | **Macros**: ['DEFAULT_TO_METALLICITY = -1']
 
+| **Parameter**: nucsyn_solver
+| **Description**: Choose the solver used in nuclear burning. 0 = KAPS_RENTROP is a Kaps-Rentrop scheme (fast, not great for stiff problems), 1 = LSODA (Adams/BSF switcher), 2 = CVODE library (https://computing.llnl.gov/projects/sundials. Default 0. 
+| **Parameter input type**: Unsigned integer
+| **Default value**: 0
+| **Macros**: ['NUCSYN_SOLVER_KAPS_RENTROP = 0', 'NUCSYN_SOLVER_LSODA = 1', 'NUCSYN_SOLVER_CVODE = 2', 'NUCSYN_SOLVER_NUMBER = 3', 'NUCSYN_SOLVER_KAPS_RENTROP = 0', 'NUCSYN_SOLVER_LSODA = 1', 'NUCSYN_SOLVER_CVODE = 2', 'NUCSYN_SOLVER_NUMBER = 3']
+| **Extra**: 0
+
 | **Parameter**: initial_abundance_mix
 | **Description**: initial abundance mixture: 0=AG89, 1=Karakas 2002, 2=Lodders 2003, 3=Asplund 2005 (not available?), 4=Garcia Berro, 5=Grevesse Noels 1993
 | **Parameter input type**: Unsigned integer
-| **Default value**: NULL
+| **Default value**: 0
 | **Macros**: ['NUCSYN_INIT_ABUND_MIX_AG89 = 0', 'NUCSYN_INIT_ABUND_MIX_KARAKAS2002 = 1', 'NUCSYN_INIT_ABUND_MIX_LODDERS2003 = 2', 'NUCSYN_INIT_ABUND_MIX_ASPLUND2005 = 3', 'NUCSYN_INIT_ABUND_MIX_GARCIABERRO = 4', 'NUCSYN_INIT_ABUND_MIX_GREVESSE_NOELS_1993 = 5', 'NUCSYN_INIT_ABUND_MIX_ASPLUND2009 = 6', 'NUCSYN_INIT_ABUND_MIX_KOBAYASHI2011_ASPLUND2009 = 7', 'NUCSYN_INIT_ABUND_MIX_LODDERS2010 = 8']
 | **Extra**: 0
 
@@ -1876,32 +1971,32 @@ Section: nucsyn
 | **Parameter**: init_abunds_only
 | **Description**: If True, outputs only the initial abundances, then exits.
 | **Parameter input type**: True|False
-| **Default value**: NULL
+| **Default value**: False
 
 | **Parameter**: initial_abunds_only
 | **Description**: If True, outputs only the initial abundances, then exits.
 | **Parameter input type**: True|False
-| **Default value**: NULL
+| **Default value**: False
 
 | **Parameter**: no_thermohaline_mixing
 | **Description**: If True, disables thermohaline mixing.
 | **Parameter input type**: True|False
-| **Default value**: NULL
+| **Default value**: False
 
 | **Parameter**: lithium_GB_post_Heflash
 | **Description**: Sets the lithium abundances after the helium flash. Requires NUCSYN and LITHIUM_TABLES.
 | **Parameter input type**: Float
-| **Default value**: NULL
+| **Default value**: 0
 
 | **Parameter**: lithium_GB_post_1DUP
 | **Description**: Sets the lithium abundance after first dredge up. Requires NUCSYN and LITHIUM_TABLES.
 | **Parameter input type**: Float
-| **Default value**: NULL
+| **Default value**: 0
 
 | **Parameter**: lithium_hbb_multiplier
 | **Description**: Multiplies the lithium abundances on the AGB during HBB (based on Karakas/Fishlock et al models).Requires NUCSYN and LITHIUM_TABLES.
 | **Parameter input type**: Float
-| **Default value**: NULL
+| **Default value**: 1
 
 | **Parameter**: angelou_lithium_decay_function
 | **Description**: Functional form which describes Li7 decay. Requires NUCSYN and NUCSYN_ANGELOU_LITHIUM. Choices are: 0 expoential (see angelou_lithium_decay_time).
@@ -2027,10 +2122,21 @@ Section: nucsyn
 Section: output
 ---------------
 
+| **Parameter**: david_logging_function
+| **Description**: Function to choose which kind of information gets logged Requires DAVID. Choices are: 0= None, >0 for custom logging functions
+| **Parameter input type**: Integer
+| **Default value**: 0
+| **Extra**: Ignore
+
 | **Parameter**: cf_amanda_log
 | **Description**: Enable logging to compare to Amanda's models.
 | **Parameter input type**: True|False
-| **Default value**: NULL
+| **Default value**: False
+
+| **Parameter**: save_pre_events_stardata
+| **Description**: Enable this to save a copy of stardata to stardata->pre_events_stardata just before an event.
+| **Parameter input type**: True|False
+| **Default value**: False
 
 | **Parameter**: disable_end_logging
 | **Description**: Disable the logging that happens at the end of the evolution.
@@ -2060,7 +2166,7 @@ Section: output
 | **Parameter**: legacy_yields
 | **Description**: Turn on ensemble legacy yield output.
 | **Parameter input type**: True|False
-| **Default value**: NULL
+| **Default value**: False
 
 | **Parameter**: ensemble_defer
 | **Description**: Defer ensemble output.
@@ -2068,7 +2174,7 @@ Section: output
 | **Default value**: False
 
 | **Parameter**: ensemble_dt
-| **Description**: When doing ensemble calculations, data is stored and/or output every ensemble_dt Myr. See also ensemble, ensemble_logdt, ensemble_startlogtime.
+| **Description**: When doing ensemble calculations, data are stored and/or output every ensemble_dt Myr. See also ensemble, ensemble_logdt, ensemble_startlogtime.
 | **Parameter input type**: Float
 | **Default value**: 1
 
@@ -2098,19 +2204,29 @@ Section: output
 | **Default value**: False
 
 | **Parameter**: EMP_logg_maximum
-| **Description**: Maximum logg that EMP stars are allowed to have. See Izzard et al 2009. See also CEMP_cfe_minimum, EMP_minimum_age.
+| **Description**: Maximum logg that EMP stars are allowed to have. See Izzard et al 2009. See also CEMP_cfe_minimum, NEMP_nfe_minimum, EMP_minimum_age.
 | **Parameter input type**: Float
-| **Default value**: NULL
+| **Default value**: 4
 
 | **Parameter**: EMP_minimum_age
-| **Description**: Minimum age that EMP stars are required to have. See Izzard et al 2009. See also CEMP_cfe_minimum, EMP_logg_maximum.
+| **Description**: Minimum age that EMP stars are required to have. See Izzard et al 2009. See also CEMP_cfe_minimum, NEMP_nfe_minimum, EMP_logg_maximum.
 | **Parameter input type**: Float
-| **Default value**: NULL
+| **Default value**: 10
+
+| **Parameter**: EMP_feh_maximum
+| **Description**: Maximum [Fe/H] that an EMP stars may have. See Izzard et al 2009. See also CEMP_cfe_minimum, NEMP_nfe_minimum, EMP_logg_maximum, EMP_minimum_age. Default -2.0.
+| **Parameter input type**: Float
+| **Default value**: -2
 
 | **Parameter**: CEMP_cfe_minimum
-| **Description**: Minimum [C/Fe] that CEMP stars are required to have. See Izzard et al 2009. See also EMP_logg_maximum, EMP_minimum_age.
+| **Description**: Minimum [C/Fe] that CEMP stars are required to have. See Izzard et al 2009. See also NEMP_cfe_minimum, EMP_logg_maximum, EMP_minimum_age. Default 0.7.
 | **Parameter input type**: Float
-| **Default value**: NULL
+| **Default value**: 0.7
+
+| **Parameter**: NEMP_cfe_minimum
+| **Description**: Minimum [N/Fe] that NEMP stars are required to have. See Izzard et al 2009, Pols et al. 2012. See also CEMP_cfe_minimum, EMP_logg_maximum, EMP_minimum_age. Default 1.0.
+| **Parameter input type**: Float
+| **Default value**: 1
 
 | **Parameter**: thick_disc_start_age
 | **Description**: Lookback time for the start of the thick disc star formation, e.g. 13e3 Myr. Units = Myr.
@@ -2154,6 +2270,12 @@ Section: output
 | **Default value**: /tmp/c_log.dat
 | **Extra**: 
 
+| **Parameter**: stopfile
+| **Description**: File which, when it exists, will stop the current binary_c repeat run.
+| **Parameter input type**: String
+| **Default value**: 
+| **Extra**: 
+
 | **Parameter**: stardata_dump_filename
 | **Description**: Location of the stardata dump file.
 | **Parameter input type**: String
@@ -2228,12 +2350,6 @@ Section: output
 | **Default value**: NULL
 | **Extra**: Ignore
 
-| **Parameter**: david_logging_function
-| **Description**: Function to choose which kind of information gets logged Requires DAVID. Choices are: 0= None, 1=
-| **Parameter input type**: Integer
-| **Default value**: 0
-| **Extra**: Ignore
-
 Section: input
 --------------
 
@@ -2243,12 +2359,68 @@ Section: input
 | **Default value**: NULL
 | **Extra**: 
 
+| **Parameter**: MINT_data_cleanup
+| **Description**: Activate checks on incoming data to try to account for problems. Will make data-loading slower, but may fix a few things.
+| **Parameter input type**: True|False
+| **Default value**: NULL
+| **Extra**: 
+
 | **Parameter**: MINT_MS_rejuvenation
-| **Description**: Turn on or off (hydrogen) main-sequence rejuvenation. :
+| **Description**: Turn on or off (hydrogen) main-sequence rejuvenation.
+| **Parameter input type**: True|False
+| **Default value**: NULL
+| **Extra**: 
+
+| **Parameter**: MINT_remesh
+| **Description**: Turn on or off MINT's remeshing.
+| **Parameter input type**: True|False
+| **Default value**: NULL
+| **Extra**: 
+
+| **Parameter**: MINT_disable_grid_load_warnings
+| **Description**: Use this to explicitly disable MINT's warnings when loading a grid with, e.g., missing or too much data.
+| **Parameter input type**: True|False
+| **Default value**: NULL
+| **Extra**: 
+
+| **Parameter**: MINT_Kippenhahn
+| **Description**: Turn on or off MINT's Kippenhahn diagrams. If 0, off, if 1, output star 1 (index 0), if 2 output star 2 (index 1). Default 0.
+| **Parameter input type**: Integer
+| **Default value**: NULL
+| **Extra**: 
+
+| **Parameter**: MINT_Kippenhahn_stellar_type
+| **Description**: Stellar type selector for Kippenhahn plots. Set to -1 to ignore, otherwise the stellar type number for which Kippenhahn plot data should be output.
+| **Parameter input type**: Integer
+| **Default value**: NULL
+| **Macros**: ['LOW_MASS_MS = 0', 'MS = 1', 'HG = 2', 'GIANT_BRANCH = 3', 'CHeB = 4', 'EAGB = 5', 'TPAGB = 6', 'HeMS = 7', 'HeHG = 8', 'HeGB = 9', 'HeWD = 10', 'COWD = 11', 'ONeWD = 12', 'NS = 13', 'BH = 14', 'MASSLESS_REMNANT = 15']
+| **Extra**: 
+
+| **Parameter**: MINT_Kippenhahn_companion_stellar_type
+| **Description**: Companion stellar type selector for Kippenhahn plots. Set to -1 to ignore, otherwise the stellar type number for the companion for which Kippenhahn plot data should be output.
+| **Parameter input type**: Integer
+| **Default value**: NULL
+| **Macros**: ['LOW_MASS_MS = 0', 'MS = 1', 'HG = 2', 'GIANT_BRANCH = 3', 'CHeB = 4', 'EAGB = 5', 'TPAGB = 6', 'HeMS = 7', 'HeHG = 8', 'HeGB = 9', 'HeWD = 10', 'COWD = 11', 'ONeWD = 12', 'NS = 13', 'BH = 14', 'MASSLESS_REMNANT = 15']
+| **Extra**: 
+
+| **Parameter**: MINT_nuclear_burning
+| **Description**: Turn on or off MINT's nuclear burning algorithm.
 | **Parameter input type**: True|False
 | **Default value**: NULL
 | **Extra**: 
 
+| **Parameter**: MINT_minimum_shell_mass
+| **Description**: Minimum shell mass in MINT's nuclear burning routines.
+| **Parameter input type**: Float
+| **Default value**: NULL
+| **Extra**: 
+
+| **Parameter**: MINT_maximum_shell_mass
+| **Description**: Maximum shell mass in MINT's nuclear burning routines. :
+| **Parameter input type**: Float
+| **Default value**: NULL
+| **Extra**: 
+
 Section: i/o
 ------------
 
@@ -2358,6 +2530,12 @@ Section: misc
 | **Default value**: NULL
 | **Extra**: Ignore
 
+| **Parameter**: argopts
+| **Description**: Display argument options. Usage: --argopts <argument>.
+| **Parameter input type**: *
+| **Default value**: NULL
+| **Extra**: Ignore
+
 | **Parameter**: help_all
 | **Description**: Display all help pages.
 | **Parameter input type**: *
diff --git a/docs/source/grid_options_descriptions.rst b/docs/source/grid_options_descriptions.rst
index 7c1c75c4f19882a77b78468b995bc5dd9594dac9..0de26933c993e37d8e811fa684ae9cdb3daf7434 100644
--- a/docs/source/grid_options_descriptions.rst
+++ b/docs/source/grid_options_descriptions.rst
@@ -37,12 +37,16 @@ The following options are meant to be changed by the user.
 
 | **log_runtime_systems**: Whether to log the runtime of the systems . Each systems run by the thread is logged to a file and is stored in the tmp_dir. (1 file per thread). Don't use this if you are planning to run alot of systems. This is mostly for debugging and finding systems that take long to run. Integer, default = 0. if value is 1 then the systems are logged
 
+| **max_queue_size**: Maximum size of the queue that is used to feed the processes. Don't make this too big! Default: 1000. Input: int
+
 | **modulo**: No description available yet
 
 | **parse_function**: Function that the user can provide to handle the output the binary_c. This function has to take the arguments (self, output). Its best not to return anything in this function, and just store stuff in the grid_options['results'] dictionary, or just output results to a file
 
 | **repeat**: Factor of how many times a system should be repeated. Consider the evolution splitting binary_c argument for supernovae kick repeating.
 
+| **run_zero_probability_system**: Whether to run the zero probability systems. Default: True. Input: boolean
+
 | **slurm**: Int flag whether to use a slurm type population evolution.
 
 | **source_file_filename**: Variable containing the source file containing lines of binary_c commandline calls. These all have to start with binary_c.
@@ -53,9 +57,90 @@ The following options are meant to be changed by the user.
 
 | **weight**: Weight factor for each system. The calculated probability is mulitplied by this. If the user wants each system to be repeated several times, then this variable should not be changed, rather change the _repeat variable instead, as that handles the reduction in probability per system. This is useful for systems that have a process with some random element in it.
 
+Moe & Distefano sampler options
+-------------------------------
+The following options are meant to be changed by the user.
+
+
+| **Mmin**: Minimum stellar mass
+
+| **multiplicity_model**: 
+multiplicity model (as a function of log10M1)
+
+You can use 'Poisson' which uses the system multiplicty
+given by Moe and maps this to single/binary/triple/quad
+fractions.
+
+Alternatively, 'data' takes the fractions directly
+from the data, but then triples and quadruples are
+combined (and there are NO quadruples).
+
+
+| **multiplicity_modulator**: 
+[single, binary, triple, quadruple]
+
+e.g. [1,0,0,0] for single stars only
+     [0,1,0,0] for binary stars only
+
+defaults to [1,1,0,0] i.e. singles and binaries
+
+
+| **normalize_multiplicities**: 
+'norm': normalize so the whole population is 1.0
+        after implementing the appropriate fractions
+        S/(S+B+T+Q), B/(S+B+T+Q), T/(S+B+T+Q), Q/(S+B+T+Q)
+        given a mix of multiplities, you can either (noting that
+        here (S,B,T,Q) = appropriate modulator * model(S,B,T,Q) )
+        note: if you only set one multiplicity_modulator
+        to 1, and all the others to 0, then normalizing
+        will mean that you effectively have the same number
+        of stars as single, binary, triple or quad (whichever
+        is non-zero) i.e. the multiplicity fraction is ignored.
+        This is probably not useful except for
+        testing purposes or comparing to old grids.
+
+'raw'   : stick to what is predicted, i.e.
+          S/(S+B+T+Q), B/(S+B+T+Q), T/(S+B+T+Q), Q/(S+B+T+Q)
+          without normalization
+          (in which case the total probability < 1.0 unless
+           all you use single, binary, triple and quadruple)
+
+'merge' : e.g. if you only have single and binary,
+          add the triples and quadruples to the binaries, so
+          binaries represent all multiple systems
+          ...
+          *** this is canonical binary population synthesis ***
+
+          It only takes the maximum multiplicity into account,
+          i.e. it doesn't multiply the resulting array by the multiplicity modulator again.
+          This prevents the resulting array to always be 1 if only 1 multiplicity modulator element is nonzero
+
+          Note: if multiplicity_modulator == [1,1,1,1] this
+                option does nothing (equivalent to 'raw').
+
+
+| **q_high_extrapolation_method**: Same as q_low_extrapolation_method
+
+| **q_low_extrapolation_method**: 
+q extrapolation (below 0.15) method
+    none
+    flat
+    linear2
+    plaw2
+    nolowq
+
+
+| **ranges**: 
+
+| **resolutions**: 
+
 Private options
 ---------------
 The following options are not meant to be changed by the user, as these options are used and set internally by the object itself. The description still is provided, but just for documentation purposes.
+
+
+| **_actually_evolve_system**: Whether to actually evolve the systems of just act as if. for testing. used in _process_run_population_grid
+
 | **_binary_c_config_executable**: Full path of the binary_c-config executable. This options is not used in the population object.
 
 | **_binary_c_dir**: Director where binary_c is stored. This options are not really used
@@ -92,11 +177,19 @@ The following options are not meant to be changed by the user, as these options
 
 | **_probtot**: Total probability of the population.
 
+| **_set_ms_grid**: Internal flag whether the M&S grid has been loaded
+
 | **_start_time_evolution**: Variable storing the start timestamp of the population evolution. Set by the object itself.
 
 | **_store_memaddr**: Memory adress of the store object for binary_c.
 
 | **_system_generator**: Function object that contains the system generator function. This can be from a grid, or a source file, or a montecarlo grid.
 
+| **_total_mass_run**: To count the total mass that thread/process has ran
+
+| **_total_probability_weighted_mass_run**: To count the total mass * probability for each system that thread/process has ran
+
 | **_total_starcount**: Variable storing the total amount of systems in the generator. Used and set by the population object.
 
+| **_zero_prob_stars_skipped**: Internal counter to track how many systems are skipped because they have 0 probability
+