diff --git a/.gitignore b/.gitignore
index b82f7e98c..b6122f165 100644
--- a/.gitignore
+++ b/.gitignore
@@ -118,4 +118,6 @@ fort.99
 savefig/
 
 # ignore generated docs files
-_generated/
\ No newline at end of file
+_generated/
+docs/auto_examples/
+sg_execution_times.rst
\ No newline at end of file
diff --git a/changelog.md b/changelog.md
index d6966756d..a01f15ef8 100644
--- a/changelog.md
+++ b/changelog.md
@@ -101,4 +101,14 @@ This release contains _several_ fixes to how CO core masses/remnant masses are h
     - WD chandrasekhar check use ``mc_co`` instead of ``mc``
     - Remnant flags 0-4 inclusive use ``mc_co`` instead of ``mc``
     - Got rid of ``mcx`` in ``assign_remnant`` in favour of clearer ``m_proto`` and ``m_FeNi`` to match the papers
-    - [Very minor] Fryer Rapid was using <= instead of < everywhere
\ No newline at end of file
+    - [Very minor] Fryer Rapid was using <= instead of < everywhere
+
+
+## 3.7.3
+- Additions/changes:
+    - Add new setting ``LBV_flag`` which allows one to turn off LBV winds, use Hurley+2000, or use Belcyznski+2008
+    - Change the default LBV winds to Hurley
+
+- Documentation:
+    - Start new settings gallery in the documentation
+    - Tag settings/options with the version they were added in the docs page and auto link them to release
\ No newline at end of file
diff --git a/docs/_static/cosmic-docs.css b/docs/_static/cosmic-docs.css
index 41cbbd2cc..3989e656a 100644
--- a/docs/_static/cosmic-docs.css
+++ b/docs/_static/cosmic-docs.css
@@ -125,6 +125,16 @@ html[data-theme="dark"] .sd-card {
     max-width: 100%;
 }
 
+.setting-chooser .name-cont {
+    display: flex;
+    align-items: center;
+}
+
+.setting-chooser .name {
+    display: inline-block;
+    margin-right: 0.5rem;
+}
+
 .setting-chooser .description {
     margin-bottom: 0;
 }
@@ -134,6 +144,10 @@ html[data-theme="dark"] .sd-card {
     font-weight: bold;
 }
 
+.setting-chooser .version-added {
+    font-style: italic;
+}
+
 .options {
     overflow: hidden;
 }
@@ -152,6 +166,41 @@ html[data-theme="dark"] .sd-card {
     color: var(--color-brand-primary);
 }
 
+.options .opt-badge-cont {
+    margin-left: 0.5rem;
+    font-style: italic;
+}
+
+/* base badge */
+.badge {
+  display: inline-block;
+  padding: 0.35em 0.65em;
+  font-size: 0.75em;
+  font-weight: 700;
+  line-height: 1;
+  color: #fff;
+  text-align: center;
+  white-space: nowrap;
+  vertical-align: baseline;
+  border-radius: 0.375rem;
+}
+
+/* pill style */
+.badge.rounded-pill {
+  border-radius: 50rem;
+}
+
+.badge-primary    { background-color: #0d6efd; }
+.badge-secondary  { background-color: #6c757d; }
+.badge-success    { background-color: #198754; }
+.badge-danger     { background-color: #dc3545; }
+.badge-warning    { background-color: #ffc107; color: #000; }
+.badge-info       { background-color: #0dcaf0; color: #000; }
+.badge-light      { background-color: #f8f9fa; color: #000; }
+.badge-dark       { background-color: #212529; }
+
+
+
 .setting .col-3 label {
     cursor: pointer;
 }
@@ -469,4 +518,19 @@ select.form-control {
 .btn-toggle.active {
     filter: brightness(1.1);
     font-weight: bold;
+}
+
+/* Centre gallery thumbnail labels and titles */
+.sphx-glr-thumbcontainer[tooltip]:hover::after,
+.sphx-glr-thumbnail-title {
+    text-align: center;
+}
+
+/* vertically align content in the thumbnails */
+.sphx-glr-thumbcontainer {
+    justify-content: center;
+}
+
+.link-white {
+    color: white!important;
 }
\ No newline at end of file
diff --git a/docs/_static/gallery.mplstyle b/docs/_static/gallery.mplstyle
new file mode 100644
index 000000000..d3ec8e44d
--- /dev/null
+++ b/docs/_static/gallery.mplstyle
@@ -0,0 +1,14 @@
+font.family: serif
+text.usetex: False
+figure.figsize: 12, 8
+axes.titlesize: 24
+legend.title_fontsize: 18
+legend.fontsize: 16
+axes.labelsize: 24
+xtick.labelsize: 21
+ytick.labelsize: 21
+axes.linewidth: 1.1
+xtick.major.size: 7
+xtick.minor.size: 4
+ytick.major.size: 7
+ytick.minor.size: 4
\ No newline at end of file
diff --git a/docs/conf.py b/docs/conf.py
index 804e0e22a..8009f9755 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -62,6 +62,7 @@ def setup(app):
     'numpydoc',
     'sphinx_design',
     'sphinx_copybutton',
+    'sphinx_gallery.gen_gallery',
 ]
 
 # -- Extensions ---------------------------------------------------------------
@@ -76,6 +77,15 @@ def setup(app):
 autoclass_content = 'class'
 autodoc_default_flags = ['show-inheritance', 'members', 'inherited-members']
 
+# -- sphinx_gallery -----------------------------
+
+sphinx_gallery_conf = {
+    'examples_dirs': 'settings_examples',   # path to your example scripts
+    'gallery_dirs': 'auto_examples',  # path to where to save gallery generated output
+    'download_all_examples': False,
+    'remove_config_comments': True,
+}
+
 # -- autosummary --------------------------------
 
 autosummary_generate = True
diff --git a/docs/create_settings_html.py b/docs/create_settings_html.py
index 02a782086..9d4f6e50b 100644
--- a/docs/create_settings_html.py
+++ b/docs/create_settings_html.py
@@ -10,6 +10,22 @@
 import bs4
 import json
 import pandas as pd
+from cosmic import __version__ as cosmic_version
+
+# how many versions in the past should we show as badges for when options were added?
+VERSION_CUTOFFS = {
+    "major": 0,
+    "minor": 3,
+    "patch": 1000
+}
+
+def version_is_recent(version_string):
+    version_nums = version_string.split(".")
+    cosmic_nums = cosmic_version.split(".")
+    for i, label in enumerate(["major", "minor", "patch"]):
+        if int(cosmic_nums[i]) - int(version_nums[i]) > VERSION_CUTOFFS[label]:
+            return False
+    return True
 
 # blame BS4 for me calling this soup
 main_soup = """<!-- This file should be created using create_settings_html.py -->
@@ -29,7 +45,7 @@
 settings_template = """<div class="setting">
                 <div class="row align-items-center setting-chooser">
                     <div class="col-9">
-                        <h3 class="name"><code></code></h3>
+                        <div class='name-cont'><h3 class="name"><code></code></h3><span class="version-added"></span></div>
                         <p class="description"></p>
                         <p class="default"></p>
                     </div>
@@ -43,7 +59,7 @@
             </div>"""
 
 # same as above, but for an option for a setting, this will go in the <ul> element
-option_template = """<li><code class="docutils literal notranslate"><span class="pre opt-val"></span></code>: <span class="opt-desc"></span></li>"""
+option_template = """<li><code class="docutils literal notranslate"><span class="pre opt-val"></span></code>: <span class="opt-desc"></span><span class="opt-badge-cont"></span></li>"""
 
 # read the settings file
 with open("../src/cosmic/data/cosmic-settings.json") as f:
@@ -93,6 +109,12 @@
         new_setting.select_one(".description").clear()
         new_setting.select_one(".description").append(bs4.BeautifulSoup(setting["description"],
                                                                         'html.parser'))
+        
+        # add the version added if it's there
+        if "version_added" in setting:
+            if version_is_recent(setting["version_added"]):
+                new_setting.select_one(".version-added").append(bs4.BeautifulSoup(
+                    f"""<span class="badge badge-success version-added"><a class="link-white" href="https://github.com/COSMIC-PopSynth/COSMIC/releases/tag/v{setting['version_added']}">Added in v{setting['version_added']}</a></span>""", 'html.parser'))
 
         # colour the sublinks the same as the border of the group
         new_setting.select_one(".options-expander")["style"] = "color: " + group["docs-colour"] + ";"
@@ -165,6 +187,12 @@
             new_option_expl = bs4.BeautifulSoup(option_template, 'html.parser')
             new_option_expl.select_one(".opt-val").string = str(option["name"])
             new_option_expl.select_one(".opt-desc").append(bs4.BeautifulSoup(option["description"], 'html.parser'))
+
+            if "version_added" in option:
+                if version_is_recent(option["version_added"]):
+                    new_option_expl.select_one(".opt-badge-cont").append(bs4.BeautifulSoup(
+                        f"""<span class="badge badge-success version-added"><a class="link-white" href="https://github.com/COSMIC-PopSynth/COSMIC/releases/tag/v{option['version_added']}">Added in v{option['version_added']}</a></span>""", 'html.parser'
+                    ))
             new_setting.select_one(".options").ul.append(new_option_expl)
 
         # convert the default options to a string and display it
diff --git a/docs/pages/reference_material.rst b/docs/pages/reference_material.rst
index bba77b61a..aff6befd2 100644
--- a/docs/pages/reference_material.rst
+++ b/docs/pages/reference_material.rst
@@ -27,4 +27,15 @@ Configuration and output
 
             output_info.rst
 
-        Descriptions of the output files generated by COSMIC simulations.
\ No newline at end of file
+        Descriptions of the output files generated by COSMIC simulations.
+
+    .. grid-item-card::
+        :link: gallery
+        :link-type: ref
+
+        .. toctree::
+            :maxdepth: 1
+
+            ../auto_examples/index.rst
+
+        Examples of the effects of different configuration options on the evolution of binaries.
\ No newline at end of file
diff --git a/docs/settings_examples/GALLERY_HEADER.rst b/docs/settings_examples/GALLERY_HEADER.rst
new file mode 100644
index 000000000..34880775a
--- /dev/null
+++ b/docs/settings_examples/GALLERY_HEADER.rst
@@ -0,0 +1,9 @@
+Impact of settings on evolution
+===============================
+
+This gallery shows examples of how different settings affect the evolution of binaries.
+
+.. admonition:: Under construction
+    :class: warning
+
+    This page is still being developed, but will eventually contain examples of how all settings affect the evolution of binaries. For now it's just a few.
\ No newline at end of file
diff --git a/docs/settings_examples/plot_LBV_flag.py b/docs/settings_examples/plot_LBV_flag.py
new file mode 100644
index 000000000..8b37b8411
--- /dev/null
+++ b/docs/settings_examples/plot_LBV_flag.py
@@ -0,0 +1,104 @@
+"""
+``LBV_flag``
+============
+
+This example shows the effect of the ``LBV_flag`` on the evolution of massive stars.
+
+The ``LBV_flag`` controls which prescription to use for LBV-like mass loss for stars that are above the Humphreys-Davidson limit. The plots below show the evolution of a few single stars and highlight the LBV regime in which we'd expect to see differences with this flag.
+"""
+
+#----------------------------------------------------------------------------------
+#----------------------------------------------------------------------------------
+## You'll want to edit this part locally to use your own BSEDict and style sheet!
+
+import sys
+sys.path.append("..")
+import generate_default_bsedict
+BSEDict = generate_default_bsedict.get_default_BSE_settings(to_python=True)
+
+import matplotlib.pyplot as plt
+plt.style.use("../_static/gallery.mplstyle")
+#----------------------------------------------------------------------------------
+#----------------------------------------------------------------------------------
+
+import numpy as np
+import astropy.units as u
+import astropy.constants as consts
+
+from cosmic.sample import InitialBinaryTable
+from cosmic.evolve import Evolve
+from cosmic.output import kstar_translator
+
+def LBV_limit(T_eff):
+    """ Compute the luminosity of the LBV limit given a certain temperature """
+    return np.maximum(np.sqrt(1e10 * u.Rsun**2 * u.Lsun * (4 * np.pi * consts.sigma_sb * (T_eff * u.K)**4)).to(u.Lsun), np.repeat(6e5 * u.Lsun, len(T_eff)))
+
+
+# create a small grid of single stars
+masses = [30, 33, 36, 40, 45, 50]
+grid = InitialBinaryTable.InitialBinaries(
+    m1=masses,
+    m2=np.zeros(len(masses)),
+    porb=np.ones(len(masses))*-1.0,
+    ecc=np.ones(len(masses))*0.0,
+    tphysf=np.ones(len(masses))*13700.0,
+    kstar1=np.ones(len(masses)),
+    kstar2=np.ones(len(masses)),
+    metallicity=np.ones(len(masses))*0.014*0.01
+)
+
+# loop over different flag choices
+for LBV_flag, label in zip([0, 1, 2], ["No LBV mass loss", "LBV mass loss following Hurley+2000",
+                                       "LBV mass loss following Belczynski+2008"]):
+    # evolve with updated BSEDict
+    BSEDict["LBV_flag"] = LBV_flag
+    bpp, bcm, initC, kick_info = Evolve.evolve(
+        initialbinarytable=grid, BSEDict=BSEDict,
+        dtp=0
+    )
+
+    # plot a HR diagram of the evolution of these stars, colour coded by kstar
+    fig, ax = plt.subplots(figsize=(10, 16))
+
+    for i, bin_num in enumerate(np.unique(bcm["bin_num"])):
+        log_teff = np.log10(bcm[bcm["bin_num"] == bin_num]["teff_1"].values)
+        log_lum = np.log10(bcm[bcm["bin_num"] == bin_num]["lum_1"].values)
+        kstar = bcm[bcm["bin_num"] == bin_num]["kstar_1"].values
+
+        mask = kstar < 10
+
+        ax.scatter(log_teff[mask], log_lum[mask], c=[kstar_translator[k]["colour"] for k in kstar[mask]],
+                   s=5)
+        ax.plot(log_teff[mask], log_lum[mask], color='grey', alpha=0.5, lw=0.5, zorder=-1)
+
+        # annotate the first point with the initial mass
+        ax.annotate(f'{bcm[bcm["bin_num"] == bin_num]["mass_1"].values[0]:.0f} M$_\odot$', (log_teff[0], log_lum[0]),
+                    fontsize=14, ha="right", color='black', va='top')
+        
+    # annotate the various kstar types
+    for k in bcm["kstar_1"].unique():
+        if k < 10:
+            ax.scatter([], [], color=kstar_translator[k]["colour"], label=kstar_translator[k]["short"], s=50)
+    ax.legend()
+
+    # make some fairly reasonable temperature range
+    T_eff_range = np.logspace(3.5, 5.4, 1000)
+    
+    # plot a line at the HD limit and fill the area
+    ax.plot(np.log10(T_eff_range), np.log10(LBV_limit(T_eff_range).value), color="grey", linestyle="dotted")
+    ax.fill_between(np.log10(T_eff_range), np.log10(LBV_limit(T_eff_range).value), 9, color="black", lw=2, alpha=0.03)
+    
+    # annotate the regime with a custom location
+    ax.annotate("LBV Regime", xy=(0.77, 0.93), xycoords="axes fraction", ha="center", va="center",
+                fontsize=20, color="grey", zorder=10)
+
+    ax.set(
+        xlabel=r"$\log_{10}(T_{\rm eff} \, [\rm K])$",
+        ylabel=r"$\log_{10}(L \, [\rm L_{\odot}])$",
+        title=f"LBV_flag = {LBV_flag}\n{label}",
+        ylim=(5, 6.2)
+    )
+
+    ax.invert_xaxis()
+    
+    plt.show()
\ No newline at end of file
diff --git a/docs/settings_examples/plot_pisn.py b/docs/settings_examples/plot_pisn.py
new file mode 100644
index 000000000..93da22815
--- /dev/null
+++ b/docs/settings_examples/plot_pisn.py
@@ -0,0 +1,109 @@
+"""
+``pisn``
+============
+
+This example shows the effect of the ``pisn`` on the remnant masses for stars with large CO cores.
+
+The ``pisn`` controls which prescription to use for pair-instability supernovae (PISN) for stars that are above a certain CO core mass.
+The plots below show the remnant masses for a grid of single stars with different choices of the ``pisn`` flag at a metallicity of :math:`Z = 0.01 Z_{\odot}`.
+"""
+
+# sphinx_gallery_thumbnail_number = 3
+
+#----------------------------------------------------------------------------------
+#----------------------------------------------------------------------------------
+## You'll want to edit this part locally to use your own BSEDict and style sheet!
+
+import sys
+sys.path.append("..")
+import generate_default_bsedict
+BSEDict = generate_default_bsedict.get_default_BSE_settings(to_python=True)
+
+import matplotlib.pyplot as plt
+plt.style.use("../_static/gallery.mplstyle")
+#----------------------------------------------------------------------------------
+#----------------------------------------------------------------------------------
+
+import numpy as np
+import astropy.units as u
+import astropy.constants as consts
+
+from cosmic.sample import InitialBinaryTable
+from cosmic.evolve import Evolve
+from cosmic.output import kstar_translator
+
+
+pisn_flags = [0, -1, -2, -3, -4, 45]
+labels = ["No PISN", "PISN following Spera & Mapelli 2017", "PISN following Marchant+2018",
+          "PISN following Woolsey+2019", "PISN following\nRenzo+2022/Hendriks+2023",
+          "Assume PPISN from 45-65 Msun He Cores\nPISN above that"]
+
+bpps = {}
+initCs = {}
+
+for pisn_flag, label in zip(pisn_flags, labels):
+    n_grid = 250
+
+    binary_grid = InitialBinaryTable.InitialBinaries(m1=np.geomspace(50, 150, n_grid),
+                                                    m2=np.zeros(n_grid),
+                                                    porb=np.ones(n_grid)*-1.0,
+                                                    ecc=np.ones(n_grid)*0.0,
+                                                    tphysf=np.ones(n_grid)*13700.0,
+                                                    kstar1=np.ones(n_grid),
+                                                    kstar2=np.ones(n_grid),
+                                                    metallicity=np.ones(n_grid)*0.014*0.01)
+
+
+    BSEDict['pisn'] = pisn_flag
+    bpp, _, initC, _ = Evolve.evolve(
+        initialbinarytable=binary_grid, BSEDict=BSEDict,
+        bpp_columns=['tphys', 'mass_1', 'kstar_1', 'massc_co_layer_1',
+                        'massc_he_layer_1', 'evol_type', "SN_1", "kstar_2"],
+        nproc=16
+    )
+    final_bpp = bpp.drop_duplicates(subset="bin_num", keep="last")
+
+    fig, ax = plt.subplots(figsize=(8, 5), constrained_layout=True)
+
+    sn_rows = bpp[bpp["evol_type"] == 15]
+    co_core_mass = sn_rows["massc_co_layer_1"].values
+    total_core_mass = co_core_mass + sn_rows["massc_he_layer_1"].values
+    remnant_mass = final_bpp["mass_1"].values
+    kstars = sn_rows["kstar_1"].values
+
+    ax.scatter(total_core_mass, remnant_mass, c=[kstar_translator[k]["colour"] for k in kstars])
+    ax.plot(total_core_mass, total_core_mass - 0.5, color='k', ls=':')
+
+    ax.set(
+        xlabel=r"Total core mass, $M_{\rm c, tot} \, [\rm M_{\odot}]$",
+        ylabel=r"Remnant mass, $M_r \, [\rm M_{\odot}]$",
+        xlim=(15, 62),
+        ylim=(15, 62),
+        title=label
+    )
+
+    ax.yaxis.set_major_locator(plt.MultipleLocator(10))
+    ax.yaxis.set_minor_locator(plt.MultipleLocator(5))
+
+    ax.xaxis.set_major_locator(plt.MultipleLocator(10))
+    ax.xaxis.set_minor_locator(plt.MultipleLocator(2))
+
+    if pisn_flag == 0:
+        pisn_turn_on = -100
+    elif pisn_flag == -1:
+        pisn_turn_on = 32
+    elif pisn_flag == -2:
+        pisn_turn_on = 32
+    elif pisn_flag == -3:
+        pisn_turn_on = 29.53
+    elif pisn_flag == -4:
+        idx = np.argmin(np.abs(co_core_mass - 38))
+        pisn_turn_on = total_core_mass[idx]
+    else:
+        pisn_turn_on = pisn_flag
+
+    ax.axvline(pisn_turn_on, color='k', ls='--')
+    ax.annotate("(P)PISN turns on", xy=(pisn_turn_on, ax.get_ylim()[-1] - 2), rotation=90,
+                fontsize=10, bbox=dict(boxstyle="round", fc="w"), ha='center', va='top')
+
+    plt.show()
\ No newline at end of file
diff --git a/requirements.txt b/requirements.txt
index 2c7bf7fe0..8317cf3bb 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -24,4 +24,5 @@ wheel
 furo
 sphinx-design
 sphinx-copybutton
+sphinx-gallery
 meson
diff --git a/src/cosmic/_version.py b/src/cosmic/_version.py
index 3bdbbae2e..221485aea 100644
--- a/src/cosmic/_version.py
+++ b/src/cosmic/_version.py
@@ -1 +1 @@
-__version__ = "3.7.2"
+__version__ = "3.7.3"
diff --git a/src/cosmic/data/cosmic-settings.json b/src/cosmic/data/cosmic-settings.json
index 4ca0eec6b..8804a45ff 100644
--- a/src/cosmic/data/cosmic-settings.json
+++ b/src/cosmic/data/cosmic-settings.json
@@ -544,7 +544,8 @@
                 "options": [
                     {
                         "name": -1,
-                        "description": "No wind mass loss"
+                        "description": "No wind mass loss",
+                        "version_added": "3.7.2"
                     },
                     {
                         "name": 0,
@@ -688,6 +689,28 @@
                     }
                 ]
             },
+            {
+                "name": "LBV_flag",
+                "description": "Selects the model for Luminous Blue Variable (LBV) mass loss",
+                "type": "dropdown",
+                "options-preface": "",
+                "version_added": "3.7.3",
+                "options": [
+                    {
+                        "name": 0,
+                        "description": "No LBV mass loss"
+                    },
+                    {
+                        "name": 1,
+                        "description": "LBV mass loss following <a href='https://ui.adsabs.harvard.edu/abs/2000MNRAS.315..543H/abstract'>Hurley+2000</a>",
+                        "default": true
+                    },
+                    {
+                        "name": 2,
+                        "description": "LBV mass loss following <a href='https://ui.adsabs.harvard.edu/abs/2008ApJS..174..223B/abstract'>Belczynski+2008</a>"
+                    }
+                ]
+            },
             {
                 "settings-section": "Common-envelope",
                 "settings-section-description": "Note: there are cases where a common envelope is forced regardless of the critical mass ratio for unstable mass transfer. In the following cases, a common envelope occurs regardless of the choices below:<ul><li><b>contact</b>: the stellar radii go into contact (common for similar ZAMS systems)</li><li><b>periapse contact</b>: the periapse distance is smaller than either of the stellar radii (common for highly eccentric systems)</li><li><b>core Roche overflow</b>: either of the stellar radii overflow their component's Roche radius (in this case, mass transfer from the convective core is always dynamically unstable)</li></ul>",
@@ -1130,6 +1153,7 @@
                 "name": "mm_mu_ns",
                 "description": "Sets the NS kick scaling prefactor for the Mandel & Mueller 2020 kick prescription (only relevant if kickflag = 6)",
                 "type": "number",
+                "version_added": "3.7.0",
                 "options-preface": "",
                 "options": [
                     {
@@ -1147,6 +1171,7 @@
                 "name": "mm_mu_bh",
                 "description": "Sets the BH kick scaling prefactor for the Mandel & Mueller 2020 kick prescription (only relevant if kickflag = 6)",
                 "type": "number",
+                "version_added": "3.7.0",
                 "options-preface": "",
                 "options": [
                     {
@@ -1190,11 +1215,13 @@
                     },
                     {
                         "name": 5,
-                        "description": "follows the prescription from <a class='reference external' href='https://scixplorer.org/abs/2020MNRAS.499.3214M/abstract'>Mandel &amp; Mueller 2020</a>. The method is stochastic. NOTE: We recommend that this prescription be used in conjunction with the Mandel &amp; Mueller 2020 kick prescription (kickflag = 6) and with mxns set to 2.0 Msun and rembar_massloss set to 0.1 Msun."
+                        "description": "follows the prescription from <a class='reference external' href='https://scixplorer.org/abs/2020MNRAS.499.3214M/abstract'>Mandel &amp; Mueller 2020</a>. The method is stochastic. NOTE: We recommend that this prescription be used in conjunction with the Mandel &amp; Mueller 2020 kick prescription (kickflag = 6) and with mxns set to 2.0 Msun and rembar_massloss set to 0.1 Msun.",
+                        "version_added": "3.7.0"
                     },
                     {
                         "name": 6,
-                        "description": "follows the prescription from <a class='reference external' href='https://scixplorer.org/abs/2025A%26A...700A..20M/abstract'>Maltsev+2025</a>, with additions from <a class='reference external' href='https://scixplorer.org/abs/2025arXiv251007573W/abstract'>Willcox+2025</a> for the definition of BH masses."
+                        "description": "follows the prescription from <a class='reference external' href='https://scixplorer.org/abs/2025A%26A...700A..20M/abstract'>Maltsev+2025</a>, with additions from <a class='reference external' href='https://scixplorer.org/abs/2025arXiv251007573W/abstract'>Willcox+2025</a> for the definition of BH masses.",
+                        "version_added": "3.7.0"
                     }
                 ]
             },
@@ -1202,6 +1229,7 @@
                 "name": "fryer_mass_limit",
                 "description": "Sets the mass limit to use when applying a remnant mass prescription from Fryer+12 (remnantflag=3 or 4).",
                 "type": "dropdown",
+                "version_added": "3.7.2",
                 "options-preface": "",
                 "options": [
                     {
@@ -1274,6 +1302,7 @@
                 "name": "maltsev_mode",
                 "description": "Sets the mode to use Maltsev remnant mass prescription (only relevant if remnantflag = 6)",
                 "type": "number",
+                "version_added": "3.7.0",
                 "options-preface": "",
                 "options": [
                     {
@@ -1295,6 +1324,7 @@
                 "name": "maltsev_fallback",
                 "description": "Sets the fallback fraction for partial fallback in the Maltsev remnant mass prescription (only relevant if remnantflag = 6). This fallback fraction, \\(f_{\\rm fallback}\\), is used to calculate the remnant mass by calculating \\(M_{\\rm rem} = M_{\\rm NS} + f_{\\rm fallback} (M_{\\rm core, total} - M_{\\rm NS})\\), where this prescription assumes that \\(M_{\\rm NS} = 1.4 \\, \\mathrm{M}_\\odot\\).",
                 "type": "number",
+                "version_added": "3.7.0",
                 "options-preface": "",
                 "options": [
                     {
@@ -1312,6 +1342,7 @@
                 "name": "maltsev_pf_prob",
                 "description": "Sets the probability of partial fallback leading to BH formation occurring in the Maltsev remnant mass prescription when the CO core mass is in the partial fallback regime of \\(M_2 \\le M_{\\rm CO} \\le M_3\\) (only relevant if remnantflag = 6)",
                 "type": "number",
+                "version_added": "3.7.0",
                 "options-preface": "",
                 "options": [
                     {
@@ -1350,7 +1381,8 @@
                     },
                     {
                         "name": -4,
-                        "description": "Uses the top-down approach from <a class='reference external' href='https://ui.adsabs.harvard.edu/abs/2022RNAAS...6...25R/abstract'>Renzo+2022</a>, with adjustable parameters from <a class='reference external' href='https://scixplorer.org/abs/2023MNRAS.526.4130H/abstract'>Hendriks+2023</a>"
+                        "description": "Uses the top-down approach from <a class='reference external' href='https://ui.adsabs.harvard.edu/abs/2022RNAAS...6...25R/abstract'>Renzo+2022</a>, with adjustable parameters from <a class='reference external' href='https://scixplorer.org/abs/2023MNRAS.526.4130H/abstract'>Hendriks+2023</a>",
+                        "version_added": "3.7.2"
                     },
                     {
                         "name": "positive values",
@@ -1362,6 +1394,7 @@
                 "name": "ppi_co_shift",
                 "description": "Shifts the CO core mass thresholds for pulsational pair instability and pair instability SNe by a constant amount in COSMIC's implementation of the Renzo+2022/Hendriks+2023 prescription (only relevant if pisn = -4). Default threshold is 38 Msun.",
                 "type": "number",
+                "version_added": "3.7.2",
                 "options-preface": "Note the default threshold is 38 Msun.",
                 "options": [
                     {
@@ -1383,6 +1416,7 @@
                 "name": "ppi_extra_ml",
                 "description": "Sets the amount of extra mass loss in Msun that is applied to stars that go through pulsational pair instability in COSMIC's implementation of the Renzo+2022/Hendriks+2023 prescription (only relevant if pisn = -4). Default is 0 Msun.",
                 "type": "number",
+                "version_added": "3.7.2",
                 "options-preface": "",
                 "options": [
                     {
diff --git a/src/cosmic/evolve.py b/src/cosmic/evolve.py
index ccf7e780e..a27a61fbd 100644
--- a/src/cosmic/evolve.py
+++ b/src/cosmic/evolve.py
@@ -107,7 +107,7 @@
                                   'beta', 'xi', 'acc2', 'epsnov',
                                   'eddfac', 'gamma', 'don_lim', 'acc_lim',
                                   'bdecayfac', 'bconst', 'ck',
-                                  'windflag', 'qcflag', 'eddlimflag',
+                                  'windflag', 'qcflag', 'eddlimflag', 'LBV_flag',
                                   'fprimc_array', 'dtp', 'randomseed',
                                   'bhspinflag', 'bhspinmag', 'rejuv_fac', 'rejuvflag', 'htpmb',
                                   'ST_cr', 'ST_tide', 'rembar_massloss', 'zsun', 'kickflag']
@@ -572,6 +572,7 @@ def _evolve_single_system(f):
         _evolvebin.windvars.epsnov = f["epsnov"]
         _evolvebin.windvars.eddfac = f["eddfac"]
         _evolvebin.windvars.gamma = f["gamma"]
+        _evolvebin.windvars.lbv_flag = f["LBV_flag"]
         _evolvebin.flags.bdecayfac = f["bdecayfac"]
         _evolvebin.magvars.bconst = f["bconst"]
         _evolvebin.magvars.ck = f["ck"]
diff --git a/src/cosmic/src/const_bse.h b/src/cosmic/src/const_bse.h
index 1c411ff8e..997e97fed 100644
--- a/src/cosmic/src/const_bse.h
+++ b/src/cosmic/src/const_bse.h
@@ -26,8 +26,9 @@
       COMMON /METVARS/ zsun
       REAL*8 neta,bwind,hewind,beta,xi,acc2,epsnov
       REAL*8 eddfac,gamma
+      INTEGER LBV_flag
       COMMON /WINDVARS/ neta,bwind,hewind,beta,xi,acc2,epsnov,
-     &                  eddfac,gamma
+     &                  eddfac,gamma,LBV_flag
       REAL*8 alpha1,lambdaf
       REAL*8 qcrit_array(16)
       COMMON /CEVARS/ qcrit_array,alpha1,lambdaf
diff --git a/src/cosmic/src/hrdiag.f b/src/cosmic/src/hrdiag.f
index ebba5922c..08eddb509 100644
--- a/src/cosmic/src/hrdiag.f
+++ b/src/cosmic/src/hrdiag.f
@@ -54,7 +54,7 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars,
       real*8 rzamsf,rtmsf,ralphf,rbetaf,rgammf,rhookf
       real*8 rgbf,rminf,ragbf,rzahbf,rzhef,rhehgf,rhegbf,rpertf
       real*8 mctmsf,mcgbtf,mcgbf,mcheif,mcagbf,lzahbf
-*      real*8 mrem
+      logical stripped_during_hrdiag
       external thookf,tblf
       external lalphf,lbetaf,lnetaf,lhookf,lgbtf,lmcgbf,lzhef,lpertf
       external rzamsf,rtmsf,ralphf,rbetaf,rgammf,rhookf
@@ -81,6 +81,8 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars,
 *
 * Make evolutionary changes to stars that have not reached KW > 5.
 *
+      ! track whether a star stripped during hrdiag
+      stripped_during_hrdiag = .false.
       mch = 1.44d0 !set here owing to AIC ECSN model.
 *
       mass0 = mass
@@ -202,8 +204,8 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars,
                   kw = 7
                   CALL star(kw,mass,mt,tm,tn,tscls,lums,GB,zpars)
 
-                  ! return so that bpp logs the stellar type change
-                  return
+                  ! flag to return so that bpp logs the stellar type change
+                  stripped_during_hrdiag = .true.
                else
 *
 * Zero-age helium white dwarf.
@@ -276,8 +278,8 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars,
                kw = 7
                CALL star(kw,mass,mt,tm,tn,tscls,lums,GB,zpars)
 
-               ! return so that bpp logs the stellar type change
-               return
+               ! flag to return so that bpp logs the stellar type change
+               stripped_during_hrdiag = .true.
             else
 *
 * Zero-age helium white dwarf.
@@ -412,8 +414,8 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars,
             CALL star(kw,mass,mt,tm,tn,tscls,lums,GB,zpars)
             aj = xx*tm
 
-            ! return so that bpp logs the stellar type change
-            return
+            ! flag to return so that bpp logs the stellar type change
+            stripped_during_hrdiag = .true.
          else
             kw = 4
          endif
@@ -457,8 +459,8 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars,
                endif
                aj = MAX(aj,tm)
 
-               ! return so that bpp logs the stellar type change
-               return
+               ! flag to return so that bpp logs the stellar type change
+               stripped_during_hrdiag = .true.
             else
                kw = 5
             endif
@@ -546,6 +548,10 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars,
          mc_he(kidx) = mt - mc
          mc_co(kidx) = mc
 
+         if(stripped_during_hrdiag)then
+            return
+         endif
+
          ! core mass for a He star to become a COWD (Hurley 2000, Eq. 89)
          mtc = MIN(mt,1.45d0*mt-0.31d0)
 
diff --git a/src/cosmic/src/mlwind.f b/src/cosmic/src/mlwind.f
index 0a7b4fe09..df1f16449 100644
--- a/src/cosmic/src/mlwind.f
+++ b/src/cosmic/src/mlwind.f
@@ -2,11 +2,11 @@
       real*8 FUNCTION mlwind(kw,lum,r,mt,mc,rl,z)
       IMPLICIT NONE
       INCLUDE 'const_bse.h'
-      integer kw,testflag
+      integer kw
       real*8 lum,r,mt,mc,rl,z,teff,alpha
       real*8 dml,dms,dmt,p0,x,mew,lum0,kap
-      real*8 MLalpha
-      external MLalpha
+      real*8 MLalpha,LBV_winds
+      external MLalpha,LBV_winds
       parameter(lum0=7.0d+04,kap=-0.5d0)
 *
 *      windflag = 0 !BSE=0, startrack08=1, vink=2, vink+LBV for all
@@ -56,11 +56,7 @@ real*8 FUNCTION mlwind(kw,lum,r,mt,mc,rl,z)
                   dms = MAX(dml,dms)
                endif
 * LBV-like mass loss beyond the Humphreys-Davidson limit.
-               x = 1.0d-5*r*sqrt(lum)
-               if(lum.gt.6.0d+05.and.x.gt.1.d0)then
-                  dml = 0.1d0*(x-1.d0)**3*(lum/6.0d+05-1.d0)
-                  dms = dms + dml
-               endif
+               dms = dms + LBV_winds(lum,r,mt,kw,z)
             endif
          endif
 *
@@ -115,11 +111,7 @@ real*8 FUNCTION mlwind(kw,lum,r,mt,mc,rl,z)
                   dms = MAX(dml,dms)
                endif
 * LBV-like mass loss beyond the Humphreys-Davidson limit.
-               x = 1.0d-5*r*sqrt(lum)
-               if(lum.gt.6.0d+05.and.x.gt.1.d0)then
-                  dml = 0.1d0*(x-1.d0)**3*(lum/6.0d+05-1.d0)
-                  dms = dms + dml
-               endif
+               dms = dms + LBV_winds(lum,r,mt,kw,z)
             endif
          endif
 *
@@ -146,7 +138,6 @@ real*8 FUNCTION mlwind(kw,lum,r,mt,mc,rl,z)
             dms = 9.6d-15*x*(r**0.81d0)*(lum**1.24d0)*(mt**0.16d0)
             alpha = 0.5d0
             dms = dms*(z/zsun)**(alpha)
-            testflag = 1
          endif
          if(kw.ge.2.and.kw.le.6)then
 * 'Reimers' mass loss
@@ -175,7 +166,6 @@ real*8 FUNCTION mlwind(kw,lum,r,mt,mc,rl,z)
      &            1.339d0*LOG10(mt/30.d0) - 1.601d0*LOG10(1.3d0/2.d0) +
      &            alpha*LOG10(z/zsun) + 1.07d0*LOG10(teff/2.0d+04)
             dms = 10.d0**dms
-            testflag = 2
          elseif(teff.gt.25000.)then
 *        Although Vink et al. formulae  are only defined until Teff=50000K,
 *        we follow the Dutch prescription of MESA, and extend to higher Teff
@@ -184,27 +174,21 @@ real*8 FUNCTION mlwind(kw,lum,r,mt,mc,rl,z)
      &            alpha*LOG10(z/zsun) +0.933d0*LOG10(teff/4.0d+04) -
      &            10.92d0*(LOG10(teff/4.0d+04)**2)
        dms = 10.d0**dms
-       testflag = 2
          endif
 
          if((windflag.eq.3.or.kw.ge.2).and.kw.le.6)then
 * LBV-like mass loss beyond the Humphreys-Davidson limit.
 * Optional flag (windflag=3) to use for every non-degenerate star
 * past the limit, rather than just for giant, evolved stars
-            x = 1.0d-5*r*sqrt(lum)
-            if(lum.gt.6.0d+05.and.x.gt.1.d0)then
-               if(eddlimflag.eq.0) alpha = 0.d0
-               if(eddlimflag.eq.1) alpha = MLalpha(mt,lum,kw)
-               dms = 1.5d0*1.0d-04*((z/zsun)**alpha)
-               testflag = 3
-            endif
+
+* TW: This previously overwrote the other mass loss, changed it to a sum
+            dms = dms + LBV_winds(lum,r,mt,kw,z)
          elseif(kw.ge.7.and.kw.le.9)then !WR (naked helium stars)
 * If naked helium use Hamann & Koesterke (1998) WR winds reduced by factor of
 * 10 (Yoon & Langer 2005), with Vink & de Koter (2005) metallicity dependence
             if(eddlimflag.eq.0) alpha = 0.86d0
             if(eddlimflag.eq.1) alpha = MLalpha(mt,lum,kw)
             dms = 1.0d-13*(lum**1.5d0)*((z/zsun)**alpha)
-            testflag = 4
          endif
 *
          mlwind = dms
@@ -233,11 +217,7 @@ real*8 FUNCTION mlwind(kw,lum,r,mt,mc,rl,z)
                   dms = MAX(dml,dms)
                endif
 * LBV-like mass loss beyond the Humphreys-Davidson limit.
-               x = 1.0d-5*r*sqrt(lum)
-               if(lum.gt.6.0d+05.and.x.gt.1.d0)then
-                  dml = 0.1d0*(x-1.d0)**3*(lum/6.0d+05-1.d0)
-                  dms = dms + dml
-               endif
+               dms = dms + LBV_winds(lum,r,mt,kw,z)
             endif
          endif
 *
@@ -248,3 +228,41 @@ real*8 FUNCTION mlwind(kw,lum,r,mt,mc,rl,z)
       return
       end
 ***
+
+
+      real*8 FUNCTION LBV_winds(lum,r,mt,kw,z)
+*        Calculate mass loss from LBV-like winds
+      IMPLICIT NONE
+      INCLUDE 'const_bse.h'
+      real*8 lum,r,mt,z,alpha,x
+      integer kw
+      real*8 MLalpha
+      external MLalpha
+
+      x = 1.0d-5*r*sqrt(lum)
+      alpha = 0.d0
+
+      ! if the star is beyond the Humphreys-Davidson limit, apply LBVs
+      if(lum.gt.6.0d+05.and.x.gt.1.d0.and.LBV_flag.ne.0)then
+         if(LBV_flag.eq.1)then
+            ! use Hurley+2000 LBV-like mass loss (Section 7.1, equation is unnumbered oof)
+            LBV_winds = 0.1d0*(x-1.d0)**3*(lum/6.0d+05-1.d0)
+         elseif(LBV_flag.eq.2)then
+            ! use StarTrack (Belczynski+08) LBV-like mass loss
+            ! adjust metallicity dependence near eddington limit if eddlimflag is set
+            if(eddlimflag.eq.1) alpha = MLalpha(mt,lum,kw)
+
+            ! use constant LBV-like mass loss
+            LBV_winds = 1.5d0*1.0d-04*((z/zsun)**alpha)
+         else
+            ! this should never happen, throw error
+            write(*,*) "Error: LBV_flag must be 0, 1, or 2. Exiting."
+            stop
+         endif
+      else
+         ! otherwise no mass loss
+         LBV_winds = 0.d0
+      endif
+
+      return
+      end
diff --git a/src/cosmic/tests/data/Params.ini b/src/cosmic/tests/data/Params.ini
index b4a11f9c2..6fa2996e2 100644
--- a/src/cosmic/tests/data/Params.ini
+++ b/src/cosmic/tests/data/Params.ini
@@ -147,6 +147,8 @@ xi=0.5
 ; default=1.5
 acc2=1.5
 
+LBV_flag=2
+
 
 
 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
diff --git a/src/cosmic/tests/data/initial_conditions_for_testing.hdf5 b/src/cosmic/tests/data/initial_conditions_for_testing.hdf5
index adf91ea81..802853af5 100644
Binary files a/src/cosmic/tests/data/initial_conditions_for_testing.hdf5 and b/src/cosmic/tests/data/initial_conditions_for_testing.hdf5 differ
diff --git a/src/cosmic/tests/data/kick_initial_conditions.h5 b/src/cosmic/tests/data/kick_initial_conditions.h5
index 3794e1791..8e6e428b7 100644
Binary files a/src/cosmic/tests/data/kick_initial_conditions.h5 and b/src/cosmic/tests/data/kick_initial_conditions.h5 differ