diff --git a/doc/papers/Caughlan_Fowler_1988.pdf b/doc/papers/Caughlan_Fowler_1988.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..c6ac1cd199a1109b1f49d30c623edffb06a124ed
Binary files /dev/null and b/doc/papers/Caughlan_Fowler_1988.pdf differ
diff --git a/doc/papers/Fowler_Caughlan_Zimmerman_1967.pdf b/doc/papers/Fowler_Caughlan_Zimmerman_1967.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..f633fde0b2173519b7b26511786ee084f7137775
Binary files /dev/null and b/doc/papers/Fowler_Caughlan_Zimmerman_1967.pdf differ
diff --git a/doc/papers/NACRE.pdf b/doc/papers/NACRE.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..2328e57e8a4eb743c6179c729bc7bb1b7e63eb16
Binary files /dev/null and b/doc/papers/NACRE.pdf differ
diff --git a/src/MINT/MINT_CHeB_TA.c b/src/MINT/MINT_CHeB_TA.c
new file mode 100644
index 0000000000000000000000000000000000000000..3f1c97304da7b876b5641e75d451d2f6638ecd3f
--- /dev/null
+++ b/src/MINT/MINT_CHeB_TA.c
@@ -0,0 +1,46 @@
+#include "../binary_c.h"
+
+No_empty_translation_unit_warning;
+
+#ifdef MINT
+
+void MINT_CHeB_TA(struct stardata_t * const stardata,
+                  const double mass,
+                  double * const CHeB_lifetime,
+                  double * const CHeB_TA_luminosity)
+{
+    /*
+     * Estimate CHeB terminal-age data:
+     * 1) the CHeB lifetime (Myr)
+     * 2) the CHeB TA luminosity (Lsun)
+     *
+     * Unit: Myr
+     *
+     * If *luminosity is non-NULL, set it also
+     */
+    double params[1] = {
+        log10(mass),
+    };
+    double result[2];
+    /*
+     * Interpolate
+     * params column is: log10(mass)
+     * result columns are : log10(age(Myr)), log10(luminosity)
+     */
+    Interpolate(stardata->store->MINT_tables[MINT_TABLE_CHeB_TA],
+                params,
+                result,
+                FALSE);
+
+    if(CHeB_lifetime != NULL)
+    {
+        *CHeB_lifetime = exp10(result[0]);
+    }
+    if(CHeB_TA_luminosity != NULL)
+    {
+        *CHeB_TA_luminosity = exp10(result[1]);
+    }
+
+    return;
+}
+#endif // MINT
diff --git a/src/MINT/MINT_CHeB_ZA.c b/src/MINT/MINT_CHeB_ZA.c
new file mode 100644
index 0000000000000000000000000000000000000000..65bf067ec92e1514bdb439f1bcdab4b8e5d64fbe
--- /dev/null
+++ b/src/MINT/MINT_CHeB_ZA.c
@@ -0,0 +1,40 @@
+#include "../binary_c.h"
+
+No_empty_translation_unit_warning;
+
+#ifdef MINT
+
+void MINT_CHeB_ZA(struct stardata_t * const stardata,
+                  const double mass,
+                  double * const CHeB_ZA_luminosity)
+{
+    /*
+     * Estimate CHeB zero-age data:
+     * 1) the CHeB ZA luminosity (Lsun)
+     *
+     * Unit: Myr
+     *
+     * If *luminosity is non-NULL, set it also
+     */
+    double params[1] = {
+        log10(mass),
+    };
+    double result[1];
+    /*
+     * Interpolate
+     * params column is: log10(mass)
+     * result columns are : log10(luminosity)
+     */
+    Interpolate(stardata->store->MINT_tables[MINT_TABLE_CHeB_ZA],
+                params,
+                result,
+                FALSE);
+
+    if(CHeB_ZA_luminosity != NULL)
+    {
+        *CHeB_ZA_luminosity = exp10(result[1]);
+    }
+
+    return;
+}
+#endif // MINT
diff --git a/src/MINT/MINT_CHeB_get_chebyshev_data.c b/src/MINT/MINT_CHeB_get_chebyshev_data.c
new file mode 100644
index 0000000000000000000000000000000000000000..988e66df3b126a76be8d1c92780652aa5143dfab
--- /dev/null
+++ b/src/MINT/MINT_CHeB_get_chebyshev_data.c
@@ -0,0 +1,25 @@
+#include "../binary_c.h"
+No_empty_translation_unit_warning;
+
+
+#ifdef MINT
+
+/*
+ * Get the data on the Chebyshev coordinates for a CHeB star of
+ * known mass and central hydrogen abundance.
+ */
+
+void MINT_CHeB_get_chebyshev_data(struct stardata_t * RESTRICT const stardata MAYBE_UNUSED,
+                                  struct star_t * const star,
+                                  double * const result_cheb)
+{
+    Interpolate(stardata->store->MINT_tables[MINT_TABLE_CHeB],
+                ((double[]){
+                    log10(star->mass),
+                    star->mint->XHec
+                }), /* log(mass), Xc */
+                result_cheb,
+                FALSE);
+}
+
+#endif // MINT
diff --git a/src/MINT/MINT_CHeB_set_chebyshev_data.c b/src/MINT/MINT_CHeB_set_chebyshev_data.c
new file mode 100644
index 0000000000000000000000000000000000000000..18b99868e7491fb554438b670e7fcedf395774be
--- /dev/null
+++ b/src/MINT/MINT_CHeB_set_chebyshev_data.c
@@ -0,0 +1,88 @@
+#include "../binary_c.h"
+No_empty_translation_unit_warning;
+
+
+#ifdef MINT
+#include "MINT.h"
+
+/*
+ * Set T, rho, etc. which are interpolated
+ * on the Chebyshev mass grid
+ */
+
+int MINT_CHeB_set_chebyshev_data(struct stardata_t * RESTRICT const stardata MAYBE_UNUSED,
+                                 struct star_t * const newstar,
+                                 struct mint_shell_t * const shell,
+                                 const double * const coords,
+                                 const double * const result_cheb,
+                                 const unsigned int nin,
+                                 const unsigned int * const map,
+                                 const Boolean * const available)
+{
+    const double qcoord = shell->m / newstar->mass; /* relative mass coord */
+    const unsigned int ncheb = RCHash_nest_get_entry(
+        stardata->store->MINT_headers[MINT_TABLE_CHeB]->hash,
+        "Chebyshev grid",
+        "number of masses"
+        )->value.value.int_data;
+
+    /*
+     * find the next Chebyshev coordinate just outside this shell
+     */
+    unsigned int n = nin>0 ? nin-1 : 0;
+    while(coords[n] < qcoord &&
+          n < ncheb - 1)
+    {
+        n++;
+    }
+    if(n==0) n=1;
+    if(unlikely(n>=ncheb)) n = ncheb - 1; /* should not happen */
+
+    /*
+     * Interpolation in the Chebyshev grid:
+     * either linear or power law
+     */
+#define _interpolate_linear(F,INDEX)            \
+    ((1.0-(F))*result_cheb[(INDEX)+n-1] +       \
+     (F)*result_cheb[(INDEX)+n])
+#define _interpolate_powerlaw(F,INDEX)                  \
+    (exp10(                                             \
+        ((1.0-(F))*log10(result_cheb[(INDEX)+n-1]) +    \
+         (F)*log10(result_cheb[(INDEX)+n]))             \
+        ))
+#define _interpolate(F,INDEX)                   \
+    (_interpolate_powerlaw((F),(INDEX)))
+
+/*
+ * __interpolate returns the data or 0.0 if not available
+ */
+#define __interpolate(INDEX)                    \
+    (available[(INDEX)] == TRUE                 \
+     ? (_interpolate(f,map[(INDEX)]))           \
+     : 0.0)
+
+    /* interpolate between Chebyshev coordinate n-1 and n */
+    const double f = Limit_range((qcoord - coords[n-1]) / (coords[n] - coords[n-1]),
+                                 0.0,1.0);
+
+    shell->T = __interpolate(MINT_CHeB_CHEBYSHEV_TEMPERATURE);
+
+    shell->T = exp10(log10(shell->T)+0.1);
+    shell->rho = __interpolate(MINT_CHeB_CHEBYSHEV_DENSITY);
+    shell->total_pressure = __interpolate(MINT_CHeB_CHEBYSHEV_TOTAL_PRESSURE);
+    shell->gas_pressure = __interpolate(MINT_CHeB_CHEBYSHEV_GAS_PRESSURE);
+    shell->radius = __interpolate(MINT_CHeB_CHEBYSHEV_RADIUS);
+    shell->radiation_pressure = shell->total_pressure - shell->gas_pressure;
+    shell->gamma1 = __interpolate(MINT_CHeB_CHEBYSHEV_GAMMA1);
+    shell->pressure_scale_height = __interpolate(MINT_CHeB_CHEBYSHEV_PRESSURE_SCALE_HEIGHT);
+
+    /*
+     * Prevent table problems
+     */
+    shell->rho = Max(1e-20,shell->rho);
+    shell->T = Max(300.0,shell->T);
+
+    return n;
+}
+
+#endif // MINT
diff --git a/src/MINT/MINT_Kippenhahn.c b/src/MINT/MINT_Kippenhahn.c
index f15011fb5073897e7a7eec222ef8ccb00897d5c0..2cefd067af59a0e427d752919e0b2f4cc2c1e1ad 100644
--- a/src/MINT/MINT_Kippenhahn.c
+++ b/src/MINT/MINT_Kippenhahn.c
@@ -55,10 +55,11 @@ static void _outkipp(struct stardata_t * const stardata,
 {
     Foreach_shell(shell)
     {
-        Printf("%s %g %g %g %g ",
+        Printf("%s %g %g %g %g %g ",
                header,
                stardata->model.time,
                shell->m,
+               shell->dm,
                shell->T,
                shell->rho);
 #ifdef NUCSYN
@@ -75,7 +76,7 @@ static void _outkipp(struct stardata_t * const stardata,
 
         Printf("\n");
     }
-//    Printf("%s\n",header);
+    Printf("%s\n",header);
 }
 
 #endif // MINT
diff --git a/src/MINT/MINT_data_tables.def b/src/MINT/MINT_data_tables.def
index 24ff408f6a20883803ba07fb96c7dd0c761536e7..bbb3c3ec0abbc1ca07de2fbad128ee642612176b 100644
--- a/src/MINT/MINT_data_tables.def
+++ b/src/MINT/MINT_data_tables.def
@@ -36,6 +36,8 @@
     X(ONeWD,NULL,NULL,FALSE)                                            \
     X(MS_TAMS,NULL,&MINT_make_TAMS_table,TRUE)                          \
     X(MS_ZAMS,NULL,&MINT_make_ZAMS_table,TRUE)                          \
+    X(CHeB_ZA,NULL,&MINT_make_CHeB_ZA_table,TRUE)                       \
+    X(CHeB_TA,NULL,&MINT_make_CHeB_TA_table,TRUE)                       \
     X(ZAMS_COMPOSITION,"ZAMS/zams_composition_ZZ.dat",&MINT_load_ZAMS_composition_grid,FALSE) \
     X(NUMBER,NULL,NULL,FALSE)
 
diff --git a/src/MINT/MINT_init_stars.c b/src/MINT/MINT_init_stars.c
index 79114924cac0c3657d4fb0c61e7cb3edb04a7e55..42121d8ddf05495b0919cd0c332d08b6a07d889b 100644
--- a/src/MINT/MINT_init_stars.c
+++ b/src/MINT/MINT_init_stars.c
@@ -19,7 +19,10 @@ void MINT_init_stars(struct stardata_t * const stardata)
 #endif //NUCSYN
     Foreach_star(star)
     {
-        star->mint->nshells = 200;
+        star->mint->nshells =
+            stardata->preferences->MINT_nshells > 0 ?
+            stardata->preferences->MINT_nshells : 200;
+
         Foreach_shell(shell)
         {
             shell->dm = star->mass / star->mint->nshells;
diff --git a/src/MINT/MINT_load_CHeB_grid.c b/src/MINT/MINT_load_CHeB_grid.c
index 88dcf362ae45e2dba232b6c44efa6233f1aac2e7..6847a284667df8768725e156cc03ae4a9c44d648 100644
--- a/src/MINT/MINT_load_CHeB_grid.c
+++ b/src/MINT/MINT_load_CHeB_grid.c
@@ -16,7 +16,6 @@ Boolean MINT_load_CHeB_grid(struct stardata_t * const stardata)
 {
     const Boolean vb = FALSE; /* set vb=TRUE for lots of verbose output */
 
-
     /*
      * Make lists of parameter and data item names
      * we require in our table of scalars
@@ -29,7 +28,7 @@ Boolean MINT_load_CHeB_grid(struct stardata_t * const stardata)
         { MINT_CHeB_DATA_ITEMS_LIST };
 #undef X
 
-    return MINT_Load_Table(
+    const Boolean ret = MINT_Load_Table(
         stardata,
         MINT_TABLE_CHeB,
         NULL,
@@ -39,6 +38,66 @@ Boolean MINT_load_CHeB_grid(struct stardata_t * const stardata)
         MINT_CHeB_parameter_actions,
         MINT_CHeB_data_actions,
         vb);
+
+
+    /*
+     * Check derivative
+     */
+    struct data_table_t * table = stardata->store->MINT_tables[MINT_TABLE_CHeB];
+    double Ywas = 1.0;
+    double twas = 0.0;
+    double Yinit = -1.0;
+
+    size_t line = table->nlines - 1;
+    while(line != (size_t)-1)
+    {
+        double * const data = table->data + line * (table->nparam + table->ndata);
+        const double Y = data[1];
+        if(Yinit < 0.0) Yinit = Y;
+        const double t = data[4];
+        const double dYdt_table = -exp10(data[21]);
+
+        const double dY = Y - Ywas;
+        const double dt = t - twas;
+
+        if(vb)
+            printf("Line %zu/%zu : t=%g (was %g) dt = %g, Y = %g dY = %g ",
+                   line,
+                   table->nlines,
+                   t,
+                   twas,
+                   dt,
+                   Y,
+                   dY);
+        if(dt > TINY)
+        {
+            const double dYdt =
+                line > 0 ?
+                ((Y - Ywas)/(t-twas)) :
+                0.0;
+            if(vb)
+                printf("dY/dt table = %g mine = %g (tscls %g %g)",
+                       dYdt_table,
+                       dYdt,
+                       -Yinit/dYdt_table,
+                       -Yinit/dYdt
+                );
+            twas = t;
+            Ywas = Y;
+
+            /* replace? */
+            //data[21] = log10(-dYdt);
+        }
+        else if(dt<-TINY)
+        {
+            line = 0;
+        }
+        if(vb) printf("\n");
+        line--;
+    }
+
+//    Flexit;
+    return ret;
 }
 
 #endif // MINT
diff --git a/src/MINT/MINT_macros.h b/src/MINT/MINT_macros.h
index 87d6f1a881f59ba63ba96f631e7477e9dec92c25..24b48459fe265afa294abf7e2486a00c3dd042bb 100644
--- a/src/MINT/MINT_macros.h
+++ b/src/MINT/MINT_macros.h
@@ -42,7 +42,7 @@
 /*
  * Number of shells in each star
  */
-#define MINT_MAX_NSHELLS 1000
+#define MINT_MAX_NSHELLS 10000
 #define MINT_MIN_NSHELLS 10
 
 /*
diff --git a/src/MINT/MINT_make_CHeB_TA_table.c b/src/MINT/MINT_make_CHeB_TA_table.c
new file mode 100644
index 0000000000000000000000000000000000000000..21066c42a673be39370368cd6a63ab01e239b2bb
--- /dev/null
+++ b/src/MINT/MINT_make_CHeB_TA_table.c
@@ -0,0 +1,94 @@
+#include "../binary_c.h"
+No_empty_translation_unit_warning;
+
+#ifdef MINT
+#include "MINT.h"
+
+#define vprint(...) printf(__VA_ARGS__);fflush(stdout);
+
+Boolean MINT_make_CHeB_TA_table(struct stardata_t * const stardata)
+{
+    /*
+     * Terminal-age CHeB (CHeB_TA) data
+     *
+     * Return TRUE on success, FALSE on failure.
+     */
+    const Boolean vb = FALSE;
+    double logM;
+    const int N = 1000;
+    const double min = log10(0.01);
+    const double max = log10(1000.0);
+    const double dlogM = (max - min) / (N - 1);
+
+    /* one parameter : logM */
+    const int nparams = 1;
+    /* two data items: age and luminosity at the end of CHeB */
+    const int ndata = 2;
+
+    double * table_lifetimes = Malloc(sizeof(double) * N * (nparams + ndata));
+
+    if(table_lifetimes == NULL)
+    {
+        Exit_binary_c(BINARY_C_MALLOC_FAILED,
+                      "Failed to malloc for table_liftimes\n");
+    }
+
+    double * const table_lifetimes_start = table_lifetimes;
+    if(vb)
+    {
+        vprint("PRE %p : nparams %zu ndata %zu nlines %zu : table_lifetimes is at %p\n",
+               (void*)stardata->store->MINT_tables[MINT_TABLE_CHeB_TA],
+               stardata->store->MINT_tables[MINT_TABLE_CHeB]->nparam,
+               stardata->store->MINT_tables[MINT_TABLE_CHeB]->ndata,
+               stardata->store->MINT_tables[MINT_TABLE_CHeB]->nlines,
+               (void*)table_lifetimes);
+    }
+
+    double * result = Malloc(MINT_result_size(MINT_TABLE_CHeB));
+
+    for(logM = min; logM < max+TINY; logM += dlogM)
+    {
+        Interpolate(stardata->store->MINT_tables[MINT_TABLE_CHeB],
+                    ((double[]){logM,1.0}), /* log10(mass), Yc=1.0 (or nearest) */
+                    result,
+                    FALSE);
+        const double start_age = result[MINT_CHeB_AGE];
+
+        Interpolate(stardata->store->MINT_tables[MINT_TABLE_CHeB],
+                    ((double[]){logM,0.0}), /* log10(mass), Yc=0.0 (or nearest) */
+                    result,
+                    FALSE);
+
+        const double lifetime = result[MINT_CHeB_AGE] - start_age; // yr
+
+        /*
+         * columns are:
+         * log10(M/Msun),
+         * log10(TAMS age (Myr)), log10(TAMS luminosity)
+         */
+        *(table_lifetimes++) = logM;
+        *(table_lifetimes++) = log10(lifetime) - 6.0; // log10(Myr)
+        *(table_lifetimes++) = log10(result[MINT_CHeB_LUMINOSITY]);
+
+        if(vb)vprint("Lifetimes table M=%g lifetime=%g-%g=%g Myr L(TAMS)=%g %s\n",
+                     exp10(logM),
+                     result[MINT_CHeB_AGE],
+                     start_age,
+                     result[MINT_CHeB_AGE] - start_age,
+                     exp10(result[MINT_MS_LUMINOSITY]),
+                     L_SOLAR
+            );
+
+    }
+
+    NewDataTable_from_Pointer(table_lifetimes_start,
+                              stardata->store->MINT_tables[MINT_TABLE_CHeB_TA],
+                              nparams,
+                              ndata,
+                              N);
+    Safe_free(result);
+
+    return TRUE;
+}
+
+#endif // MINT
diff --git a/src/MINT/MINT_make_CHeB_ZA_table.c b/src/MINT/MINT_make_CHeB_ZA_table.c
new file mode 100644
index 0000000000000000000000000000000000000000..08e7ab9701bd6a3f5f33f6184f6a626b5bafb881
--- /dev/null
+++ b/src/MINT/MINT_make_CHeB_ZA_table.c
@@ -0,0 +1,75 @@
+#include "../binary_c.h"
+No_empty_translation_unit_warning;
+
+#ifdef MINT
+#include "MINT.h"
+
+#define vprint(...) printf(__VA_ARGS__);fflush(stdout);
+
+Boolean MINT_make_CHeB_ZA_table(struct stardata_t * const stardata)
+{
+    /*
+     * Zero-age CHeB (CHeB_ZA) data
+     *
+     * Return TRUE on success, FALSE on failure.
+     */
+    const Boolean vb = FALSE;
+    double logM;
+    const int N = 1000;
+    const double min = log10(0.01);
+    const double max = log10(1000.0);
+    const double dlogM = (max - min) / (N - 1);
+
+    /* one parameter : logM */
+    const int nparams = 1;
+    /* one data item: luminosity at the CHeB ZA */
+    const int ndata = 1;
+
+    double * table_CHeB_ZA = Malloc(sizeof(double) * N * (nparams + ndata));
+
+    if(table_CHeB_ZA == NULL)
+    {
+        Exit_binary_c(BINARY_C_MALLOC_FAILED,
+                      "Failed to malloc for table_liftimes\n");
+    }
+
+    double * const table_CHeB_ZA_start = table_CHeB_ZA;
+    if(vb)
+    {
+        vprint("PRE %p : nparams %zu ndata %zu nlines %zu : table_CHeB_ZA is at %p\n",
+               (void*)stardata->store->MINT_tables[MINT_TABLE_CHeB_ZA],
+               stardata->store->MINT_tables[MINT_TABLE_CHeB]->nparam,
+               stardata->store->MINT_tables[MINT_TABLE_CHeB]->ndata,
+               stardata->store->MINT_tables[MINT_TABLE_CHeB]->nlines,
+               (void*)table_CHeB_ZA);
+    }
+
+    double * result = Malloc(MINT_result_size(MINT_TABLE_CHeB));
+
+    for(logM = min; logM < max+TINY; logM += dlogM)
+    {
+        Interpolate(stardata->store->MINT_tables[MINT_TABLE_CHeB],
+                    ((double[]){logM,1.0}), /* log10(mass), Yc=1.0 (or nearest) */
+                    result,
+                    FALSE);
+
+        /*
+         * columns are:
+         * 1 log10(M/Msun)
+         * 2 log10(CHeB ZA luminosity)
+         */
+        *(table_CHeB_ZA++) = logM;
+        *(table_CHeB_ZA++) = log10(result[MINT_CHeB_LUMINOSITY]);
+    }
+
+    NewDataTable_from_Pointer(table_CHeB_ZA_start,
+                              stardata->store->MINT_tables[MINT_TABLE_CHeB_ZA],
+                              nparams,
+                              ndata,
+                              N);
+    Safe_free(result);
+
+    return TRUE;
+}
+
+#endif // MINT
diff --git a/src/MINT/MINT_merge_shells.c b/src/MINT/MINT_merge_shells.c
new file mode 100644
index 0000000000000000000000000000000000000000..c264728b49351d0eb47b5efb7704d796fcdfefaa
--- /dev/null
+++ b/src/MINT/MINT_merge_shells.c
@@ -0,0 +1,45 @@
+#include "../binary_c.h"
+No_empty_translation_unit_warning;
+
+#ifdef MINT
+
+void MINT_merge_shells(struct stardata_t * const stardata MAYBE_UNUSED,
+                       struct star_t * const star MAYBE_UNUSED,
+                       struct mint_shell_t * const from,
+                       struct mint_shell_t * const into)
+{
+    /*
+     * Merge shell "from" with shell "into" leaving the results in
+     * "into"
+     */
+#define _combine(V) into->V=((from->V * from->dm + into->V * into->dm)/(from->dm + into->dm)); from->V = 0.0;
+
+    /*
+     * merge most variables by averaging by mass
+     */
+    _combine(T);
+    _combine(rho);
+    _combine(gas_pressure);
+    _combine(radiation_pressure);
+    _combine(total_pressure);
+    _combine(radius);
+
+#ifdef NUCSYN
+    for(Isotope i=0; i<ISOTOPE_ARRAY_SIZE;i++)
+    {
+        _combine(X[i]);
+    }
+#endif // NUCSYN
+
+    /* if one is convective, assume both are */
+    into->convective = from->convective || into->convective;
+
+    /* combine mass */
+    into->dm += from->dm;
+    from->dm = -1.0;
+
+    /* reduce the shell count */
+    star->mint->nshells--;
+}
+
+#endif // MINT
diff --git a/src/MINT/MINT_nucsyn_CHeB_pre.c b/src/MINT/MINT_nucsyn_CHeB_pre.c
index 3e995d3c5f060c5934a587fed3a9034115b2f0f0..ebd9244649190cc2c63d3290653fb140917bf7d6 100644
--- a/src/MINT/MINT_nucsyn_CHeB_pre.c
+++ b/src/MINT/MINT_nucsyn_CHeB_pre.c
@@ -8,74 +8,78 @@ void MINT_nucsyn_CHeB_pre(struct stardata_t * RESTRICT const stardata,
                           struct star_t * RESTRICT const star)
 {
     /*
-     * Disable shellular nuclear burning
+     * Enable shellular nuclear burning
      */
-    star->do_burn = FALSE;
+    star->do_burn = TRUE;
 
-    /*
-     * Set the (convective) helium-burning core mass.
-     *
-     * BSE does not give us this, so assume it is half
-     * the core mass (~true for M=2.5 Z=0.02)
-     */
-    const double mc_He = star->core_mass[CORE_He] * 0.5;
+    if(star->do_burn == FALSE)
+    {
+        /*
+         * Set the (convective) helium-burning core mass.
+         *
+         * BSE does not give us this, so assume it is half
+         * the core mass (~true for M=2.5 Z=0.02)
+         */
+        const double mc_He = star->convective_core_mass;
 
-    /* maximum possible helium content */
-    const double max_Y = 1.0 - stardata->common.metallicity;
+        /* maximum possible helium content */
+        const double max_Y = 1.0 - stardata->common.metallicity;
 
-    /*
-     * He-burning lifetime from BSE
-     */
-    double GB[GB_ARRAY_SIZE];
-    double lums[LUMS_ARRAY_SIZE];
-    double tscls[TSCLS_ARRAY_SIZE];
-    stellar_timescales(star->stellar_type,
-                       star->phase_start_mass,
-                       star->mass,
-                       &star->tm,
-                       &star->tn,
-                       tscls,
-                       lums,
-                       GB,
-                       stardata,
-                       star);
-    const double dY_dt = -max_Y / (tscls[T_HE_BURNING]*1e6);
-    const double dY = stardata->model.dt * dY_dt;
+        /*
+         * He-burning lifetime from BSE
+         */
+        double GB[GB_ARRAY_SIZE];
+        double lums[LUMS_ARRAY_SIZE];
+        double tscls[TSCLS_ARRAY_SIZE];
+        stellar_timescales(star->stellar_type,
+                           star->phase_start_mass,
+                           star->mass,
+                           &star->tm,
+                           &star->tn,
+                           tscls,
+                           lums,
+                           GB,
+                           stardata,
+                           star);
+        const double dY_dt = -max_Y / (tscls[T_HE_BURNING]*1e6);
+        const double dY = stardata->model.dt * dY_dt;
 
-    Foreach_shell(shell)
-    {
-        if(shell->m < mc_He)
+
+        Foreach_shell(shell)
         {
-            /*
-             * Evolve helium in the core.
-             *
-             * We assume all the core is burning helium,
-             * such that at the beginning of CHeB Y=1-Z,
-             * and at the end Y=0.
-             *
-             * Shells which are newly ingested in the core
-             * have X>1, so in those we immediately
-             * convert H to He, and CNO to N, prior to
-             * helium burning.
-             */
-            if(shell->X[XH1] > 0.0)
+            if(shell->m < mc_He)
             {
-                shell->X[XHe4] += shell->X[XH1];
-                shell->X[XH1] = 0.0;
-                shell->X[XN14] += shell->X[XC12] + shell->X[XC13] + shell->X[XN15] + shell->X[XO16] + shell->X[XO17] + shell->X[XO18];
-                shell->X[XC12] = 0.0;
-                shell->X[XC13] = 0.0;
-                shell->X[XN15] = 0.0;
-                shell->X[XO16] = 0.0;
-                shell->X[XO17] = 0.0;
-                shell->X[XO18] = 0.0;
+                /*
+                 * Evolve helium in the core.
+                 *
+                 * We assume all the core is burning helium,
+                 * such that at the beginning of CHeB Y=1-Z,
+                 * and at the end Y=0.
+                 *
+                 * Shells which are newly ingested in the core
+                 * have X>1, so in those we immediately
+                 * convert H to He, and CNO to N, prior to
+                 * helium burning.
+                 */
+                if(shell->X[XH1] > 0.0)
+                {
+                    shell->X[XHe4] += shell->X[XH1];
+                    shell->X[XH1] = 0.0;
+                    shell->X[XN14] += shell->X[XC12] + shell->X[XC13] + shell->X[XN15] + shell->X[XO16] + shell->X[XO17] + shell->X[XO18];
+                    shell->X[XC12] = 0.0;
+                    shell->X[XC13] = 0.0;
+                    shell->X[XN15] = 0.0;
+                    shell->X[XO16] = 0.0;
+                    shell->X[XO17] = 0.0;
+                    shell->X[XO18] = 0.0;
+                }
+                shell->X[XHe4] += dY;
+                shell->X[XHe4] = Max(0.0,shell->X[XHe4]);
+                shell->X[XC12] -= 0.5 * dY;
+                shell->X[XO16] -= 0.5 * dY;
+                shell->X[XNe22] += shell->X[XN14];
+                shell->X[XN14] = 0.0;
             }
-            shell->X[XHe4] += dY;
-            shell->X[XHe4] = Max(0.0,shell->X[XHe4]);
-            shell->X[XC12] -= 0.5 * dY;
-            shell->X[XO16] -= 0.5 * dY;
-            shell->X[XNe22] += shell->X[XN14];
-            shell->X[XN14] = 0.0;
         }
     }
 }
diff --git a/src/MINT/MINT_nucsyn_functions.def b/src/MINT/MINT_nucsyn_functions.def
index a912358a131c03e840277c08bd7e2beab9a03e04..fba70c65aa3b166be031a6d14d9ccc497e763c79 100644
--- a/src/MINT/MINT_nucsyn_functions.def
+++ b/src/MINT/MINT_nucsyn_functions.def
@@ -19,17 +19,17 @@
         X(MAIN_SEQUENCE, &MINT_nucsyn_MS_pre, TRUE, &MINT_nucsyn_MS_post) \
         X(HERZTSPRUNG_GAP, &MINT_nucsyn_HG_pre, FALSE, NULL)            \
         X(GIANT_BRANCH, &MINT_nucsyn_GB_pre, FALSE, NULL)               \
-        X(CORE_HELIUM_BURNING, &MINT_nucsyn_CHeB_pre, FALSE, NULL)      \
+        X(CORE_HELIUM_BURNING, &MINT_nucsyn_CHeB_pre, TRUE, NULL)       \
         X(EAGB, &MINT_nucsyn_EAGB_pre, FALSE, NULL)                     \
-        X(TPAGB, &MINT_nucsyn_TPAGB_pre, FALSE, NULL)                   \
-        X(HeMS, &MINT_nucsyn_HeMS_pre, FALSE, NULL)                     \
-        X(HeHG, &MINT_nucsyn_HeHG_pre, FALSE, NULL)                     \
-        X(HeGB, &MINT_nucsyn_HeGB_pre, FALSE, NULL)                     \
-        X(HeWD, &MINT_nucsyn_HeWD_pre, FALSE, NULL)                     \
-        X(COWD, &MINT_nucsyn_COWD_pre, FALSE, NULL)                     \
-        X(ONeWD, &MINT_nucsyn_ONeWD_pre, FALSE, NULL)                   \
-        X(NS, &MINT_nucsyn_NS_pre, FALSE, NULL)                         \
-        X(BH, &MINT_nucsyn_BH_pre, FALSE, NULL)
+    X(TPAGB, &MINT_nucsyn_TPAGB_pre, FALSE, NULL)                       \
+    X(HeMS, &MINT_nucsyn_HeMS_pre, FALSE, NULL)                         \
+    X(HeHG, &MINT_nucsyn_HeHG_pre, FALSE, NULL)                         \
+    X(HeGB, &MINT_nucsyn_HeGB_pre, FALSE, NULL)                         \
+    X(HeWD, &MINT_nucsyn_HeWD_pre, FALSE, NULL)                         \
+    X(COWD, &MINT_nucsyn_COWD_pre, FALSE, NULL)                         \
+    X(ONeWD, &MINT_nucsyn_ONeWD_pre, FALSE, NULL)                       \
+    X(NS, &MINT_nucsyn_NS_pre, FALSE, NULL)                             \
+    X(BH, &MINT_nucsyn_BH_pre, FALSE, NULL)
 
 
 #endif // MINT_NUCSYN_FUNCTIONS_DEF
diff --git a/src/MINT/MINT_prototypes.h b/src/MINT/MINT_prototypes.h
index 1ec0d4c3ac2f3fc0a3bca4d2efd069b5b2629b08..1e0d6b24ceff910db10f305db533cccb88319fc4 100644
--- a/src/MINT/MINT_prototypes.h
+++ b/src/MINT/MINT_prototypes.h
@@ -52,7 +52,8 @@ Boolean MINT_load_ZAMS_composition_grid(struct stardata_t * const stardata);
 Boolean MINT_make_ZAMS_table(struct stardata_t * const stardata);
 Boolean MINT_make_TAMS_table(struct stardata_t * const stardata);
 Boolean MINT_load_CHeB_grid(struct stardata_t * const stardata);
-
+Boolean MINT_make_CHeB_ZA_table(struct stardata_t * const stardata);
+Boolean MINT_make_CHeB_TA_table(struct stardata_t * const stardata);
 
 /*
  * Stellar structure
@@ -112,16 +113,45 @@ Boolean MINT_rejuvenate_MS(struct stardata_t * const stardata,
                            struct star_t * star_was,
                            struct star_t * star_is);
 
+/*
+ * Stellar-evolution phase functions
+ */
 double MINT_MS_lifetime(struct stardata_t * const stardata,
                         struct star_t * const star);
+void MINT_CHeB_ZA(struct stardata_t * const stardata,
+                  const double mass,
+                  double * const CHeB_ZA_luminosity);
+void MINT_CHeB_TA(struct stardata_t * const stardata,
+                  const double mass,
+                  double * const CHeB_lifetime,
+                  double * const CHeB_TA_luminosity);
+/*
+ * Initial conditions
+ */
 double MINT_initial_XHc(struct stardata_t * const stardata);
 
 void MINT_make_massless_remnant(struct stardata_t * const stardata,
                                 struct star_t * const star);
 
+/*
+ * Meshing
+ */
 void MINT_remesh(struct stardata_t * const stardata,
                  struct star_t * const star);
-
+void MINT_split_shell(struct stardata_t * const stardata MAYBE_UNUSED,
+                      struct star_t * const star MAYBE_UNUSED,
+                      const int k);
+void MINT_merge_shells(struct stardata_t * const stardata MAYBE_UNUSED,
+                       struct star_t * const star MAYBE_UNUSED,
+                       struct mint_shell_t * const from,
+                       struct mint_shell_t * const into);
+void MINT_resolve_core_boundary(struct stardata_t * RESTRICT const stardata MAYBE_UNUSED,
+                                struct star_t * const star,
+                                const double band_centre_mass,
+                                const double band_thickness_mass,
+                                const double thinnest_dm,
+                                const double default_dm,
+                                const Boolean bail_above_core);
 /*
  * Derivative application function
  */
@@ -149,6 +179,7 @@ void MINT_shellular_burning(struct stardata_t * RESTRICT const stardata,
                             double * const convective_sigmav,
                             double * const Xconv,
                             double * const mconv);
+
 struct mint_shell_t * MINT_mixed_shell(struct stardata_t * const stardata,
                                        struct star_t * const star,
                                        const double m_low,
@@ -218,10 +249,12 @@ void MINT_convection(struct stardata_t * const stardata MAYBE_UNUSED,
 double MINT_MS_non_eq_conv_coremass(struct stardata_t * const stardata,
                                     struct star_t * star);
 
+/*
+ * Chebyshev grid functions
+ */
 void MINT_MS_get_chebyshev_data(struct stardata_t * RESTRICT const stardata MAYBE_UNUSED,
                                 struct star_t * const newstar,
                                 double * const result_cheb);
-
 int MINT_MS_set_chebyshev_data(struct stardata_t * RESTRICT const stardata MAYBE_UNUSED,
                                struct star_t * const newstar,
                                struct mint_shell_t * const shell,
@@ -230,6 +263,18 @@ int MINT_MS_set_chebyshev_data(struct stardata_t * RESTRICT const stardata MAYBE
                                const unsigned int nin,
                                const unsigned int * const map,
                                const Boolean * const available);
+void MINT_CHeB_get_chebyshev_data(struct stardata_t * RESTRICT const stardata MAYBE_UNUSED,
+                                  struct star_t * const star,
+                                  double * const result_cheb);
+int MINT_CHeB_set_chebyshev_data(struct stardata_t * RESTRICT const stardata MAYBE_UNUSED,
+                                 struct star_t * const newstar,
+                                 struct mint_shell_t * const shell,
+                                 const double * const coords,
+                                 const double * const result_cheb,
+                                 const unsigned int nin,
+                                 const unsigned int * const map,
+                                 const Boolean * const available);
+
 
 double MINT_ZAMS_core_mass_from_luminosity(struct stardata_t * const stardata,
                                            struct star_t * star,
diff --git a/src/MINT/MINT_remesh.c b/src/MINT/MINT_remesh.c
index f16fb12328a9e683e40b891808c4df4df5fc72cf..951751605b7fd14c91eafcac0e69bd02b3af5d94 100644
--- a/src/MINT/MINT_remesh.c
+++ b/src/MINT/MINT_remesh.c
@@ -2,10 +2,6 @@
 No_empty_translation_unit_warning;
 
 #ifdef MINT
-static void merge_shells(struct stardata_t * const stardata MAYBE_UNUSED,
-                         struct star_t * const star MAYBE_UNUSED,
-                         struct mint_shell_t * const from,
-                         struct mint_shell_t * const into);
 
 /*
  * Remesh MINT data
@@ -75,10 +71,10 @@ void MINT_remesh(struct stardata_t * const stardata MAYBE_UNUSED,
                             inner :
                             outer;
 
-                        merge_shells(stardata,
-                                     star,
-                                     with,
-                                     shell);
+                        MINT_merge_shells(stardata,
+                                          star,
+                                          with,
+                                          shell);
 
                         /*
                          * We merge with the inner or outer shell
@@ -111,7 +107,6 @@ void MINT_remesh(struct stardata_t * const stardata MAYBE_UNUSED,
                                     shell,
                                     sizeof(struct mint_shell_t) * (star->mint->nshells - k));
                         }
-                        star->mint->nshells--;
                         if(vb)printf("now have %d shells\n",star->mint->nshells);
                     }
                     surface = &star->mint->shells[star->mint->nshells-1];
@@ -129,40 +124,5 @@ void MINT_remesh(struct stardata_t * const stardata MAYBE_UNUSED,
     if(vb)printf("done remesh\n");
 }
 
-static void merge_shells(struct stardata_t * const stardata MAYBE_UNUSED,
-                         struct star_t * const star MAYBE_UNUSED,
-                         struct mint_shell_t * const from,
-                         struct mint_shell_t * const into)
-{
-    /*
-     * Merge shell "from" with shell "into" leaving the results in
-     * "into"
-     */
-#define _combine(V) into->V=((from->V * from->dm + into->V * into->dm)/(from->dm + into->dm)); from->V = 0.0;
-
-    /*
-     * merge most variables by averaging by mass
-     */
-    _combine(T);
-    _combine(rho);
-    _combine(gas_pressure);
-    _combine(radiation_pressure);
-    _combine(total_pressure);
-    _combine(radius);
-
-#ifdef NUCSYN
-    for(Isotope i=0; i<ISOTOPE_ARRAY_SIZE;i++)
-    {
-        _combine(X[i]);
-    }
-#endif // NUCSYN
-
-    /* if one is convective, assume both are */
-    into->convective = from->convective || into->convective;
-
-    /* combine mass */
-    into->dm += from->dm;
-    from->dm = -1.0;
-}
 
 #endif // MINT
diff --git a/src/MINT/MINT_resolve_core_boundary.c b/src/MINT/MINT_resolve_core_boundary.c
new file mode 100644
index 0000000000000000000000000000000000000000..6bfab71c5992874f62570c903031400f045ba094
--- /dev/null
+++ b/src/MINT/MINT_resolve_core_boundary.c
@@ -0,0 +1,96 @@
+#include "../binary_c.h"
+No_empty_translation_unit_warning;
+
+
+#ifdef MINT
+
+/*
+ * Subroutine to take a band (in mass) in the star
+ * and resolve it well, while keeping to a default
+ * shell thickness outside it.
+ *
+ * band_centre_mass is the band's central mass co-ordinate,
+ *     e.g. the convective core boundary
+ *
+ * band_thickness_mass is the thickness of the band
+ *     that requires the extra resolution
+ *
+ * thinnest_dm is how thin the shells in the band should be
+ *
+ * default_dm is the thickness of the shells that do
+ *     not have extra resolution.
+ *
+ * bail_above_core means once we've looped above the band,
+ * we return immediately. Usually you want this to be TRUE
+ * unless you (for some reason) want to keep scanning the rest
+ * of the star.
+ */
+
+void MINT_resolve_core_boundary(struct stardata_t * RESTRICT const stardata MAYBE_UNUSED,
+                                struct star_t * const star,
+                                const double band_centre_mass,
+                                const double band_thickness_mass,
+                                const double thinnest_dm,
+                                const double default_dm,
+                                const Boolean bail_above_core)
+{
+    int k = 0;
+    const double half_thickness = 0.5 * band_thickness_mass;
+    const double band_low = band_centre_mass - half_thickness;
+    const double band_high = band_centre_mass + half_thickness;
+
+    while(k < star->mint->nshells)
+    {
+        struct mint_shell_t * const shell =
+            &star->mint->shells[k];
+
+        if(k>0 &&
+           shell->m < band_low &&
+           (shell-1)->dm + shell->dm <= default_dm)
+        {
+            /*
+             * Inside the core band, if the shell
+             * is thin, merge shell into shell-1
+             */
+            MINT_merge_shells(stardata,
+                              star,
+                              shell,
+                              shell-1);
+
+            /*
+             * Move all the shells at shell (and above)
+             * down one
+             */
+            memmove(shell,
+                    shell+1,
+                    sizeof(struct mint_shell_t) * (star->mint->nshells - k - 1));
+        }
+        else if(shell->dm > thinnest_dm &&
+                fabs(shell->m - band_centre_mass) < half_thickness)
+        {
+            /*
+             * Split the shell in the band that's not thin enough
+             */
+            MINT_split_shell(stardata,
+                             star,
+                             k);
+        }
+        else
+        {
+            /*
+             * Do nothing: move to next shell
+             */
+            if(bail_above_core == TRUE &&
+               shell->m > band_high)
+            {
+                /*
+                 * Bail out now we're above the band
+                 */
+                return;
+            }
+            k++;
+        }
+    }
+}
+
+#endif//MINT
diff --git a/src/MINT/MINT_shellular_burning.c b/src/MINT/MINT_shellular_burning.c
index 3f22969c3e239aa24502a3bffecd6c8e6528f9b5..07dc1525e09a080de272ce24b46a5e34b7923c82 100644
--- a/src/MINT/MINT_shellular_burning.c
+++ b/src/MINT/MINT_shellular_burning.c
@@ -17,10 +17,33 @@ void MINT_shellular_burning(struct stardata_t * RESTRICT const stardata,
                             double * const mconv)
 {
     /*
-     * Nuclear burning in shells for MINT
+     * If slowmix is TRUE, actually burn all the layers
+     * one by one (which is slow!) and then mix them.
+     * Normally you don't want to do this.
+     */
+    const Boolean slowmix = FALSE;
+
+    /*
+     * Detect central and surface shell
+     */
+    //const Boolean is_central_shell = k == 0 ;
+    const Boolean is_surface_shell = k == star->mint->nshells-1;
+
+    /*
+     * Locate a convection zone outer edge.
      */
+    const Boolean convective_outer_edge =
+        (prevshell->convective == TRUE &&
+         (shell->convective == FALSE ||
+          is_surface_shell));
+
+    //struct mint_shell_t * const nextshell =
+    //    is_surface_shell ? shell : (shell+1);
 
-    Number_density Nshell[ISOTOPE_ARRAY_SIZE];
+    /*
+     * Nuclear burning in shells for MINT
+     */
+    Number_density * Nshell = New_isotope_array;
     Reaction_rate sigmav[SIGMAV_SIZE];
 
     /*
@@ -39,7 +62,7 @@ void MINT_shellular_burning(struct stardata_t * RESTRICT const stardata,
      */
     nucsyn_set_sigmav(stardata,
                       stardata->store,
-                      log10(shell->T),
+                      log10(shell->T)-0.1,
                       sigmav,
                       stardata->preferences->reaction_rate_multipliers
 #ifdef RATES_OF_AMANDA_NE22PG_FIX
@@ -47,122 +70,128 @@ void MINT_shellular_burning(struct stardata_t * RESTRICT const stardata,
 #endif
         );
 
-    if(shell->convective == TRUE)
-    {
-        /*
-         * Calculate contribution to
-         * convective sigmav
-         */
-        const double fshell = shell->dm * shell->rho;
-        for(Reaction_rate_index j=0; j<SIGMAV_SIZE;j++)
-        {
-            convective_sigmav[j] += fshell * sigmav[j];
-        }
-    }
-
-    if(shell->convective == TRUE)
-    {
-        /*
-         * Calculate contribution to the mass of each particle
-         * in the convective zone
-         */
-        for(Isotope j=0; j<ISOTOPE_ARRAY_SIZE; j++)
-        {
-            Xconv[j] += shell->dm * shell->X[j];
-        }
-        *mconv += shell->dm;
-    }
-
-    if(prevshell->convective == TRUE &&
-       (shell->convective == FALSE ||
-        k == star->mint->nshells-1))
+    if(slowmix == FALSE)
     {
-        /*
-         * Reached the outer edge of the convection zone,
-         * or were convective and have reached the surface.
-         *
-         * Finish computation of its mean <sigmav> and
-         * do nuclear burning in the convection zone.
-         */
-        const double mconv_squared = Pow2(*mconv);
-        for(Reaction_rate_index j=0;j<SIGMAV_SIZE;j++)
-        {
-            convective_sigmav[j] /= mconv_squared;
-        }
-
-        /*
-         * Finalise the Xconv calculation
-         */
-        for(Isotope j=0; j<ISOTOPE_ARRAY_SIZE; j++)
+        if(shell->convective == TRUE)
         {
-            Xconv[j] /= *mconv;
+            /*
+             * Calculate contribution to
+             * convective sigmav
+             */
+            const double fshell = shell->dm * shell->rho;
+            for(Reaction_rate_index j=0; j<SIGMAV_SIZE; j++)
+            {
+                convective_sigmav[j] +=
+                    (j == SIGMAV_3He4 ? (shell->rho) : 1.0) *
+                    fshell * sigmav[j];
+            }
+
+            /*
+             * Calculate contribution to the mass of each particle
+             * in the convective zone
+             */
+            for(Isotope j=0; j<ISOTOPE_ARRAY_SIZE; j++)
+            {
+                Xconv[j] += shell->dm * shell->X[j];
+            }
+            *mconv += shell->dm;
         }
 
-        /*
-         * Convert mass fraction to total number of particles
-         */
-        Abundance Nconv[ISOTOPE_ARRAY_SIZE];
-        X_to_N(stardata->store->imnuc,
-               *mconv,
-               Nconv,
-               Xconv,
-               ISOTOPE_ARRAY_SIZE);
-
-        Nconv[Xe] = nucsyn_free_electron_density(Nconv,
-                                                 stardata->store->atomic_number);
-
-        /*
-         * Use the mean convective_sigmav to do the nuclear burning
-         * in the convection zone.
-         */
-        nucsyn_burning_cycles(Nconv,
-                              convective_sigmav,
-                              stardata->model.dt, /* dt in years */
-                              star->mint->shells[0].T,
-                              star->mint->shells[0].rho,
-                              stardata);
-        Nconv[Xe] = 0.0;
-
-        nancheck("MINT burn (convective N) ",Nconv,ISOTOPE_ARRAY_SIZE);
-
-        /*
-         * Convert number of particles back to mass
-         * fractions and set in each shell in the
-         * convective zone
-         */
-        N_to_X(stardata->store->mnuc,
-               *mconv,
-               Nconv,
-               Xconv,
-               ISOTOPE_ARRAY_SIZE);
-        nancheck("MINT burn (convective X) ",Xconv,ISOTOPE_ARRAY_SIZE);
-
 
-        for(Shell_index nn=0; nn<k; nn++)
+        if(convective_outer_edge == TRUE)
         {
-            memcpy(star->mint->shells[nn].X,
+            /*
+             * Reached the outer edge of the convection zone,
+             * or were convective and have reached the surface.
+             *
+             * Finish computation of its mean <sigmav> and
+             * do nuclear burning in the convection zone.
+             */
+            const double mconv_squared = Pow2(*mconv);
+            for(Reaction_rate_index j=0;j<SIGMAV_SIZE;j++)
+            {
+                convective_sigmav[j] /=
+                    (j == SIGMAV_3He4 ? *mconv : 1.0) *
+                    mconv_squared;
+            }
+
+            /*
+             * Finalise the Xconv calculation
+             */
+            for(Isotope j=0; j<ISOTOPE_ARRAY_SIZE; j++)
+            {
+                Xconv[j] /= *mconv;
+            }
+
+            /*
+             * Convert mass fraction to total number of particles
+             */
+            Abundance * Nconv = New_isotope_array;
+            X_to_N(stardata->store->imnuc,
+                   *mconv,
+                   Nconv,
                    Xconv,
-                   ISOTOPE_MEMSIZE);
-        }
-
-        /*
-         * Clear the convective variables
-         * ready for the next convective zone
-         * (if there is one)
-         */
-        *mconv = 0.0;
-        for(Reaction_rate_index j=0;j<SIGMAV_SIZE;j++)
-        {
-            convective_sigmav[j] = 0.0;
-        }
-
-        for(Isotope j=0; j<ISOTOPE_ARRAY_SIZE; j++)
-        {
-            Xconv[j] = 0.0;
+                   ISOTOPE_ARRAY_SIZE);
+
+            Nconv[Xe] = nucsyn_free_electron_density(
+                Nconv,
+                stardata->store->atomic_number);
+
+            /*
+             * Use the mean convective_sigmav to do the nuclear burning
+             * in the convection zone.
+             */
+            nucsyn_burning_cycles(Nconv,
+                                  convective_sigmav,
+                                  stardata->model.dt, /* dt in years */
+                                  star->mint->shells[0].T,
+                                  star->mint->shells[0].rho, // irrelevant
+                                  stardata);
+            Nconv[Xe] = 0.0;
+
+            nancheck("MINT burn (convective N) ",Nconv,ISOTOPE_ARRAY_SIZE);
+
+            /*
+             * Convert number of particles back to mass
+             * fractions and set in each shell in the
+             * convective zone
+             */
+            N_to_X(stardata->store->mnuc,
+                   *mconv,
+                   Nconv,
+                   Xconv,
+                   ISOTOPE_ARRAY_SIZE);
+            nancheck("MINT burn (convective X) ",Xconv,ISOTOPE_ARRAY_SIZE);
+
+            Safe_free(Nconv);
+
+            for(Shell_index nn=0; nn<k; nn++)
+            {
+                memcpy(star->mint->shells[nn].X,
+                       Xconv,
+                       ISOTOPE_MEMSIZE);
+            }
+
+            /*
+             * Clear the convective variables
+             * ready for the next convective zone
+             * (if there is one)
+             */
+            *mconv = 0.0;
+            for(Reaction_rate_index j=0;j<SIGMAV_SIZE;j++)
+            {
+                convective_sigmav[j] = 0.0;
+            }
+
+            for(Isotope j=0; j<ISOTOPE_ARRAY_SIZE; j++)
+            {
+                Xconv[j] = 0.0;
+            }
         }
     }
 
-    if(shell->convective == FALSE)
+    if(slowmix == TRUE ||
+       shell->convective == FALSE)
     {
         /*
          * Radiative shell, normal burn
@@ -181,10 +210,41 @@ void MINT_shellular_burning(struct stardata_t * RESTRICT const stardata,
                Nshell,
                shell->X,
                ISOTOPE_ARRAY_SIZE);
-
         nancheck("MINT burn (radiative X) ",shell->X,ISOTOPE_ARRAY_SIZE);
+    }
 
+    if(slowmix == TRUE &&
+       convective_outer_edge == TRUE)
+    {
+        /*
+         * homogenize convective zone
+         */
+        double m = 0.0;
+        Abundance * X = New_clear_isotope_array;
+        for(Shell_index nn=0; nn<k; nn++)
+        {
+            struct mint_shell_t * const s = &star->mint->shells[nn];
+            for(Isotope j=0; j<ISOTOPE_ARRAY_SIZE; j++)
+            {
+                X[j] += s->X[j] * s->dm;
+            }
+            m += s->dm;
+        }
+        for(Isotope j=0; j<ISOTOPE_ARRAY_SIZE; j++)
+        {
+            X[j] /= m;
+        }
+        for(Shell_index nn=0; nn<k; nn++)
+        {
+            memcpy(star->mint->shells[nn].X,
+                   X,
+                   ISOTOPE_MEMSIZE);
+        }
+        Safe_free(X);
+        fprintf(stderr,"slowmix\n");
     }
+
+    Safe_free(Nshell);
 }
 
 #endif // MINT && NUCSYN
diff --git a/src/MINT/MINT_split_shell.c b/src/MINT/MINT_split_shell.c
new file mode 100644
index 0000000000000000000000000000000000000000..554d81bae5bbbdd3a622755a303fa20e04c0cc00
--- /dev/null
+++ b/src/MINT/MINT_split_shell.c
@@ -0,0 +1,28 @@
+#include "../binary_c.h"
+No_empty_translation_unit_warning;
+
+#ifdef MINT
+
+void MINT_split_shell(struct stardata_t * const stardata MAYBE_UNUSED,
+                      struct star_t * const star MAYBE_UNUSED,
+                      const int k)
+{
+    /*
+     * Split shell k into two identical shells
+     * and increase the number of shells
+     */
+    if(star->mint->nshells < MINT_MAX_NSHELLS - 2)
+    {
+        struct mint_shell_t * const this = &star->mint->shells[k];
+        memmove(this+1,
+                this,
+                sizeof(struct mint_shell_t) * (star->mint->nshells - k));
+        this->dm *= 0.5;
+        (this+1)->m += this->dm;
+        (this+1)->dm *= 0.5;
+        star->mint->nshells++;
+    }
+}
+
+
+#endif // MINT
diff --git a/src/MINT/MINT_stellar_structure.c b/src/MINT/MINT_stellar_structure.c
index 15c4aac1930e5e1da7ecfb48795c9eb2698004ec..9c7f40fa6c00f57fa85375693289e06fd1aab936 100644
--- a/src/MINT/MINT_stellar_structure.c
+++ b/src/MINT/MINT_stellar_structure.c
@@ -23,7 +23,6 @@ int MINT_stellar_structure(struct stardata_t * const stardata,
                            const Caller_id caller_id MAYBE_UNUSED)
 {
     const Boolean vb = FALSE;
-
     /*
      * Trim shells, e.g. because of mass loss
      */
@@ -83,6 +82,11 @@ int MINT_stellar_structure(struct stardata_t * const stardata,
     /*
      * Call BSE before MINT?
      */
+    if(vb)printf("MINT t=%g st = %d action = %u\n",
+                 stardata->model.time,
+                 newstar->stellar_type,
+                 call_bse_action[newstar->stellar_type]);
+
     if(call_bse_action[newstar->stellar_type] == MINT_CALL_BSE_BEFORE ||
        call_bse_action[newstar->stellar_type] == MINT_CALL_BSE_BOTH)
     {
diff --git a/src/MINT/MINT_stellar_structure_CHeB.c b/src/MINT/MINT_stellar_structure_CHeB.c
index 1288fcc00eab30458e62d02bc3bf752ae1b0f69d..81ad9755594dbf1efdcc8138d86d98afdf561db4 100644
--- a/src/MINT/MINT_stellar_structure_CHeB.c
+++ b/src/MINT/MINT_stellar_structure_CHeB.c
@@ -7,16 +7,15 @@ No_empty_translation_unit_warning;
 /*
  * MINT function to calculate stellar structure
  * during core helium burning.
- *
- * Note that prevstar can be NULL, in which case
- * we cannot compare to the previous star to do (say)
- * rejuvenation.
  */
 
 Stellar_type MINT_stellar_structure_CHeB(struct stardata_t * RESTRICT const stardata MAYBE_UNUSED,
-                                       struct star_t * const prevstar MAYBE_UNUSED,
-                                       struct star_t * const newstar)
+                                         struct star_t * const prevstar MAYBE_UNUSED,
+                                         struct star_t * const newstar)
 {
+//#define OLDCODE
+#ifdef OLDCODE
+
     /*
      * Make the envelope cool
      */
@@ -45,6 +44,355 @@ Stellar_type MINT_stellar_structure_CHeB(struct stardata_t * RESTRICT const star
             (shell->m < mc_He) ? TRUE :
             FALSE;
     }
+#else
+
+    const Boolean vb = FALSE;
+    /******************************/
+    /** Core-helium burning star **/
+    /******************************/
+    Dprint("CHeB star : MINT algorithm\n");
+
+    newstar->derivative[DERIVATIVE_STELLAR_He_CORE_MASS] = 0.0;
+    newstar->derivative[DERIVATIVE_STELLAR_GB_CORE_MASS] = 0.0;
+    newstar->derivative[DERIVATIVE_STELLAR_CO_CORE_MASS] = 0.0;
+
+    /*
+     * Set stellar type
+     */
+    newstar->stellar_type = CHeB;
+
+    /*
+     * Booleans to test for data availability
+     */
+    const Boolean * const available = stardata->store->MINT_data_available[MINT_TABLE_CHeB];
+
+    if(available == NULL)
+    {
+        Exit_binary_c(BINARY_C_ALLOC_FAILED,
+                      "stardata->store->MINT_data_available[MINT_TABLE_CHeB] is NULL. This should not happen.");
+    }
+
+    /*
+     * Map between MINT's wanted columns and columns available
+     * in the data table
+     */
+    const unsigned int * const map = stardata->store->MINT_data_map[MINT_TABLE_CHeB];
+
+    if(map == NULL)
+    {
+        Exit_binary_c(BINARY_C_ALLOC_FAILED,
+                      "stardata->store->MINT_data_map[MINT_TABLE_CHeB] is NULL. This should not happen.");
+    }
+
+    /*
+     * Space for interpolations
+     */
+    double * result, * result_cheb, * result_conv;
+    if(stardata->tmpstore->MINT_CHeB_result == NULL)
+    {
+        result =
+            stardata->tmpstore->MINT_CHeB_result =
+            Malloc(MINT_result_size(MINT_TABLE_CHeB));
+        result_cheb =
+            stardata->tmpstore->MINT_CHeB_result_cheb =
+            Malloc(MINT_result_size(MINT_TABLE_CHeB));
+        result_conv =
+            stardata->tmpstore->MINT_CHeB_result_conv =
+            Malloc(MINT_result_size(MINT_TABLE_CHeB));
+
+        if(result == NULL ||
+           result_cheb == NULL ||
+           result_conv == NULL)
+        {
+            Exit_binary_c(BINARY_C_ALLOC_FAILED,
+                          "Failed to allocate memory for MINT MS interpolation tables");
+        }
+    }
+    else
+    {
+        result =
+            stardata->tmpstore->MINT_CHeB_result;
+        result_cheb =
+            stardata->tmpstore->MINT_CHeB_result_cheb;
+        result_conv =
+            stardata->tmpstore->MINT_CHeB_result_conv;
+    }
+
+    /*
+     * Macro to get data from the interpolation, if it is
+     * available, otherwise return 0.0. You can check
+     * available[VARTYPE] yourself if you want to apply
+     * a more complicated algorithm.
+     */
+#undef get_data
+#define get_data(VARTYPE)                                           \
+    ( available[(VARTYPE)]==TRUE ? result[map[(VARTYPE)]] : 0.0 )
+
+    /*
+     * Use table to look up stellar structure, e.g. dY/dt,L,R,
+     * as a function of phase start mass and central hydrogen abundance
+     */
+    Loop_forever
+    {
+        /*
+         * Interpolate to find stellar structure
+         */
+        Interpolate(stardata->store->MINT_tables[MINT_TABLE_CHeB],
+                    ((double[]){
+                        log10(newstar->mass),
+                        newstar->mint->XHec
+                    }),
+                    result,
+                    FALSE);
+
+#if defined DEBUG && DEBUG == 1
+        {
+            char * s = NULL;
+            for(size_t i=0;
+                i<MINT_result_size(MINT_TABLE_CHeB)/sizeof(double);
+                i++)
+            {
+                reasprintf(&s,
+                           "%s %zu=%g",
+                           s == NULL ? "" : s,
+                           i,
+                           result[i]);
+            }
+            Dprint("Interpolate:%s\n",
+                   s);
+            Safe_free(s);
+        }
+#endif
+
+        /*
+         * Hence set the stellar properties
+         */
+        newstar->luminosity = exp10(get_data(MINT_CHeB_LUMINOSITY));
+        newstar->radius = exp10(get_data(MINT_CHeB_RADIUS));
+
+        if(Is_zero(newstar->age))
+        {
+            newstar->rzams = newstar->radius;
+        }
+
+        /*
+         * Zero burning derivatives
+         */
+        newstar->derivative[DERIVATIVE_STELLAR_CENTRAL_HYDROGEN] = 0.0;
+        newstar->derivative[DERIVATIVE_STELLAR_CENTRAL_HELIUM] = 0.0;
+        newstar->derivative[DERIVATIVE_STELLAR_CENTRAL_CARBON] = 0.0;
+
+        /*
+         * burn He to C
+         */
+        const double dycdt = -exp10(get_data(MINT_CHeB_FIRST_DERIVATIVE_CENTRAL_HELIUM));
+
+        newstar->derivative[DERIVATIVE_STELLAR_CENTRAL_HELIUM] = dycdt;
+        newstar->derivative[DERIVATIVE_STELLAR_CENTRAL_CARBON] = -0.8*dycdt;
+        newstar->derivative[DERIVATIVE_STELLAR_CENTRAL_OXYGEN] = -0.2*dycdt;
+        // TODO: C12 + alpha -> O16 + gamma properly
+
+        /*
+         * Rate of change of convective core mass. We calculate this
+         * by calculating the core mass at a slightly higher Xc, hence
+         * we have dMc/dYc, then multiply by dYc/dt to calculate
+         * dMc/dt.
+         */
+        {
+            const double dYc = 1e-7;
+            Interpolate(stardata->store->MINT_tables[MINT_TABLE_CHeB],
+                        ((const double[]){
+                            log10(newstar->mass),
+                            newstar->mint->XHec + dYc
+                        }),
+                        result_conv,
+                        FALSE);
+            const double dmc =
+                newstar->mass * result_conv[MINT_CHeB_CONVECTIVE_CORE_MASS_OVERSHOOT] -
+                newstar->convective_core_mass;
+            newstar->convective_core_mass_time_derivative =
+                dmc/dYc * dycdt;
+        }
+
+        newstar->core_radius = 0.0;
+        newstar->convective_core_mass = newstar->mass * get_data(MINT_CHeB_CONVECTIVE_CORE_MASS_OVERSHOOT);
+        newstar->convective_core_radius = newstar->radius * get_data(MINT_CHeB_CONVECTIVE_CORE_RADIUS_OVERSHOOT);
+
+        // TODO : update
+
+        // helium core growth
+        // newstar->core_mass[CORE_He] = // assumed constant in CHeB
+
+        // assume CO core never shrinks
+        newstar->core_mass[CORE_CO] = Max(newstar->convective_core_mass,
+                                          newstar->core_mass[CORE_CO]);
+
+        newstar->menv = newstar->mass * get_data(MINT_CHeB_CONVECTIVE_ENVELOPE_MASS);
+        newstar->renv = newstar->radius * (1.0 - get_data(MINT_CHeB_CONVECTIVE_ENVELOPE_RADIUS));
+
+        newstar->mint->mean_molecular_weight_star = get_data(MINT_CHeB_MEAN_MOLECULAR_WEIGHT_AVERAGE);
+        newstar->mint->mean_molecular_weight_core = get_data(MINT_CHeB_MEAN_MOLECULAR_WEIGHT_CORE);
+
+        /*
+         * Set the phase start time if it hasn't already been set
+         */
+        if(newstar->phase_start[CHeB].time < 0.0)
+        {
+            newstar->phase_start[CHeB].time = stardata->model.time;
+        }
+
+        /*
+         * To calculate the "age" we use XHec as a (linear) proxy
+         */
+        double phase_lifetime;
+        MINT_CHeB_TA(stardata,
+                     newstar->phase_start_mass,
+                     &phase_lifetime,
+                     NULL);
+
+        const double helium_available = 1.0 - stardata->common.metallicity;
+        newstar->age =
+            newstar->phase_start[CHeB].time +
+            phase_lifetime * (1.0 - newstar->mint->XHec / helium_available );
+        newstar->epoch = stardata->model.time - newstar->age;
+        if(vb)printf("CHeB = %g + %g * (1 - %g / %g) = %g : EPOCH %g\n",
+                     newstar->phase_start[CHeB].time,
+                     phase_lifetime,
+                     newstar->mint->XHec,
+                     helium_available,
+                     newstar->age,
+                     newstar->epoch);
+
+
+        newstar->k2 = get_data(MINT_CHeB_K2);
+        newstar->E2 =
+            available[MINT_CHeB_TIDAL_E2] ?
+            get_data(MINT_CHeB_TIDAL_E2) :
+            E2(stardata,newstar);
+
+        newstar->E_Zahn = get_data(MINT_CHeB_TIDAL_E_FOR_LAMBDA);
+
+        if(newstar->mint->XHec < 0.9)
+        {
+            //Flexit;
+        }
+
+        if(newstar->mint->XHec > 0.1)
+        {
+            /*
+             * increase shell resolution around the burning core
+             */
+            MINT_resolve_core_boundary(stardata,
+                                       newstar,
+                                       newstar->convective_core_mass, // band centre
+                                       newstar->phase_start_mass * 0.005, // band thickness
+                                       1e-3, // thinnest shell dm
+                                       newstar->phase_start_mass / stardata->preferences->MINT_nshells, // default shell dm
+                                       TRUE);
+        }
+
+        /*
+         * Interpolate to get Chebyshev data for shells
+         */
+        MINT_CHeB_get_chebyshev_data(stardata,
+                                     newstar,
+                                     result_cheb);
+
+        /*
+         * Chebyshev coordinates
+         */
+        const double * const coords = RCHash_nest_get_entry(
+            stardata->store->MINT_headers[MINT_TABLE_CHeB]->hash,
+            "Chebyshev grid",
+            "masses"
+            )->value.value.double_array_data;
+
+        int n_cheb = 0;
+
+        Foreach_shell(newstar,shell)
+        {
+            /*
+             * Set stellar shell data from Chebyshev data
+             */
+            n_cheb = MINT_CHeB_set_chebyshev_data(stardata,
+                                                  newstar,
+                                                  shell,
+                                                  coords,
+                                                  result_cheb,
+                                                  n_cheb,
+                                                  map,
+                                                  available);
+
+            /*
+             * Identify instantly-mixed region from the
+             * convective core mass : note this needs
+             * updating to include the overshooting region
+             * and the surface convection zone
+             */
+            shell->convective = (shell->m < newstar->convective_core_mass) ? TRUE : FALSE;
+        } /* loop over shells */
+
+        break;
+    }
+
+
+    if(newstar->mint->XHec <= 0.0)
+    {
+        /*
+         * Finished CHeB
+         */
+        newstar->stellar_type = EAGB;
+        newstar->mint->XHec = 0.0;
+    }
+
+    if(vb)
+    {
+        printf("MINT CHeB star %d: stellar_type=%d (Z=%g, M=%g, XHec=%g) at t=%g = L=%g R=%g dYHec/dt=%g\n",
+               newstar->starnum,
+               newstar->stellar_type,
+               stardata->common.metallicity,
+               newstar->mass,
+               newstar->mint->XHec,
+               stardata->model.time,
+               newstar->luminosity,
+               newstar->radius,
+               newstar->derivative[DERIVATIVE_STELLAR_CENTRAL_HELIUM]);
+
+        printf("MINT CHeB time remaining %g Myr\n",
+               - 1e-6 * newstar->mint->XHec /newstar->derivative[DERIVATIVE_STELLAR_CENTRAL_HELIUM]
+            );
+
+    }
+
+#ifdef PRE_MAIN_SEQUENCE
+    if(stardata->preferences->pre_main_sequence == TRUE)
+    {
+        double f = preMS_radius_factor(newstar->mass, 1e6*newstar->age);
+        newstar->radius *= f;
+        newstar->PMS = Boolean_(f>1.00001);
+        Dprint("preMS f=%g r=%g roche_radius=%g roche_radius_at_periastron=%g PMS? %d\n",
+               f,
+               newstar->radius,
+               newstar->roche_radius,
+               newstar->roche_radius_at_periastron,
+               newstar->PMS);
+    }
+#endif
+
+    Dprint("MS radius %g\n",newstar->radius);
+
+#ifdef BSE
+    /*
+     * BSE-MINT compatibility
+     */
+#endif // BSE
+
+    Dprint("return stellar type %d, XHe C %g\n",
+           newstar->stellar_type,
+           newstar->mint->XHec);
+
+
+#endif // OLDCODE
 
     return newstar->stellar_type;
 }
diff --git a/src/MINT/MINT_stellar_structure_MS.c b/src/MINT/MINT_stellar_structure_MS.c
index 59979fbc89b53253547395cf3c7c9e3fdd4508f8..2196f0d90fcfac6ad0d62910b848acab8f3a5633 100644
--- a/src/MINT/MINT_stellar_structure_MS.c
+++ b/src/MINT/MINT_stellar_structure_MS.c
@@ -26,7 +26,6 @@ Stellar_type MINT_stellar_structure_MS(struct stardata_t * RESTRICT const starda
     /* MS stars have no He, CO or GB core */
     newstar->core_mass[CORE_He] = 0.0;
     newstar->core_mass[CORE_CO] = 0.0;
-    newstar->core_mass[CORE_He] = 0.0;
 
     newstar->derivative[DERIVATIVE_STELLAR_He_CORE_MASS] = 0.0;
     newstar->derivative[DERIVATIVE_STELLAR_GB_CORE_MASS] = 0.0;
@@ -150,7 +149,7 @@ Stellar_type MINT_stellar_structure_MS(struct stardata_t * RESTRICT const starda
     Loop_forever
     {
         /*
-         * Interpolate) to find stellar structure
+         * Interpolate to find stellar structure
          */
         Interpolate(stardata->store->MINT_tables[MINT_TABLE_MS],
                     ((double[]){
@@ -299,14 +298,10 @@ Stellar_type MINT_stellar_structure_MS(struct stardata_t * RESTRICT const starda
                      stardata->model.dt);
 
         newstar->k2 = get_data(MINT_MS_K2);
-        if(available[MINT_MS_TIDAL_E2])
-        {
-            newstar->E2 = get_data(MINT_MS_TIDAL_E2);
-        }
-        else
-        {
-            newstar->E2 = E2(stardata,newstar);
-        }
+        newstar->E2 =
+            available[MINT_MS_TIDAL_E2] ?
+            get_data(MINT_MS_TIDAL_E2) :
+            E2(stardata,newstar);
         newstar->E_Zahn = get_data(MINT_MS_TIDAL_E_FOR_LAMBDA);
 
         /*
@@ -357,7 +352,7 @@ Stellar_type MINT_stellar_structure_MS(struct stardata_t * RESTRICT const starda
     if(vb)
     {
 
-        if(0)printf("MINT%d(%g,%g,%g) at t=%g = L=%g R=%g dXHc/dt=%g\n",
+        if(0)printf("MINT MS %d(%g,%g,%g) at t=%g = L=%g R=%g dXHc/dt=%g\n",
                     newstar->starnum,
                     log10(stardata->common.metallicity),
                     newstar->phase_start_mass,
diff --git a/src/MINT/MINT_stellar_structure_functions.def b/src/MINT/MINT_stellar_structure_functions.def
index e11a22a393de3b4431084a1db0b2c4b9f48f8a4a..955100e9fdf73d5b238fc5dba365ae2ce1638541 100644
--- a/src/MINT/MINT_stellar_structure_functions.def
+++ b/src/MINT/MINT_stellar_structure_functions.def
@@ -28,7 +28,7 @@ enum {
         X(MAIN_SEQUENCE, &MINT_stellar_structure_MS, FALSE)             \
         X(HERZTSPRUNG_GAP, &MINT_stellar_structure_HG, AFTER)           \
         X(GIANT_BRANCH, &MINT_stellar_structure_GB, BEFORE)             \
-        X(CORE_HELIUM_BURNING, &MINT_stellar_structure_CHeB, BEFORE)    \
+        X(CORE_HELIUM_BURNING, &MINT_stellar_structure_CHeB, FALSE)     \
         X(EAGB, &MINT_stellar_structure_EAGB, BEFORE)                   \
         X(TPAGB, &MINT_stellar_structure_TPAGB, BEFORE)                 \
         X(HeMS, &MINT_stellar_structure_HeMS, BEFORE)                   \
diff --git a/src/binary_c_derivatives.def b/src/binary_c_derivatives.def
index 4678f9432659979ef2f667c9dae6e5fec5a32c21..ea0ef2d500ea12f18a4f25924f037170e4191e58 100644
--- a/src/binary_c_derivatives.def
+++ b/src/binary_c_derivatives.def
@@ -19,80 +19,81 @@
  */
 
 #define DERIVATIVES_STELLAR_LIST                                        \
-    X( MASS, "M (Total)", &star->mass, &derivative_check_stellar_mass,STELLAR )  \
-        X( MASS_WIND_LOSS, "M Wind Loss", NULL, NULL ,STELLAR )                  \
-        X( MASS_WIND_GAIN, "M Wind Gain", NULL, NULL ,STELLAR )                  \
-        X( RADIUS, "Radius", NULL, NULL ,STELLAR )                               \
-        X( LUMINOSITY,"Luminosity", NULL, NULL  ,STELLAR )                       \
-        X( CORE_MASS, "Core Mass", NULL, NULL ,STELLAR )                         \
-        X( TEMPERATURE, "Stellar temperature", NULL, NULL ,STELLAR )             \
-        X( MASS_RLOF_LOSS,"M RLOF Loss", NULL, NULL  ,STELLAR )                  \
-        X( MASS_RLOF_GAIN, "M RLOF Gain", NULL, NULL ,STELLAR )                  \
-        X( MASS_NOVA, "M Nova", NULL, NULL ,STELLAR )                            \
-        X( MASS_RLOF_TRANSFER, "M RLOF Transfer", NULL, NULL ,STELLAR )          \
-        X( MASS_DISC_GAIN, "M Disc Gain", NULL, NULL ,STELLAR )                  \
-        X( MASS_DISC_LOSS, "M Disc Loss", NULL, NULL ,STELLAR )                  \
-        X( MASS_DECRETION_DISC , "M Decretion Disc Loss", NULL, NULL,STELLAR )   \
-        X( MASS_CBDISC_GAIN , "M CBdisc Gain", NULL, NULL,STELLAR )              \
-        X( MASS_NONCONSERVATIVE_LOSS , "M Non-conservative Loss", NULL, NULL ,STELLAR ) \
-        X( MASS_IRRADIATIVE_LOSS, "M Irradiative Loss", NULL, NULL ,STELLAR )    \
-        X( MASS_ARTIFICIAL, "M artificial", NULL, NULL ,STELLAR )                \
-        X( ANGMOM_RLOF_LOSS, "J RLOF Loss", NULL, NULL ,STELLAR )                \
-        X( ANGMOM_RLOF_GAIN, "J RLOF Gain", NULL, NULL,STELLAR )                 \
-        X( ANGMOM_WIND_LOSS, "J Wind Loss", NULL, NULL ,STELLAR )                \
-        X( ANGMOM_WIND_GAIN, "J Wind Gain", NULL, NULL ,STELLAR )                \
-        X( ANGMOM_TIDES, "J Tides", NULL, NULL ,STELLAR )                        \
-        X( ANGMOM_MAGNETIC_BRAKING,"J Magnetic Braking", NULL, NULL  ,STELLAR )  \
-        X( ANGMOM_DECRETION_DISC, "J Decretion Disc", NULL, NULL ,STELLAR )      \
-        X( ANGMOM_CBDISC_GAIN, "J CBdisc gain", NULL, NULL ,STELLAR )            \
-        X( ANGMOM_NOVA, "J Nova", NULL, NULL ,STELLAR )                          \
-        X( ANGMOM_NONCONSERVATIVE_LOSS, "J non-conservative", NULL, NULL ,STELLAR ) \
-        X( ANGMOM_ARTIFICIAL,"J artificial", NULL, NULL  ,STELLAR )              \
-        X( ANGMOM_CORE_ENVELOPE_COUPLING, "J core-envelope coupling", NULL, NULL ,STELLAR ) \
-        X( ANGMOM, "J (Total)", &star->angular_momentum, &derivative_check_stellar_angmom ,STELLAR ) \
-        X( ECCENTRICITY,  "eccentricity", NULL, NULL,STELLAR )                   \
-        X( ANGULAR_VELOCITY_TIDES, "angular velocity (tides)", NULL, NULL ,STELLAR ) \
-        X( He_CORE_MASS, "He core mass", &star->core_mass[CORE_He], &derivative_check_stellar_he_core_mass ,STELLAR ) \
-        X( GB_CORE_MASS, "GB core mass" , NULL, NULL,STELLAR )                   \
-        X( CO_CORE_MASS, "CO core mass", &star->core_mass[CORE_CO], &derivative_check_stellar_CO_core_mass,STELLAR ) \
-        X( ONe_CORE_MASS, "ONe core mass", &star->core_mass[CORE_ONe], NULL, STELLAR ) \
-        X( Si_CORE_MASS, "Si core mass", &star->core_mass[CORE_Si], NULL, STELLAR ) \
-        X( Fe_CORE_MASS, "Fe core mass", &star->core_mass[CORE_Fe], NULL, STELLAR ) \
-        X( NEUTRON_CORE_MASS, "Neutron core mass", &star->core_mass[CORE_NEUTRON], NULL, STELLAR ) \
-        X( BH_CORE_MASS, "Black-hole core mass", &star->core_mass[CORE_BLACK_HOLE], NULL, STELLAR ) \
-        X( CENTRAL_HYDROGEN, "Central hydrogen", &star->mint->XHc, &derivative_check_stellar_central_hydrogen ,STELLAR ) \
-        X( CENTRAL_HELIUM, "Central helium", &star->mint->XHec, &derivative_check_stellar_central_helium ,STELLAR ) \
-        X( CENTRAL_CARBON,"Central carbon", NULL, NULL  ,STELLAR )               \
-        X( He_CORE_MASS_NO_TDUP, "He core mass (no 3DUP)", &star->core_mass_no_3dup, NULL ,STELLAR ) \
-        X( NUM_THERMAL_PULSES, "Number of thermal pulses", &star->num_thermal_pulses, &derivative_check_num_thermal_pulses,STELLAR ) \
-        X( NUM_THERMAL_PULSES_SINCE_MCMIN, "Number of thermal pulses since mcmin", &star->num_thermal_pulses_since_mcmin, &derivative_check_num_thermal_pulses,STELLAR ) \
-        X( H_LAYER_MASS,"nova H layer mass", &star->dm_novaH, &derivative_check_dm_novaH,STELLAR ) \
-        X( ROCHE_RADIUS, "Roche radius", NULL, NULL ,STELLAR )                   \
-        X( NUM_NOVAE,"Number of novae", &star->num_novae, &derivative_check_num_novae,STELLAR ) \
-        X (NUMBER, "Number", NULL, NULL,STELLAR )
+    X( MASS, "M (Total)", &star->mass, &derivative_check_stellar_mass,STELLAR ) \
+        X( MASS_WIND_LOSS, "M Wind Loss", NULL, NULL ,STELLAR )         \
+        X( MASS_WIND_GAIN, "M Wind Gain", NULL, NULL ,STELLAR )         \
+        X( RADIUS, "Radius", NULL, NULL ,STELLAR )                      \
+        X( LUMINOSITY,"Luminosity", NULL, NULL  ,STELLAR )              \
+        X( CORE_MASS, "Core Mass", NULL, NULL ,STELLAR )                \
+    X( TEMPERATURE, "Stellar temperature", NULL, NULL ,STELLAR )        \
+    X( MASS_RLOF_LOSS,"M RLOF Loss", NULL, NULL  ,STELLAR )             \
+    X( MASS_RLOF_GAIN, "M RLOF Gain", NULL, NULL ,STELLAR )             \
+    X( MASS_NOVA, "M Nova", NULL, NULL ,STELLAR )                       \
+    X( MASS_RLOF_TRANSFER, "M RLOF Transfer", NULL, NULL ,STELLAR )     \
+    X( MASS_DISC_GAIN, "M Disc Gain", NULL, NULL ,STELLAR )             \
+    X( MASS_DISC_LOSS, "M Disc Loss", NULL, NULL ,STELLAR )             \
+    X( MASS_DECRETION_DISC , "M Decretion Disc Loss", NULL, NULL,STELLAR ) \
+    X( MASS_CBDISC_GAIN , "M CBdisc Gain", NULL, NULL,STELLAR )         \
+    X( MASS_NONCONSERVATIVE_LOSS , "M Non-conservative Loss", NULL, NULL ,STELLAR ) \
+    X( MASS_IRRADIATIVE_LOSS, "M Irradiative Loss", NULL, NULL ,STELLAR ) \
+    X( MASS_ARTIFICIAL, "M artificial", NULL, NULL ,STELLAR )           \
+    X( ANGMOM_RLOF_LOSS, "J RLOF Loss", NULL, NULL ,STELLAR )           \
+    X( ANGMOM_RLOF_GAIN, "J RLOF Gain", NULL, NULL,STELLAR )            \
+    X( ANGMOM_WIND_LOSS, "J Wind Loss", NULL, NULL ,STELLAR )           \
+    X( ANGMOM_WIND_GAIN, "J Wind Gain", NULL, NULL ,STELLAR )           \
+    X( ANGMOM_TIDES, "J Tides", NULL, NULL ,STELLAR )                   \
+    X( ANGMOM_MAGNETIC_BRAKING,"J Magnetic Braking", NULL, NULL  ,STELLAR ) \
+    X( ANGMOM_DECRETION_DISC, "J Decretion Disc", NULL, NULL ,STELLAR ) \
+    X( ANGMOM_CBDISC_GAIN, "J CBdisc gain", NULL, NULL ,STELLAR )       \
+    X( ANGMOM_NOVA, "J Nova", NULL, NULL ,STELLAR )                     \
+    X( ANGMOM_NONCONSERVATIVE_LOSS, "J non-conservative", NULL, NULL ,STELLAR ) \
+    X( ANGMOM_ARTIFICIAL,"J artificial", NULL, NULL  ,STELLAR )         \
+    X( ANGMOM_CORE_ENVELOPE_COUPLING, "J core-envelope coupling", NULL, NULL ,STELLAR ) \
+    X( ANGMOM, "J (Total)", &star->angular_momentum, &derivative_check_stellar_angmom ,STELLAR ) \
+    X( ECCENTRICITY,  "eccentricity", NULL, NULL,STELLAR )              \
+    X( ANGULAR_VELOCITY_TIDES, "angular velocity (tides)", NULL, NULL ,STELLAR ) \
+    X( He_CORE_MASS, "He core mass", &star->core_mass[CORE_He], &derivative_check_stellar_he_core_mass ,STELLAR ) \
+    X( GB_CORE_MASS, "GB core mass" , NULL, NULL,STELLAR )              \
+    X( CO_CORE_MASS, "CO core mass", &star->core_mass[CORE_CO], &derivative_check_stellar_CO_core_mass,STELLAR ) \
+    X( ONe_CORE_MASS, "ONe core mass", &star->core_mass[CORE_ONe], NULL, STELLAR ) \
+    X( Si_CORE_MASS, "Si core mass", &star->core_mass[CORE_Si], NULL, STELLAR ) \
+    X( Fe_CORE_MASS, "Fe core mass", &star->core_mass[CORE_Fe], NULL, STELLAR ) \
+    X( NEUTRON_CORE_MASS, "Neutron core mass", &star->core_mass[CORE_NEUTRON], NULL, STELLAR ) \
+    X( BH_CORE_MASS, "Black-hole core mass", &star->core_mass[CORE_BLACK_HOLE], NULL, STELLAR ) \
+    X( CENTRAL_HYDROGEN, "Central hydrogen", &star->mint->XHc, &derivative_check_stellar_central_hydrogen ,STELLAR ) \
+    X( CENTRAL_HELIUM, "Central helium", &star->mint->XHec, &derivative_check_stellar_central_helium ,STELLAR ) \
+    X( CENTRAL_CARBON,"Central carbon", &star->mint->XCc, NULL  ,STELLAR )          \
+    X( CENTRAL_OXYGEN,"Central oxygen", &star->mint->XOc, NULL  ,STELLAR )          \
+    X( He_CORE_MASS_NO_TDUP, "He core mass (no 3DUP)", &star->core_mass_no_3dup, NULL ,STELLAR ) \
+    X( NUM_THERMAL_PULSES, "Number of thermal pulses", &star->num_thermal_pulses, &derivative_check_num_thermal_pulses,STELLAR ) \
+    X( NUM_THERMAL_PULSES_SINCE_MCMIN, "Number of thermal pulses since mcmin", &star->num_thermal_pulses_since_mcmin, &derivative_check_num_thermal_pulses,STELLAR ) \
+    X( H_LAYER_MASS,"nova H layer mass", &star->dm_novaH, &derivative_check_dm_novaH,STELLAR ) \
+    X( ROCHE_RADIUS, "Roche radius", NULL, NULL ,STELLAR )              \
+    X( NUM_NOVAE,"Number of novae", &star->num_novae, &derivative_check_num_novae,STELLAR ) \
+    X (NUMBER, "Number", NULL, NULL,STELLAR )
 
 #define SYSTEM_DERIVATIVES_LIST                                         \
     X(ORBIT_ANGMOM , "J (total)", &stardata->common.orbit.angular_momentum, NULL, SYSTEM ) \
-        X(ORBIT_SEMI_MAJOR_AXIS, "a", NULL, NULL, SYSTEM )                       \
+        X(ORBIT_SEMI_MAJOR_AXIS, "a", NULL, NULL, SYSTEM )              \
         X(ORBIT_ECCENTRICITY, "e (total)", &stardata->common.orbit.eccentricity, NULL, SYSTEM ) \
         X(ORBIT_ANGMOM_GRAVITATIONAL_RADIATION, "J Gravitational Radiation", NULL, NULL, SYSTEM ) \
-        X(ORBIT_ANGMOM_WIND_LOSS, "J Wind Loss", NULL, NULL, SYSTEM )            \
-        X(ORBIT_ANGMOM_WIND_GAIN, "J Wind Gain", NULL, NULL, SYSTEM )            \
-        X(ORBIT_ANGMOM_RLOF_LOSS, "J RLOF Loss", NULL, NULL, SYSTEM )            \
-        X(ORBIT_ANGMOM_RLOF_GAIN, "J RLOF Gain", NULL, NULL, SYSTEM )            \
-        X(ORBIT_ANGMOM_CBDISC, "J Circumbinary Disc", NULL, NULL, SYSTEM )       \
-        X(ORBIT_ANGMOM_TIDES, "J tides", NULL, NULL, SYSTEM )                    \
-        X(ORBIT_ANGMOM_NONCONSERVATIVE_LOSS, "J nonconservative", NULL, NULL, SYSTEM ) \
-        X(ORBIT_ANGMOM_NOVA, "J nova", NULL, NULL, SYSTEM )                      \
-        X(ORBIT_ANGMOM_ORBITING_OBJECTS, "J orbiting objects", NULL, NULL, SYSTEM ) \
-        X(ORBIT_ANGMOM_ARTIFICIAL , "J artificial", NULL, NULL, SYSTEM )         \
-        X(ORBIT_ECCENTRICITY_GRAVITATIONAL_RADIATION , "e Gravitational Radiation", NULL, NULL, SYSTEM ) \
-        X(ORBIT_ECCENTRICITY_TIDES, "e Tides", NULL, NULL, SYSTEM )              \
-        X(ORBIT_ECCENTRICITY_WINDS , "e Winds", NULL, NULL, SYSTEM )             \
-        X(ORBIT_ECCENTRICITY_CBDISC , "e Circumbinary Disc", NULL, NULL, SYSTEM ) \
-        X(SYSTEM_CBDISC_MASS_LOSS, "Mdot Circumbinary Disc", NULL, NULL, SYSTEM ) \
-        X(SYSTEM_TEST , "Test", &stardata->common.test, NULL, SYSTEM )           \
-        X(SYSTEM_NUMBER, "Number", NULL, NULL, SYSTEM )
+        X(ORBIT_ANGMOM_WIND_LOSS, "J Wind Loss", NULL, NULL, SYSTEM )   \
+        X(ORBIT_ANGMOM_WIND_GAIN, "J Wind Gain", NULL, NULL, SYSTEM )   \
+    X(ORBIT_ANGMOM_RLOF_LOSS, "J RLOF Loss", NULL, NULL, SYSTEM )       \
+    X(ORBIT_ANGMOM_RLOF_GAIN, "J RLOF Gain", NULL, NULL, SYSTEM )       \
+    X(ORBIT_ANGMOM_CBDISC, "J Circumbinary Disc", NULL, NULL, SYSTEM )  \
+    X(ORBIT_ANGMOM_TIDES, "J tides", NULL, NULL, SYSTEM )               \
+    X(ORBIT_ANGMOM_NONCONSERVATIVE_LOSS, "J nonconservative", NULL, NULL, SYSTEM ) \
+    X(ORBIT_ANGMOM_NOVA, "J nova", NULL, NULL, SYSTEM )                 \
+    X(ORBIT_ANGMOM_ORBITING_OBJECTS, "J orbiting objects", NULL, NULL, SYSTEM ) \
+    X(ORBIT_ANGMOM_ARTIFICIAL , "J artificial", NULL, NULL, SYSTEM )    \
+    X(ORBIT_ECCENTRICITY_GRAVITATIONAL_RADIATION , "e Gravitational Radiation", NULL, NULL, SYSTEM ) \
+    X(ORBIT_ECCENTRICITY_TIDES, "e Tides", NULL, NULL, SYSTEM )         \
+    X(ORBIT_ECCENTRICITY_WINDS , "e Winds", NULL, NULL, SYSTEM )        \
+    X(ORBIT_ECCENTRICITY_CBDISC , "e Circumbinary Disc", NULL, NULL, SYSTEM ) \
+    X(SYSTEM_CBDISC_MASS_LOSS, "Mdot Circumbinary Disc", NULL, NULL, SYSTEM ) \
+    X(SYSTEM_TEST , "Test", &stardata->common.test, NULL, SYSTEM )      \
+    X(SYSTEM_NUMBER, "Number", NULL, NULL, SYSTEM )
 
 #define ORBITING_OBJECT_DERIVATIVES_LIST                                \
     X(MASS, "Mass", NULL, NULL, ORBITING_OBJECT)                        \
diff --git a/src/binary_c_structures.h b/src/binary_c_structures.h
index e1295cbbd55d42ac8d103a74d975dbcc98f95dba..6130d5bc995c59c3d6d427323d0eb6eaa2abf3cc 100644
--- a/src/binary_c_structures.h
+++ b/src/binary_c_structures.h
@@ -496,6 +496,9 @@ struct tmpstore_t {
     double * MINT_MS_result;
     double * MINT_MS_result_cheb;
     double * MINT_MS_result_conv;
+    double * MINT_CHeB_result;
+    double * MINT_CHeB_result_cheb;
+    double * MINT_CHeB_result_conv;
 #endif//MINT
     unsigned int kaps_rentrop_GSL_n;
     gsl_matrix * kaps_rentrop_GSL_m;
@@ -1208,6 +1211,7 @@ struct preferences_t {
     int MINT_Kippenhahn;
     Stellar_type MINT_Kippenhahn_stellar_type;
     Stellar_type MINT_Kippenhahn_companion_stellar_type;
+    int MINT_nshells;
 #endif // MINT
 
     int float_overflow_checks;
@@ -1557,7 +1561,7 @@ struct mint_shell_t {
     double radius Rchash_var; /* r */
     double gamma1 Rchash_var; /* gamma1 (from equation of state) */
     double pressure_scale_height Rchash_var; /* -dlnr/dlnP */
-    Boolean convective Rchash_var;
+    Boolean convective Rchash_var; /* TRUE if convective */
 };
 
 struct mint_t {
@@ -1567,8 +1571,8 @@ struct mint_t {
      */
     Abundance XHc Rchash_var; /* central hydrogen abundance */
     Abundance XHec Rchash_var; /* central helium abundance */
-    Abundance XCc Rchash_var; /* central helium abundance */
-    Abundance XOc Rchash_var; /* central helium abundance */
+    Abundance XCc Rchash_var; /* central carbon abundance */
+    Abundance XOc Rchash_var; /* central oxygen abundance */
 
     /*
      * molecular weights of the whole star
diff --git a/src/binary_star_functions/binary_star_functions_prototypes.h b/src/binary_star_functions/binary_star_functions_prototypes.h
index 0c3230af46fc62b8041bf905ec1f75375ee3954b..f5607fce25e74c88c179f7c6daabbfde77d941b7 100644
--- a/src/binary_star_functions/binary_star_functions_prototypes.h
+++ b/src/binary_star_functions/binary_star_functions_prototypes.h
@@ -359,16 +359,6 @@ void merge_cores(const struct star_t * const a,
 double accretion_radius(struct stardata_t * const stardata,
                         struct star_t * const star);
 
-Pure_function
-double Hut_f5(const double ecc);
-Pure_function
-double Hut_f4(const double ecc);
-Pure_function
-double Hut_f3(const double ecc);
-Pure_function
-double Hut_f2(const double ecc);
-Pure_function
-double Hut_f1(const double ecc);
 
 double wind_orbital_angular_momentum_derivative(
     const struct stardata_t * RESTRICT const stardata,
diff --git a/src/logging/log_every_timestep.c b/src/logging/log_every_timestep.c
index a2bb49119e26b4178f01e7d3b9c870c402c67668..e4420343b66d11f5aaeae551e24cfc1c213c928a 100644
--- a/src/logging/log_every_timestep.c
+++ b/src/logging/log_every_timestep.c
@@ -3054,4 +3054,5 @@ void log_every_timestep(struct stardata_t * RESTRICT const stardata)
 #ifdef ORBITING_OBJECTS
     log_orbiting_objects(stardata);
 #endif// ORBITING_OBJECTS
+
 }
diff --git a/src/memory/free_tmpstore.c b/src/memory/free_tmpstore.c
index d807888daf7406b06cdc074419a81aeee562ea8a..7322c337fcf46e610ab482e519c02fe8c6a4415b 100644
--- a/src/memory/free_tmpstore.c
+++ b/src/memory/free_tmpstore.c
@@ -21,6 +21,9 @@ void free_tmpstore(struct tmpstore_t * RESTRICT const tmpstore,
         Tmp_free(MINT_MS_result);
         Tmp_free(MINT_MS_result_cheb);
         Tmp_free(MINT_MS_result_conv);
+        Tmp_free(MINT_CHeB_result);
+        Tmp_free(MINT_CHeB_result_cheb);
+        Tmp_free(MINT_CHeB_result_conv);
 #endif // MINT
         Tmp_free(random_system_argstring);
         Tmp_free(unresolved_magnitudes);
diff --git a/src/nucsyn/nucsyn_burn_helium.c b/src/nucsyn/nucsyn_burn_helium.c
index e0ba2a650e3f0ca3230a3b3b211711beec437774..2d7695d00dc901bbf17cb975eaab9375235f988f 100644
--- a/src/nucsyn/nucsyn_burn_helium.c
+++ b/src/nucsyn/nucsyn_burn_helium.c
@@ -69,6 +69,9 @@ double nucsyn_burn_helium(struct stardata_t * RESTRICT const stardata,
     return err;
 }
 
+#define HE_FAC1 (3.0/6.0)
+#define HE_FAC2 (1.0/6.0)
+
 static void helium_derivatives(const double * RESTRICT const y,
                                double * RESTRICT const dydt,
                                const double * RESTRICT const sigmav,
@@ -79,20 +82,20 @@ static void helium_derivatives(const double * RESTRICT const y,
     const double c = sigmav[SIGMAV_O16_He4] * y[3] * y[1];
     const double d = sigmav[SIGMAV_Ne20_He4] * y[4] * y[1];
 
-    dydt[1] = -3.0*a -b -c -d;
-    dydt[2] = a - b;
+    dydt[1] = -a*HE_FAC1 -b -c -d;
+    dydt[2] = a*HE_FAC2 - b;
     dydt[3] = b - c;
     dydt[4] = c - d;
     dydt[5] = d;
 
     /*
     dydt[1] =
-        -3.0*sigmav[SIGMAV_3He4] * Pow3(N[XHe4])
+        -1/2*sigmav[SIGMAV_3He4] * Pow3(N[XHe4])
         -sigmav[SIGMAV_C12_He4] * N[XC12] * N[XHe4]
         -sigmav[SIGMAV_O16_He4] * N[XO16] * N[XHe4];
         -sigmav[SIGMAV_Ne20_He4] * N[XNe20] * N[XHe4];
     dydt[2] =
-        + sigmav[SIGMAV_3He4] * Pow3(N[XHe4])
+        +1/6*sigmav[SIGMAV_3He4] * Pow3(N[XHe4])
         - sigmav[SIGMAV_C12_He4] * N[XC12] * N[XHe4];
     dydt[3] =
         + sigmav[SIGMAV_C12_He4] * N[XC12] * N[XHe4]
@@ -109,7 +112,7 @@ static void helium_Jacobian(double *const*const RESTRICT J,
                             const double * RESTRICT const N,
                             const double * RESTRICT const sigmav)
 {
-    const double dummy0 = 3.0 * Pow2(N[XHe4]) * sigmav[SIGMAV_3He4];
+    const double dummy0 = sigmav[SIGMAV_3He4] * 3.0 * Pow2(N[XHe4]) ;
     const double dummy1 = sigmav[SIGMAV_C12_He4] * N[XC12];
     const double dummy2 = sigmav[SIGMAV_O16_He4] * N[XO16];
     const double dummy3 = sigmav[SIGMAV_Ne20_He4] * N[XNe20];
@@ -118,7 +121,7 @@ static void helium_Jacobian(double *const*const RESTRICT J,
     const double dummy6 = sigmav[SIGMAV_Ne20_He4] * N[XHe4];
 
     J[1][1] =
-        -3.0 * dummy0
+        -dummy0*HE_FAC1
         -dummy1
         -dummy2
         -dummy3;
@@ -128,7 +131,7 @@ static void helium_Jacobian(double *const*const RESTRICT J,
     J[1][5] = 0.0;
 
     J[2][1] =
-        + dummy0
+        + dummy0*HE_FAC2
         - dummy1;
     J[2][2] = - dummy4;
     J[2][3] = 0.0;
diff --git a/src/nucsyn/nucsyn_hbb.c b/src/nucsyn/nucsyn_hbb.c
index ecb6b35ed77c23b027cafcb1dea5510e810ba269..6c11183963df7b7786b9578d52fee9560218d489 100644
--- a/src/nucsyn/nucsyn_hbb.c
+++ b/src/nucsyn/nucsyn_hbb.c
@@ -35,11 +35,12 @@ void nucsyn_hbb(const double T, /* (log10) temperature */
     Number_density * N = New_isotope_array;
 
     const double timemax = 100;
-    const double dtime = timemax/100.0;
-    const double logtemp = log10(1e8);
-    const double testrho = 1e3;
+    const double dtime = timemax/1000.0;
+    const double logtemp = log10(3e8);
+    const double testrho = 1e2;
 
     /* start with a lot of helium */
+    Clear_isotope_array(X);
     X[XHe4] = 1.0;
     X_to_N(imnuc,testrho,N,X,ISOTOPE_ARRAY_SIZE);
 
@@ -48,6 +49,7 @@ void nucsyn_hbb(const double T, /* (log10) temperature */
     {
         nucsyn_burn(logtemp,testrho,dtime,N,star,stardata);
         N_to_X(mnuc,testrho,N,X,ISOTOPE_ARRAY_SIZE);
+        time += dtime;
         printf("HEBURN %g %g %g %g %g %g %g\n",
                time,
                X[XHe4],
@@ -57,7 +59,6 @@ void nucsyn_hbb(const double T, /* (log10) temperature */
                X[XMg24],
                XXsum(stardata,X)
             );
-        time += dtime;
     }
     Safe_free(X);
     Safe_free(N);
diff --git a/src/nucsyn/nucsyn_prototypes.h b/src/nucsyn/nucsyn_prototypes.h
index 4e9a3e45c202db18b3bf33c92b423bff2f32d152..bee1dbb1e931e96c8609f514e9cf23a78a09d3ff 100644
--- a/src/nucsyn/nucsyn_prototypes.h
+++ b/src/nucsyn/nucsyn_prototypes.h
@@ -729,6 +729,21 @@ double Hot_function Pure_function nucsyn_nucleon_density(
 void nucsyn_burn_CNO_completely(Abundance * const X);
 
 
+Constant_function double nucsyn_reaclib(
+    const double t9, // temperature (/1e9)
+    const double i9, // 1.0/t9
+    const double tm13, // pow(t9,-0.3333333333)
+    const double t913, // pow(t9,0.333333333333333333333333)
+    const double t953, // pow(t9,1.66666666666666666666)
+    const double lt9, // (natural) log(t9)
+    const double ra, // reaclib parameters
+    const double rb,
+    const double rc,
+    const double rd,
+    const double re,
+    const double rf,
+    const double rg);
+
 
 #endif /* NUCSYN */
 #endif /* NUCSYN_PROTOYPES */
diff --git a/src/nucsyn/nucsyn_reaclib.c b/src/nucsyn/nucsyn_reaclib.c
index 66bfccc457b16f5da780897bdec65fea2120ae58..1da3839d5d81609afb424160402554489920446e 100644
--- a/src/nucsyn/nucsyn_reaclib.c
+++ b/src/nucsyn/nucsyn_reaclib.c
@@ -2,33 +2,20 @@
 No_empty_translation_unit_warning;
 
 
-Constant_function inline double nucsyn_reaclib(const double  t9, // temperature (/1e9)
-				 const double  i9, // 1.0/t9
-				 const double  tm13, // pow(t9,-0.3333333333)
-				 const double  t913, // pow(t9,0.333333333333333333333333)
-				 const double  t953, // pow(t9,1.66666666666666666666)
-				 const double  lt9, // (natural) log(t9)
-				 const double  ra, // reaclib parameters
-				 const double  rb,
-				 const double  rc,
-				 const double  rd,
-				 const double  re,
-				 const double  rf,
-				 const double  rg) ;
-
-Constant_function inline double nucsyn_reaclib(const double  t9, // temperature (/1e9)
-                                        const double  i9, // 1.0/t9
-                                        const double  tm13, // pow(t9,-0.3333333333)
-                                        const double  t913, // pow(t9,0.333333333333333333333333)
-                                        const double  t953, // pow(t9,1.66666666666666666666)
-                                        const double  lt9, // (natural) log(t9)
-                                        const double  ra, // reaclib parameters
-                                        const double   rb,
-                                        const double   rc,
-                                        const double   rd,
-                                        const double   re,
-                                        const double   rf,
-                                        const double   rg) 
+Constant_function double nucsyn_reaclib(
+    const double t9, // temperature (/1e9)
+    const double i9, // 1.0/t9
+    const double tm13, // pow(t9,-0.3333333333)
+    const double t913, // pow(t9,0.333333333333333333333333)
+    const double t953, // pow(t9,1.66666666666666666666)
+    const double lt9, // (natural) log(t9)
+    const double ra, // reaclib parameters
+    const double rb,
+    const double rc,
+    const double rd,
+    const double re,
+    const double rf,
+    const double rg)
 {
     /* function to evaluate reaclib-style cross sections */
     return exp(ra+
diff --git a/src/nucsyn/nucsyn_set_sigmav.c b/src/nucsyn/nucsyn_set_sigmav.c
index d179d6ba5e3eb7906e14a67fd01779bbb9112e39..f621224232cab26f171cd2ec0cbca551481d4168 100644
--- a/src/nucsyn/nucsyn_set_sigmav.c
+++ b/src/nucsyn/nucsyn_set_sigmav.c
@@ -10,6 +10,8 @@ No_empty_translation_unit_warning;
  */
 //#define REACTION_RATES_TABLE
 
+#define Reaclib(...) nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,__VA_ARGS__)
+
 void nucsyn_set_sigmav(struct stardata_t * const stardata,
                        struct store_t * const store,
                        const double temp, /* log temperature */
@@ -68,17 +70,19 @@ void nucsyn_set_sigmav(struct stardata_t * const stardata,
          * different then don't reset the sigmav array
          */
 
-
-/* speed advances (not pretty, but functional) */
-        const double  tm13 = 1.0/cbrt(t9); // t9^(-1/3)
-        const double  t23 = Pow2(tm13);    // t9^(-2/3)
-        const double  t15 = 1.0/Pow1p5(t9);// t9^(-3/2)
-        const double  t92 = Pow2(t9);      // t9^2
-        const double  t93 = t92*t9;        // t9^3
-        const double  t94 = t93*t9;        // t9^4
-        const double  t95 = t94*t9;        // t9^5
-        //const double  logt9 = log10(t9);
-        const double  i9 = 1.00/t9;        // t9^-1
+        /* speed advances (not pretty, but functional) */
+        const double tm13 = 1.0/cbrt(t9); // t9^(-1/3)
+        const double t23 = Pow2(tm13);    // t9^(-2/3)
+        const double t15 = 1.0/Pow1p5(t9);// t9^(-3/2)
+        const double t92 = Pow2(t9);      // t9^2
+        const double t93 = t92*t9;        // t9^3
+        const double t94 = t93*t9;        // t9^4
+        const double t95 = t94*t9;        // t9^5
+        //const double logt9 = log10(t9);
+        const double i9 = 1.00/t9;        // t9^-1
+        const double t953 = pow(t9,5.0/3.0);
+        const double lt9 = log(t9); // ln(T9)
+        const double t913 = cbrt(t9);
 
 #ifdef RATES_OF_AREND_JAN
         /* Sorry for the fortran variables */
@@ -115,9 +119,9 @@ void nucsyn_set_sigmav(struct stardata_t * const stardata,
 #endif // RATES_OF_AREND_JAN
 
 #ifdef RATES_OF_MARIA
-        const double t913 = cbrt(t9);
-        const double t953 = pow(t9,5.0/3.0);
-        const double lt9 = log(t9);
+
+
+
 #endif//RATES_OF_MARIA
 
 #ifdef RATES_OF_AMANDA
@@ -614,7 +618,10 @@ void nucsyn_set_sigmav(struct stardata_t * const stardata,
         /*
          * Maria's REACTION Ne20pgNa21 -> SIGMAV_T20
          */
-        sigmav[SIGMAV_T20]=nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.112078e+03,-0.705964e-01,-0.300779e+01,-0.149812e+03,0.893819e+02,-0.929669e+02,0.270204e+02)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.247838e+01,-0.691248,-0.409804e+01,0.119259e+02,-0.954890e+01,-0.511767,0.593970e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.863238e+01,-0.836539e+01,0.451031,-0.300658,-0.312270,0.314422e-01,0.148890e+01);
+        sigmav[SIGMAV_T20]=
+            Reaclib(0.112078e+03,-0.705964e-01,-0.300779e+01,-0.149812e+03,0.893819e+02,-0.929669e+02,0.270204e+02)+
+            Reaclib(0.247838e+01,-0.691248,-0.409804e+01,0.119259e+02,-0.954890e+01,-0.511767,0.593970e+01)+
+            Reaclib(0.863238e+01,-0.836539e+01,0.451031,-0.300658,-0.312270,0.314422e-01,0.148890e+01);
 
 #elif (defined(RATES_OF_AREND_JAN))
         sigmav[SIGMAV_T20]=(9.55e+06*exp(-19.447/T913)/(Pow2(T9)*Pow2(1.+0.0127/T923))
@@ -652,7 +659,8 @@ void nucsyn_set_sigmav(struct stardata_t * const stardata,
 
 #elif (defined(RATES_OF_MARIA))
         /* Maria's REACTION Ne21pgNa22 -> SIGMAV_T21 */
-        sigmav[SIGMAV_T21]=nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.296371e+02,-0.179861e-01,-0.163786e+02,-0.153383e+02,0.300382e+01,-0.900470,0.345694e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.437079e+02,-0.192696,0.765401e-02,-0.280351e-01,0.305416e-02,-0.269857e-03,-0.149092e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.114485e+02,-0.109787e+01,0.497670e-02,-0.154510e-01,0.147205e-02,-0.119000e-03,-0.149456e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.181545e+01,-0.140335e+01,-0.454476e-02,0.134943e-01,-0.123631e-02,0.971986e-04,-0.150486e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.768386e+01,-0.300261e+01,-0.713715,0.290262e+01,-0.260861,0.180284e-01,-0.237968e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.305967,-0.540027e+01,-0.125000e+02,0.257126e+02,-0.154104e+01,0.780952e-01,-0.118132e+02);
+        sigmav[SIGMAV_T21]=
+            Reaclib(0.296371e+02,-0.179861e-01,-0.163786e+02,-0.153383e+02,0.300382e+01,-0.900470,0.345694e+01)+Reaclib(-0.437079e+02,-0.192696,0.765401e-02,-0.280351e-01,0.305416e-02,-0.269857e-03,-0.149092e+01)+Reaclib(-0.114485e+02,-0.109787e+01,0.497670e-02,-0.154510e-01,0.147205e-02,-0.119000e-03,-0.149456e+01)+Reaclib(0.181545e+01,-0.140335e+01,-0.454476e-02,0.134943e-01,-0.123631e-02,0.971986e-04,-0.150486e+01)+Reaclib(0.768386e+01,-0.300261e+01,-0.713715,0.290262e+01,-0.260861,0.180284e-01,-0.237968e+01)+Reaclib(0.305967,-0.540027e+01,-0.125000e+02,0.257126e+02,-0.154104e+01,0.780952e-01,-0.118132e+02);
 
 #elif (defined(RATES_OF_AREND_JAN))
         sigmav[SIGMAV_T21]=(4.37e+08/T923*exp(-19.462/T913)
@@ -694,7 +702,7 @@ void nucsyn_set_sigmav(struct stardata_t * const stardata,
 
 #if (defined(RATES_OF_MARIA))
         /* Maria's REACTION Ne22pgNa23 -> SIGMAV_T22 */
-        sigmav[SIGMAV_T22]=nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.208200e+02,-0.108809e-03,-0.194097e+02,-0.372081e-01,0.739377e-02,-0.128009e-02,-0.648178)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.212579e+02,-0.410772,-0.580367e-02,0.189368e-01,-0.188510e-02,0.157242e-03,-0.150650e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.645203e+01,-0.175119e+01,-0.100753e-02,0.318931e-02,-0.309677e-03,0.253352e-04,-0.150111e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.134847e+02,-0.207387e+01,-0.109783e+01,0.137621e+02,-0.175262e+01,0.147358,-0.434513e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.101447e+02,-0.413567e+01,0.487085e+01,-0.130202e+02,0.973889,-0.663540e-01,0.368155e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.716413e+01,-0.479423e+01,-0.459213e+01,0.256030e+02,-0.219524e+01,0.130067,-0.771031e+01);
+        sigmav[SIGMAV_T22]=Reaclib(0.208200e+02,-0.108809e-03,-0.194097e+02,-0.372081e-01,0.739377e-02,-0.128009e-02,-0.648178)+Reaclib(-0.212579e+02,-0.410772,-0.580367e-02,0.189368e-01,-0.188510e-02,0.157242e-03,-0.150650e+01)+Reaclib(-0.645203e+01,-0.175119e+01,-0.100753e-02,0.318931e-02,-0.309677e-03,0.253352e-04,-0.150111e+01)+Reaclib(-0.134847e+02,-0.207387e+01,-0.109783e+01,0.137621e+02,-0.175262e+01,0.147358,-0.434513e+01)+Reaclib(0.101447e+02,-0.413567e+01,0.487085e+01,-0.130202e+02,0.973889,-0.663540e-01,0.368155e+01)+Reaclib(-0.716413e+01,-0.479423e+01,-0.459213e+01,0.256030e+02,-0.219524e+01,0.130067,-0.771031e+01);
 
 #elif (defined(RATES_OF_AMANDA))
         // Ne22(pg)Na23 NACRE rate
@@ -813,7 +821,7 @@ void nucsyn_set_sigmav(struct stardata_t * const stardata,
         sigmav[SIGMAV_T23]=(a1+a2+a3+a4+a5);
 #elif (defined(RATES_OF_MARIA))
         /* Maria's REACTION Na23paNe20 -> SIGMAV_T23 */
-        sigmav[SIGMAV_T23]=nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.156632e+03,-0.105827,0.326018e+01,-0.206279e+03,0.832450e+02,-0.499856e+02,0.403909e+02)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.525919e+03,-0.301515,0.601990e+02,-0.777241e+03,0.364372e+03,-0.315476e+03,0.147864e+03)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.119895e+03,-0.633817e-01,-0.581045e-01,0.192651,-0.190087e-01,0.154373e-02,-0.156567e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.360244e+02,-0.430617,0.911788e-02,-0.293146e-01,0.281655e-02,-0.224405e-03,-0.148985e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.776022e+01,-0.160524e+01,0.481135,-0.969025,0.395533e-01,-0.114858e-02,-0.101910e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.323173e+01,-0.199800e+01,0.328151e+01,-0.426201e+01,-0.945722e-01,0.254831e-01,0.141237e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.142176e+02,-0.323598e+01,0.589178e+01,-0.109880e+02,0.481526,-0.178061e-01,0.398717e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.132213e+02,-0.656750e+01,-0.132927e+01,0.406530e+01,-0.303075,0.182623e-01,-0.293149e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.154555e+01,-0.104906e+02,0.927879e+01,0.799835e+01,-0.143715e+01,0.114795,0.121158e+01);
+        sigmav[SIGMAV_T23]=Reaclib(0.156632e+03,-0.105827,0.326018e+01,-0.206279e+03,0.832450e+02,-0.499856e+02,0.403909e+02)+Reaclib(0.525919e+03,-0.301515,0.601990e+02,-0.777241e+03,0.364372e+03,-0.315476e+03,0.147864e+03)+Reaclib(-0.119895e+03,-0.633817e-01,-0.581045e-01,0.192651,-0.190087e-01,0.154373e-02,-0.156567e+01)+Reaclib(-0.360244e+02,-0.430617,0.911788e-02,-0.293146e-01,0.281655e-02,-0.224405e-03,-0.148985e+01)+Reaclib(-0.776022e+01,-0.160524e+01,0.481135,-0.969025,0.395533e-01,-0.114858e-02,-0.101910e+01)+Reaclib(0.323173e+01,-0.199800e+01,0.328151e+01,-0.426201e+01,-0.945722e-01,0.254831e-01,0.141237e+01)+Reaclib(0.142176e+02,-0.323598e+01,0.589178e+01,-0.109880e+02,0.481526,-0.178061e-01,0.398717e+01)+Reaclib(0.132213e+02,-0.656750e+01,-0.132927e+01,0.406530e+01,-0.303075,0.182623e-01,-0.293149e+01)+Reaclib(0.154555e+01,-0.104906e+02,0.927879e+01,0.799835e+01,-0.143715e+01,0.114795,0.121158e+01);
 
 #elif (defined(RATES_OF_AREND_JAN))
         sigmav[SIGMAV_T23]=(8.56e+09*TM23*exp(-20.766*TM13-Pow2(T9/0.131))*
@@ -863,7 +871,7 @@ void nucsyn_set_sigmav(struct stardata_t * const stardata,
             2.34e-4*t15*exp(-1.590*i9);
 #elif (defined(RATES_OF_MARIA))
         /* Maria's REACTION Na23pgMg24 -> SIGMAV_T23P */
-        sigmav[SIGMAV_T23P]=nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.321402e+02,-0.226477e-01,-0.169774e+02,-0.188483e+02,0.333428e+01,-0.862626,0.444836e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.122215e+03,-0.638582e-01,0.455271e-02,-0.174054e-01,0.192201e-02,-0.168311e-03,-0.149447e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.460634e+02,-0.430413,-0.171111e-01,0.560886e-01,-0.546019e-02,0.438306e-03,-0.151924e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.694817e+01,-0.160160e+01,-0.104266,0.129067e+01,-0.164065,0.133457e-01,-0.174908e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.444413e+01,-0.278984e+01,-0.636430e-03,0.196576e-02,-0.187741e-03,0.149812e-04,-0.150069e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.978039e+01,-0.343464e+01,0.613711e-01,-0.958165e-01,0.487948e-02,-0.299863e-03,-0.144865e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.889387e+01,-0.570055e+01,-0.367651e+01,0.254056e+02,-0.205687e+01,0.120096,-0.795560e+01);
+        sigmav[SIGMAV_T23P]=Reaclib(0.321402e+02,-0.226477e-01,-0.169774e+02,-0.188483e+02,0.333428e+01,-0.862626,0.444836e+01)+Reaclib(-0.122215e+03,-0.638582e-01,0.455271e-02,-0.174054e-01,0.192201e-02,-0.168311e-03,-0.149447e+01)+Reaclib(-0.460634e+02,-0.430413,-0.171111e-01,0.560886e-01,-0.546019e-02,0.438306e-03,-0.151924e+01)+Reaclib(-0.694817e+01,-0.160160e+01,-0.104266,0.129067e+01,-0.164065,0.133457e-01,-0.174908e+01)+Reaclib(0.444413e+01,-0.278984e+01,-0.636430e-03,0.196576e-02,-0.187741e-03,0.149812e-04,-0.150069e+01)+Reaclib(0.978039e+01,-0.343464e+01,0.613711e-01,-0.958165e-01,0.487948e-02,-0.299863e-03,-0.144865e+01)+Reaclib(-0.889387e+01,-0.570055e+01,-0.367651e+01,0.254056e+02,-0.205687e+01,0.120096,-0.795560e+01);
 
 #elif (defined(RATES_OF_AREND_JAN))
         sigmav[SIGMAV_T23P]=(2.93e+08*TM23*exp(-20.776*TM13-Pow2(T9/0.297))
@@ -938,7 +946,7 @@ void nucsyn_set_sigmav(struct stardata_t * const stardata,
 #endif// NUCSYN_HOT_SIGMAV
 #elif (defined(RATES_OF_MARIA))
         /* Maria's REACTION Mg24pgAl25 -> SIGMAV_T24 */
-        sigmav[SIGMAV_T24]=nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.199565e+02,-0.697781e-04,-0.219525e+02,-0.162580e-01,0.366352e-02,-0.321486e-03,-0.655153)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.767211e+01,-0.248344e+01,0.204845e-01,-0.574066e-01,0.527151e-02,-0.420928e-03,-0.147899e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.178763e+02,-0.441813e+01,-0.215198e+02,0.509589e+02,-0.289045e+01,0.132303,-0.220245e+02);
+        sigmav[SIGMAV_T24]=Reaclib(0.199565e+02,-0.697781e-04,-0.219525e+02,-0.162580e-01,0.366352e-02,-0.321486e-03,-0.655153)+Reaclib(0.767211e+01,-0.248344e+01,0.204845e-01,-0.574066e-01,0.527151e-02,-0.420928e-03,-0.147899e+01)+Reaclib(-0.178763e+02,-0.441813e+01,-0.215198e+02,0.509589e+02,-0.289045e+01,0.132303,-0.220245e+02);
 
 #elif (defined(RATES_OF_AREND_JAN))
         sigmav[SIGMAV_T24]=(5.60e+08/T923*exp(-22.019/T913)
@@ -1035,7 +1043,7 @@ void nucsyn_set_sigmav(struct stardata_t * const stardata,
 
 #elif (defined(RATES_OF_MARIA))
         /* Maria's REACTION Mg25pgAl26 -> SIGMAV_T25 */
-        sigmav[SIGMAV_T25]=nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.169033e+02,-0.666234,-0.130374e-02,0.455645e-02,-0.478223e-03,0.412771e-04,-0.150151e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.940794e+01,-0.107258e+01,0.682242,-0.212774e+01,0.167209,-0.116951e-01,-0.692396)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.336968e+02,-0.258601e+01,0.367150e+02,-0.727831e+02,0.395660e+01,-0.196978,0.322550e+02)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.197844e+02,-0.348095e+01,-0.847000e+01,0.416773e+02,-0.358991e+01,0.238561,-0.132509e+02)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.171000e+02,-0.666234,-0.130471e-02,0.456896e-02,-0.480198e-03,0.414636e-04,-0.150151e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.979089e+01,-0.107201e+01,0.589102,-0.183108e+01,0.143440,-0.100112e-01,-0.803649)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.381956e+02,-0.261289e+01,0.397609e+02,-0.811782e+02,0.462963e+01,-0.244607,0.354598e+02)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.840186e+01,-0.352648e+01,0.304344e-01,-0.852458e-01,0.782572e-02,-0.624801e-03,-0.146880e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.522874e+01,-0.440515e+01,0.408244e+01,0.115893e+01,0.650368e-01,-0.321815e-01,0.677192)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.185678e+02,-0.666233,-0.131744e-02,0.460061e-02,-0.482097e-03,0.415574e-04,-0.150153e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.100318e+02,-0.107603e+01,0.125711e+01,-0.399984e+01,0.319997,-0.226480e-01,0.125236e-02)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.212179e+02,-0.251613e+01,0.291324e+02,-0.528962e+02,0.242770e+01,-0.919813e-01,0.244690e+02)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.650467e+01,-0.352648e+01,0.303677e-01,-0.850962e-01,0.781329e-02,-0.623843e-03,-0.146886e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.225344e+02,-0.425986e+01,-0.115838e+02,0.471413e+02,-0.391512e+01,0.256920,-0.159977e+02);
+        sigmav[SIGMAV_T25]=Reaclib(-0.169033e+02,-0.666234,-0.130374e-02,0.455645e-02,-0.478223e-03,0.412771e-04,-0.150151e+01)+Reaclib(-0.940794e+01,-0.107258e+01,0.682242,-0.212774e+01,0.167209,-0.116951e-01,-0.692396)+Reaclib(0.336968e+02,-0.258601e+01,0.367150e+02,-0.727831e+02,0.395660e+01,-0.196978,0.322550e+02)+Reaclib(-0.197844e+02,-0.348095e+01,-0.847000e+01,0.416773e+02,-0.358991e+01,0.238561,-0.132509e+02)+Reaclib(-0.171000e+02,-0.666234,-0.130471e-02,0.456896e-02,-0.480198e-03,0.414636e-04,-0.150151e+01)+Reaclib(-0.979089e+01,-0.107201e+01,0.589102,-0.183108e+01,0.143440,-0.100112e-01,-0.803649)+Reaclib(0.381956e+02,-0.261289e+01,0.397609e+02,-0.811782e+02,0.462963e+01,-0.244607,0.354598e+02)+Reaclib(0.840186e+01,-0.352648e+01,0.304344e-01,-0.852458e-01,0.782572e-02,-0.624801e-03,-0.146880e+01)+Reaclib(0.522874e+01,-0.440515e+01,0.408244e+01,0.115893e+01,0.650368e-01,-0.321815e-01,0.677192)+Reaclib(-0.185678e+02,-0.666233,-0.131744e-02,0.460061e-02,-0.482097e-03,0.415574e-04,-0.150153e+01)+Reaclib(-0.100318e+02,-0.107603e+01,0.125711e+01,-0.399984e+01,0.319997,-0.226480e-01,0.125236e-02)+Reaclib(0.212179e+02,-0.251613e+01,0.291324e+02,-0.528962e+02,0.242770e+01,-0.919813e-01,0.244690e+02)+Reaclib(0.650467e+01,-0.352648e+01,0.303677e-01,-0.850962e-01,0.781329e-02,-0.623843e-03,-0.146886e+01)+Reaclib(-0.225344e+02,-0.425986e+01,-0.115838e+02,0.471413e+02,-0.391512e+01,0.256920,-0.159977e+02);
 
 #elif (defined(RATES_OF_AREND_JAN))
         sigmav[SIGMAV_T25]=(3.57e+09/T923*exp(-22.031/T913-Pow2(T9/0.06))
@@ -1133,7 +1141,7 @@ void nucsyn_set_sigmav(struct stardata_t * const stardata,
 #endif
 #elif (defined(RATES_OF_MARIA))
         /* Maria's REACTION Al26pgSi27 -> SIGMAV_T26P */
-        sigmav[SIGMAV_T26P]=nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.211853e+02,-0.181149e-03,-0.231640e+02,-0.806922e-01,0.126339e-01,-0.218412e-02,-0.635626)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.194622e+02,-0.799519,-0.833775e-02,0.263742e-01,-0.254759e-02,0.207772e-03,-0.150919e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.140023e+02,-0.108957e+01,-0.192068e-01,0.655138e-01,-0.672448e-02,0.570960e-03,-0.152200e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.915706e+01,-0.149593e+01,0.568974e-02,-0.184887e-01,0.183024e-02,-0.151859e-03,-0.149364e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.217872e+01,-0.218510e+01,-0.111226e-01,0.274505e-01,-0.223807e-02,0.163224e-03,-0.151072e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.615109e+02,-0.285261e+01,0.251283e+02,-0.845335e+02,0.641356e+01,-0.435046,0.306598e+02)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.185333e+02,-0.820942e+01,0.122707e+02,-0.209107e+02,0.105220e+01,-0.540533e-01,0.905230e+01);
+        sigmav[SIGMAV_T26P]=Reaclib(0.211853e+02,-0.181149e-03,-0.231640e+02,-0.806922e-01,0.126339e-01,-0.218412e-02,-0.635626)+Reaclib(-0.194622e+02,-0.799519,-0.833775e-02,0.263742e-01,-0.254759e-02,0.207772e-03,-0.150919e+01)+Reaclib(-0.140023e+02,-0.108957e+01,-0.192068e-01,0.655138e-01,-0.672448e-02,0.570960e-03,-0.152200e+01)+Reaclib(-0.915706e+01,-0.149593e+01,0.568974e-02,-0.184887e-01,0.183024e-02,-0.151859e-03,-0.149364e+01)+Reaclib(0.217872e+01,-0.218510e+01,-0.111226e-01,0.274505e-01,-0.223807e-02,0.163224e-03,-0.151072e+01)+Reaclib(0.615109e+02,-0.285261e+01,0.251283e+02,-0.845335e+02,0.641356e+01,-0.435046,0.306598e+02)+Reaclib(0.185333e+02,-0.820942e+01,0.122707e+02,-0.209107e+02,0.105220e+01,-0.540533e-01,0.905230e+01);
 
 #elif (defined(RATES_OF_AREND_JAN))
         sigmav[SIGMAV_T26P]=(6.78e+13/T923*exp((-23.261/T913)
@@ -1263,7 +1271,7 @@ void nucsyn_set_sigmav(struct stardata_t * const stardata,
 #elif (defined(RATES_OF_MARIA))
         /* Maria's REACTION Mg26pgAl27 -> SIGMAV_T26 */
 
-        sigmav[SIGMAV_T26]=nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.367326e+02,-0.285822e-01,-0.171860e+02,-0.241996e+02,0.454763e+01,-0.134721e+01,0.584008e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.264345e+02,-0.610945,-0.387739e-01,0.130991,-0.134051e-01,0.113848e-02,-0.154416e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.557003e+01,-0.106206e+01,0.564415e+01,-0.242421e+02,0.228001e+01,-0.177490,0.653801e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.266603e+02,-0.124876e+01,-0.277958e+02,0.532718e+02,-0.384878e+01,0.264868,-0.245900e+02)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.224299e+01,-0.254252e+01,0.101529,-0.263173,0.221733e-01,-0.168360e-02,-0.139841e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.799631e+01,-0.305480e+01,-0.152255e+02,0.322869e+02,-0.244178e+01,0.171841,-0.149407e+02)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.106582e+02,-0.378911e+01,0.303665,-0.780598,0.127192e+01,-0.144673,-0.136316e+01);
+        sigmav[SIGMAV_T26]=Reaclib(0.367326e+02,-0.285822e-01,-0.171860e+02,-0.241996e+02,0.454763e+01,-0.134721e+01,0.584008e+01)+Reaclib(-0.264345e+02,-0.610945,-0.387739e-01,0.130991,-0.134051e-01,0.113848e-02,-0.154416e+01)+Reaclib(0.557003e+01,-0.106206e+01,0.564415e+01,-0.242421e+02,0.228001e+01,-0.177490,0.653801e+01)+Reaclib(-0.266603e+02,-0.124876e+01,-0.277958e+02,0.532718e+02,-0.384878e+01,0.264868,-0.245900e+02)+Reaclib(0.224299e+01,-0.254252e+01,0.101529,-0.263173,0.221733e-01,-0.168360e-02,-0.139841e+01)+Reaclib(-0.799631e+01,-0.305480e+01,-0.152255e+02,0.322869e+02,-0.244178e+01,0.171841,-0.149407e+02)+Reaclib(0.106582e+02,-0.378911e+01,0.303665,-0.780598,0.127192e+01,-0.144673,-0.136316e+01);
 
 
 #elif (defined(RATES_OF_AREND_JAN))
@@ -1431,7 +1439,7 @@ void nucsyn_set_sigmav(struct stardata_t * const stardata,
 
 #if (defined(RATES_OF_MARIA))
         /* Maria's REACTION Al27pgSi28 -> SIGMAV_T27 */
-        sigmav[SIGMAV_T27]=nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.365503e+02,-0.275874e-01,-0.185636e+02,-0.234291e+02,0.429262e+01,-0.117105e+01,0.562347e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.229525e+02,-0.832147,0.892092e-02,-0.292581e-01,0.291664e-02,-0.243224e-03,-0.148998e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.188665e+02,-0.980667,0.539901e-02,-0.161642e-01,0.149001e-02,-0.117713e-03,-0.149421e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.173054e+02,-0.233495e+01,0.972645e+01,-0.248307e+02,0.148285e+01,-0.831689e-01,0.932530e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,0.901190e+01,-0.379595e+01,0.115778e+02,-0.137693e+02,-0.876306e-01,0.573784e-01,0.762694e+01)+nucsyn_reaclib(t9,i9,tm13,t913,t953,lt9,-0.258564e+02,-0.549543e+01,-0.180073e+02,0.585367e+02,-0.443583e+01,0.278242,-0.214648e+02);
+        sigmav[SIGMAV_T27]=Reaclib(0.365503e+02,-0.275874e-01,-0.185636e+02,-0.234291e+02,0.429262e+01,-0.117105e+01,0.562347e+01)+Reaclib(-0.229525e+02,-0.832147,0.892092e-02,-0.292581e-01,0.291664e-02,-0.243224e-03,-0.148998e+01)+Reaclib(-0.188665e+02,-0.980667,0.539901e-02,-0.161642e-01,0.149001e-02,-0.117713e-03,-0.149421e+01)+Reaclib(0.173054e+02,-0.233495e+01,0.972645e+01,-0.248307e+02,0.148285e+01,-0.831689e-01,0.932530e+01)+Reaclib(0.901190e+01,-0.379595e+01,0.115778e+02,-0.137693e+02,-0.876306e-01,0.573784e-01,0.762694e+01)+Reaclib(-0.258564e+02,-0.549543e+01,-0.180073e+02,0.585367e+02,-0.443583e+01,0.278242,-0.214648e+02);
 
 
 #elif (defined(RATES_OF_AREND_JAN))
@@ -1602,8 +1610,39 @@ void nucsyn_set_sigmav(struct stardata_t * const stardata,
             sigmav_Be8 *
             part3;
 
+        /*
+         * reaclib
+         */
+
+        sigmav[SIGMAV_3He4] =
+            Reaclib(-9.710520e-01,	0.000000e+00,	-3.706000e+01,	2.934930e+01,	-1.155070e+02,	-1.000000e+01,	-1.333330e+00) +
+            Reaclib(-1.178840e+01,	-1.024460e+00,	-2.357000e+01,	2.048860e+01,	-1.298820e+01,	-2.000000e+01,	-2.166670e+00) +
+            Reaclib(-2.435050e+01,	-4.126560e+00,	-1.349000e+01,	2.142590e+01,	-1.347690e+00,	8.798160e-02,	-1.316530e+01);
+
         sigmav[SIGMAV_3He4] /= Pow2(AVOGADRO_CONSTANT);
 
+        /*
+         * Reaction rate:
+         *
+         * Two particle
+         * r = lambda * N = <sigmav> N1 N2
+         *
+         * Three particle
+         * r = <sigmav> N1 N2 N3
+         *
+         * hence, if N is the number of helium, which
+         * is the rate of C12 production,
+         *
+         * r = 1/6 * NA^2 * rho^2 * Y^3 <sigmav> N*N*N
+         *
+         * and the rate of He4 destruction is 3* this.
+         *
+         * lambda = ln2/tau
+         *
+         *
+         */
+
+
         /*
          * C12 + He4 -> O16 + gamma
          * page 67
diff --git a/src/perl/scripts2/make_population_ensemble.pl b/src/perl/scripts2/make_population_ensemble.pl
index 895fc82d9c0d382c4fb7a87c35a1bcb47452bd05..5829df6bd25bf23e69c5419ce0ec9b1cd075ece9 100755
--- a/src/perl/scripts2/make_population_ensemble.pl
+++ b/src/perl/scripts2/make_population_ensemble.pl
@@ -1249,8 +1249,8 @@ sub local_setup
             # mixed single/binary : weighting is binary fraction
             'my $binary_fraction = main::binary_fraction($self,$m1);  $self->{_grid_options}{weight}= $self->{_grid_options}{multiplicity}==1 ? 1.0-$binary_fraction : $binary_fraction;   ' :
 
-                                                                                                                                                                                         # single or binary stars : weighting is just 1.0
-                                                                                                                                                                                         '  my $binary_fraction=($self->{_grid_options}{multiplicity}==1); $self->{_grid_options}{weight}=1.0; ';
+                                                                                                                            # single or binary stars : weighting is just 1.0
+           '$self->{_grid_options}{weight}=1.0; ';
 
         # Mass 1
         if($distribution_options->{m1spacing} eq 'log')
diff --git a/src/setup/cmd_line_args_list.h b/src/setup/cmd_line_args_list.h
index da8173423e5a61b38d8110c853542eeb1dfefaee..9260dfac954d945e6c8084a4fccaae86912cb6db 100644
--- a/src/setup/cmd_line_args_list.h
+++ b/src/setup/cmd_line_args_list.h
@@ -3606,6 +3606,16 @@ ARG_SECTION_INPUT,
         "",
         MINT_KIPPENHAHN_VAR,
         1.0
+    CMD_LINE_ARG_T_OPTIONAL_ARGS },
+
+{
+    ARG_SECTION_INPUT,
+        "MINT_nshells",
+        "Set the number of shells MINT uses in each star when doing nuclear burning. (200)",
+        ARG_INTEGER ,
+        "",
+        MINT_NSHELLS_VAR,
+        1.0
         CMD_LINE_ARG_T_OPTIONAL_ARGS },
 
 {
diff --git a/src/setup/cmd_line_macros.h b/src/setup/cmd_line_macros.h
index 5a5cc32d0c1a109f11935c0e53239642c84b5018..b5c8536ed1cee70a055d5b9141027444ac8ed124 100644
--- a/src/setup/cmd_line_macros.h
+++ b/src/setup/cmd_line_macros.h
@@ -36,6 +36,7 @@
 #define MINT_MINIMUM_SHELL_MASS_VAR Var(stardata->preferences->MINT_minimum_shell_mass)
 #define MINT_MAXIMUM_SHELL_MASS_VAR Var(stardata->preferences->MINT_maximum_shell_mass)
 #define MINT_FALLBACK_TO_TEST_DATA_VAR Var(stardata->preferences->MINT_fallback_to_test_data)
+#define MINT_NSHELLS_VAR Var(stardata->preferences->MINT_nshells)
 #else
 #define MINT_DATA_CLEANUP_VAR Not_used_arg
 #define MINT_METALLICITY_VAR Not_used_arg
@@ -51,6 +52,7 @@
 #define MINT_MINIMUM_SHELL_MASS_VAR Not_used_arg
 #define MINT_MAXIMUM_SHELL_MASS_VAR Not_used_arg
 #define MINT_FALLBACK_TO_TEST_DATA_VAR Not_used_arg
+#define MINT_NSHELLS_VAR Not_used_arg
 #endif//MINT
 
 #ifdef COMENV_MS_ACCRETION
diff --git a/src/setup/set_default_preferences.c b/src/setup/set_default_preferences.c
index 85d4c43b66ee98dadec6d6d09ec093bb2ea7e473..66326c959fc84cd0c5f4215371f1a12d5d83f038 100644
--- a/src/setup/set_default_preferences.c
+++ b/src/setup/set_default_preferences.c
@@ -567,6 +567,7 @@ preferences->WD_accretion_rate_novae_upper_limit_other_donor = ACCRETION_RATE_NO
     preferences->MINT_disable_grid_load_warnings = FALSE;
     preferences->MINT_use_ZAMS_profiles = TRUE;
     preferences->MINT_fallback_to_test_data = FALSE;
+    preferences->MINT_nshells = 0;
 #endif // MINT
 
     /* multiplicity is unset */
diff --git a/src/stellar_structure/stellar_structure_BSE.c b/src/stellar_structure/stellar_structure_BSE.c
index 8a8bb9caa12c71737ec3533fa7846fcfd9b2ff50..7227d0a419f2b23a75ff3b73746d0655390b2c4b 100644
--- a/src/stellar_structure/stellar_structure_BSE.c
+++ b/src/stellar_structure/stellar_structure_BSE.c
@@ -149,6 +149,7 @@ int stellar_structure_BSE_with_newstar(struct stardata_t * const stardata,
            Yesno(alloc_luminosities),
            Yesno(alloc_timescales),
            Yesno(alloc_GB));
+
     stellar_structure_BSE_given_timescales(newstar,
                                            oldstar,
                                            stardata,
diff --git a/src/stellar_structure/stellar_structure_CHeB.c b/src/stellar_structure/stellar_structure_CHeB.c
index f6c1c7aba162c6c2ac7bf700bbbe828c7ca842a4..7d682c02eda0d058e019ae4a79a3abc6af53bcf1 100644
--- a/src/stellar_structure/stellar_structure_CHeB.c
+++ b/src/stellar_structure/stellar_structure_CHeB.c
@@ -45,6 +45,17 @@ static Stellar_type stellar_structure_CHeB_implementation(struct star_t *newstar
     /*
      * Core helium burning star.
      */
+#ifdef MINT
+    if(stardata->preferences->stellar_structure_algorithm ==
+       STELLAR_STRUCTURE_ALGORITHM_MINT)
+    {
+        /* let MINT handle this phase */
+        return MINT_stellar_structure_CHeB(stardata,
+                                           NULL,
+                                           newstar);
+    }
+#endif//MINT
+
     const double * metallicity_parameters=stardata->common.metallicity_parameters;
     const double * giant_branch_parameters=stardata->common.giant_branch_parameters;
     double tau,lx,ly,rx,ry,rmin,texp,taul,tauh,tau2,tloop;
diff --git a/src/stellar_timescales/stellar_timescales_high_mass_GB.c b/src/stellar_timescales/stellar_timescales_high_mass_GB.c
index 4403819015d2d71ea03d729f18ac8ee2d2778ec4..23b6eee8130e653a225a6a2edbfc26bf9ee5fe34 100644
--- a/src/stellar_timescales/stellar_timescales_high_mass_GB.c
+++ b/src/stellar_timescales/stellar_timescales_high_mass_GB.c
@@ -2,14 +2,15 @@
 No_empty_translation_unit_warning;
 
 #ifdef BSE
-void stellar_timescales_high_mass_GB(const double mass,
+void stellar_timescales_high_mass_GB(struct stardata_t * const stardata,
+                                     const double mass,
                                      const double mp2,
                                      double * const timescales,
                                      double * const luminosities,
                                      const double * RESTRICT const giant_branch_parameters)
 {
 #ifndef USE_BSE_TIMESCALES_H
-    /* 
+    /*
      * Note that for M>metallicity_parameters(3) there is no
      * GB as the star goes from
      * HG -> CHeB -> AGB. So in effect timescales(1) refers to the time of
@@ -21,14 +22,35 @@ void stellar_timescales_high_mass_GB(const double mass,
                                     mp2,
                                     giant_branch_parameters)
         *timescales[T_BGB];
-    
+
 #endif
-    /* 
-     * This now represents the luminosity at the end of CHeB, ie. BAGB
-     */
-    luminosities[L_HE_BURNING] = luminosities[L_BAGB];
 
-    /* 
+    Boolean have_set = FALSE;
+#ifdef MINT
+    if(stardata->preferences->stellar_structure_algorithm ==
+       STELLAR_STRUCTURE_ALGORITHM_MINT)
+    {
+        have_set = TRUE;
+
+        MINT_CHeB_ZA(stardata,
+                     mass,
+                     &luminosities[L_HE_BURNING]);
+        MINT_CHeB_TA(stardata,
+                     mass,
+                     &timescales[T_HE_BURNING],
+                     NULL);
+    }
+#endif
+
+    if(have_set == FALSE)
+    {
+        /*
+         * This now represents the luminosity at the end of CHeB, ie. BAGB
+         */
+        luminosities[L_HE_BURNING] = luminosities[L_BAGB];
+    }
+
+    /*
      * We set luminosities(3) to be the luminosity at the end of the HG
      */
     luminosities[L_BGB] = luminosities[L_HE_IGNITION];
diff --git a/src/stellar_timescales/stellar_timescales_hydrogen_stars.c b/src/stellar_timescales/stellar_timescales_hydrogen_stars.c
index 6e68a95ddea50246c816f9358d7297e6afad709c..7399203779d347c409c417eece3dd4ec4d368b22 100644
--- a/src/stellar_timescales/stellar_timescales_hydrogen_stars.c
+++ b/src/stellar_timescales/stellar_timescales_hydrogen_stars.c
@@ -283,8 +283,6 @@ int stellar_timescales_hydrogen_stars(const double phase_start_mass,
 
         Dprint("compare phase_start_mass=%g vs metallicity_parameter[3]=%g\n",phase_start_mass,metallicity_parameters[ZPAR_MASS_FGB]);
 
-
-
         if(Less_or_equal(phase_start_mass,
                          metallicity_parameters[ZPAR_MASS_FGB]))
         {
@@ -302,12 +300,14 @@ int stellar_timescales_hydrogen_stars(const double phase_start_mass,
         }
         else
         {
-            stellar_timescales_high_mass_GB(phase_start_mass,
+            stellar_timescales_high_mass_GB(stardata,
+                                            phase_start_mass,
                                             metallicity_parameters[ZPAR_MASS_HE_FLASH],
                                             timescales,
                                             luminosities,
                                             giant_branch_parameters);
         }
+
 #endif // USE_BSE_TIMESCALES_H
 
         Dprint("Done compare stellar_type %d\n",stellar_type);
diff --git a/src/stellar_timescales/stellar_timescales_low_mass_GB.c b/src/stellar_timescales/stellar_timescales_low_mass_GB.c
index da607eeba0836092ac6cf51ba1e87a812d5aa624..d37b5ab5cc4f881b787215c25bb020a7f2035ac5 100644
--- a/src/stellar_timescales/stellar_timescales_low_mass_GB.c
+++ b/src/stellar_timescales/stellar_timescales_low_mass_GB.c
@@ -92,51 +92,70 @@ void stellar_timescales_low_mass_GB(const double GBp,
            timescales[T_HE_IGNITION],
            luminosities[L_HE_IGNITION]);
 
-    if(Less_or_equal(mass,metallicity_parameters[ZPAR_MASS_HE_FLASH]))
+
+
+    Boolean have_set = FALSE;
+#ifdef MINT
+    if(stardata->preferences->stellar_structure_algorithm ==
+       STELLAR_STRUCTURE_ALGORITHM_MINT)
     {
-        double mc1 =
-            mcgbf(luminosities[L_HE_IGNITION],
-                  GB,
-                  luminosities[L_LMX]);
-
-
-        luminosities[L_HE_BURNING] =
-            lzahbf(mass,
-                   mc1,
-                   metallicity_parameters[ZPAR_MASS_HE_FLASH],
-                   giant_branch_parameters);
-
-        timescales[T_HE_BURNING] =
-            thef(mass,
-                 mc1,
-                 metallicity_parameters[ZPAR_MASS_HE_FLASH],
-                 giant_branch_parameters);
-
-        Dprint("M<=Mflash : L=%30.20e timescale=%30.20e\n",
-               luminosities[L_HE_BURNING],
-               timescales[T_HE_BURNING]);
+        have_set = TRUE;
+        MINT_CHeB_ZA(stardata,
+                     mass,
+                     &luminosities[L_HE_BURNING]);
+        MINT_CHeB_TA(stardata,
+                     mass,
+                     &timescales[T_HE_BURNING],
+                     NULL);
     }
-    else
-    {
-        luminosities[L_HE_BURNING] =
-            lhef(mass,
-                 giant_branch_parameters)
-            *luminosities[L_HE_IGNITION];
-
-        timescales[T_HE_BURNING] =
-            thef(mass,
-                 1.0,
-                 metallicity_parameters[ZPAR_MASS_HE_FLASH],
-                 giant_branch_parameters)
-            *timescales[T_BGB];
-
-        Dprint("M>Mflash : L=%30.20e timescale=%30.20e\n",
-               luminosities[L_HE_BURNING],
-               timescales[T_HE_BURNING]);
+#endif
 
+    if(have_set == FALSE)
+    {
+        if(Less_or_equal(mass,metallicity_parameters[ZPAR_MASS_HE_FLASH]))
+        {
+            const double mc1 =
+                mcgbf(luminosities[L_HE_IGNITION],
+                      GB,
+                      luminosities[L_LMX]);
+
+            luminosities[L_HE_BURNING] =
+                lzahbf(mass,
+                       mc1,
+                       metallicity_parameters[ZPAR_MASS_HE_FLASH],
+                       giant_branch_parameters);
+
+            timescales[T_HE_BURNING] =
+                thef(mass,
+                     mc1,
+                     metallicity_parameters[ZPAR_MASS_HE_FLASH],
+                     giant_branch_parameters);
+
+            Dprint("M<=Mflash : L=%30.20e timescale=%30.20e\n",
+                   luminosities[L_HE_BURNING],
+                   timescales[T_HE_BURNING]);
+        }
+        else
+        {
+            luminosities[L_HE_BURNING] =
+                lhef(mass,
+                     giant_branch_parameters)
+                *luminosities[L_HE_IGNITION];
+
+            timescales[T_HE_BURNING] =
+                thef(mass,
+                     1.0,
+                     metallicity_parameters[ZPAR_MASS_HE_FLASH],
+                     giant_branch_parameters)
+                *timescales[T_BGB];
+
+            Dprint("M>Mflash : L=%30.20e timescale=%30.20e\n",
+                   luminosities[L_HE_BURNING],
+                   timescales[T_HE_BURNING]);
+
+        }
     }
 
 
-
 }
 #endif//BSE
diff --git a/src/stellar_timescales/stellar_timescales_prototypes.h b/src/stellar_timescales/stellar_timescales_prototypes.h
index 99824f0d8e4cd08ada06f9da9c011a953ee16c05..5e0de4bf43b0332a4e650c9183f28305f4b59f7d 100644
--- a/src/stellar_timescales/stellar_timescales_prototypes.h
+++ b/src/stellar_timescales/stellar_timescales_prototypes.h
@@ -68,7 +68,8 @@ int stellar_timescales_post_HG(const double tbagb,
                                struct stardata_t * const stardata,
                                struct star_t * const star,
                                const Stellar_type stellar_type);
-void stellar_timescales_high_mass_GB(const double mass,
+void stellar_timescales_high_mass_GB(struct stardata_t * const stardata,
+                                     const double mass,
                                      const double mp2,
                                      double * const timescales,
                                      double * const luminosities,
diff --git a/src/binary_star_functions/hut_polynomials.c b/src/tides/hut_polynomials.c
similarity index 100%
rename from src/binary_star_functions/hut_polynomials.c
rename to src/tides/hut_polynomials.c
diff --git a/src/tides/tides_prototypes.h b/src/tides/tides_prototypes.h
index c4a97f41cd8c76602c8eeb13706965515cb01730..41dcf193babd403e41d5764a3ca032b4a087ea76 100644
--- a/src/tides/tides_prototypes.h
+++ b/src/tides/tides_prototypes.h
@@ -28,4 +28,15 @@ void tides_generic(struct stardata_t * const stardata,
                    double * const dedt_star,
                    double * const dedt_orbit);
 
+Pure_function
+double Hut_f5(const double ecc);
+Pure_function
+double Hut_f4(const double ecc);
+Pure_function
+double Hut_f3(const double ecc);
+Pure_function
+double Hut_f2(const double ecc);
+Pure_function
+double Hut_f1(const double ecc);
+
 #endif // TIDES_PROTOTYPES_H
diff --git a/src/timestep/timestep_CHeB.c b/src/timestep/timestep_CHeB.c
index f536dad0cdb6ea63a92cd24b37f91b6d4e510944..087722251ffeb21871ceeea7b90a4a480840bd2a 100644
--- a/src/timestep/timestep_CHeB.c
+++ b/src/timestep/timestep_CHeB.c
@@ -5,24 +5,38 @@ No_empty_translation_unit_warning;
 
 void timestep_CHeB(Timestep_prototype_args)
 {
-#ifdef BSE
-    if(Use_timestep(T_HE_BURNING))
+    if(stardata->preferences->stellar_structure_algorithm == STELLAR_STRUCTURE_ALGORITHM_MINT)
     {
-        Set_timestep(*dt,stardata->preferences->timestep_multipliers[DT_LIMIT_CHeB]*tscls[T_HE_BURNING],star,DT_LIMIT_CHeB);
+#ifdef MINT
+        double phase_lifetime;
+        MINT_CHeB_TA(stardata,
+                     star->phase_start_mass,
+                     &phase_lifetime,
+                     NULL);
+
+        Set_timestep(*dt,stardata->preferences->timestep_multipliers[DT_LIMIT_CHeB]*phase_lifetime,star,DT_LIMIT_CHeB);
+#endif//MINT
     }
+    else
+    {
+#ifdef BSE
+        if(Use_timestep(T_HE_BURNING))
+        {
+            Set_timestep(*dt,stardata->preferences->timestep_multipliers[DT_LIMIT_CHeB]*tscls[T_HE_BURNING],star,DT_LIMIT_CHeB);
+        }
 #ifdef HRDIAG
-    stardata->model.hrdiag_dt_guess = HRDIAG_FRAC* tscls[T_HE_BURNING];
+        stardata->model.hrdiag_dt_guess = HRDIAG_FRAC* tscls[T_HE_BURNING];
 #ifdef HRDIAG_DTLOG
-    printf("HRDT (deltat) CHeB %g\n",stardata->model.hrdiag_dt_guess);
+        printf("HRDT (deltat) CHeB %g\n",stardata->model.hrdiag_dt_guess);
 #endif // HRDIAG_DTLOG
-    HRDIAG_DT_NEGATIVE_CHECK;
+        HRDIAG_DT_NEGATIVE_CHECK;
 #endif // HRDIAG
+    }
 
- 
     double accretion_rate = Mdot_net(star);
     if(Is_zero(accretion_rate))
     {
-        *time_remaining = 
+        *time_remaining =
             Min(tn,
                 (tscls[T_HE_BURNING] + tscls[T_HE_IGNITION]))
             - age;
@@ -38,8 +52,8 @@ void timestep_CHeB(Timestep_prototype_args)
         /*
          * Accretion CHeB star has an unknown lifetime,
          * because it depends on the accretion timescale,
-         * which is generally unknown. Assume the mass  
-         * transfer rate timescale limiter will resolve 
+         * which is generally unknown. Assume the mass
+         * transfer rate timescale limiter will resolve
          * the timescale.
          */
         *time_remaining = LONG_TIME_REMAINING;
diff --git a/tbse b/tbse
index 7472c373f70474ffeccb7cee89dff33673b7d844..ee67f9ea2b7716b23f3ae916e8a7cbfd898a883a 100755
--- a/tbse
+++ b/tbse
@@ -1948,6 +1948,9 @@ MINT_USE_ZAMS_PROFILES=True
 # turn shellular nuclear burning on or off (False)
 MINT_NUCLEAR_BURNING=True
 
+# number of shells (200)
+MINT_NSHELLS=200
+
 # turn on MINT remeshing (True)
 MINT_REMESH=True
 
@@ -2572,6 +2575,7 @@ $TIMESTEP_MULTIPLIERS \
 --MINT_disable_grid_load_warnings $MINT_DISABLE_GRID_LOAD_WARNINGS \
 --MINT_MS_rejuvenation $MINT_MS_REJUVENATION \
 --MINT_nuclear_burning $MINT_NUCLEAR_BURNING \
+--MINT_nshells $MINT_NSHELLS \
 --MINT_minimum_shell_mass $MINT_MINIMUM_SHELL_MASS \
 --MINT_maximum_shell_mass $MINT_MAXIMUM_SHELL_MASS \
 --MINT_remesh $MINT_REMESH \
diff --git a/tbse.mint b/tbse.mint
index 1ecbd41b345c886bb4aa7bedcc5fd1babc169226..5b5eb409756c8c694679747c08130d1e0feb8504 100755
--- a/tbse.mint
+++ b/tbse.mint
@@ -18,12 +18,13 @@ MAGNETIC_BRAKING_GAMMA=3.0
 
 # MINT options
 MINT_METALLICITY=0.02 # set to -1 to ignore
-MINT_KIPPENHAHN=0
-MINT_KIPPENHAHN_STELLAR_TYPE=-1
-MINT_KIPPENHAHN_COMPANION_STELLAR_TYPE=-1
-MINT_MINIMUM_SHELL_MASS=1e-3
-MINT_MAXIMUM_SHELL_MASS=1e-2
-MINT_NUCLEAR_BURNING=0
+MINT_KIPPENHAHN=1 # 0=off, 1=star1, 2=star2
+MINT_KIPPENHAHN_STELLAR_TYPE=4 # set stellar type for output, -1=always,-2=disable
+MINT_KIPPENHAHN_COMPANION_STELLAR_TYPE=-2 # see previous
+MINT_MINIMUM_SHELL_MASS=1e-4
+MINT_MAXIMUM_SHELL_MASS=1e-3
+MINT_NUCLEAR_BURNING=1
+MINT_NSHELLS=200
 MINT_REMESH=True
 MINT_MS_REJUVENATION=True
 MINT_DISABLE_GRID_LOAD_WARNINGS=False
@@ -31,9 +32,14 @@ MINT_DATA_CLEANUP=True
 MINT_USE_ZAMS_PROFILES=True
 MINT_FALLBACK_TO_TEST_DATA=True
 MINT_DIR=$HOME/data/MINT/data
+# maximum stellar type allowed : 16 means run them all
+# a negative number means this is ignores if MASSLESS_REMNANT (15)
+MAX_STELLAR_TYPE1=6
+MAX_STELLAR_TYPE2=-16
+
 
 # RUN tbse
-tbse \
+tbse $TBSEARGS \
     --max_evolution_time 15000 \
     --minimum_timestep $MINIMUM_TIMESTEP \
     --maximum_timestep $MAXIMUM_TIMESTEP \
@@ -48,6 +54,7 @@ tbse \
     --MINT_maximum_shell_mass $MINT_MAXIMUM_SHELL_MASS \
     --MINT_nuclear_burning $MINT_NUCLEAR_BURNING \
     --MINT_remesh $MINT_REMESH \
+    --MINT_nshells $MINT_NSHELLS \
     --MINT_MS_rejuvenation $MINT_MS_REJUVENATION \
     --MINT_disable_grid_load_warnings $MINT_DISABLE_GRID_LOAD_WARNINGS \
     --MINT_data_cleanup $MINT_DATA_CLEANUP \
@@ -58,5 +65,7 @@ tbse \
     --magnetic_braking_gamma $MAGNETIC_BRAKING_GAMMA \
     --tides_convective_damping  $TIDES_CONVECTIVE_DAMPING \
     --E2_prescription $E2_PRESCRIPTION \
+    --max_stellar_type_1 $MAX_STELLAR_TYPE1 \
+    --max_stellar_type_2 $MAX_STELLAR_TYPE2 \
      \
     $@