diff --git a/src/nucsyn/nucsyn_burn_helium.c b/src/nucsyn/nucsyn_burn_helium.c
new file mode 100644
index 0000000000000000000000000000000000000000..99466f49651cec8c2d252319f2f236e5d6285c64
--- /dev/null
+++ b/src/nucsyn/nucsyn_burn_helium.c
@@ -0,0 +1,165 @@
+#include "../binary_c.h"
+
+No_empty_translation_unit_warning;
+
+
+#ifdef NUCSYN
+
+static void helium_Jacobian(double *const*const RESTRICT J,
+ const double * RESTRICT const N,
+ const double * RESTRICT const sigmav);
+static void helium_derivatives(const double * RESTRICT const y,
+ double * RESTRICT const dydt,
+ const double * RESTRICT const sigmav,
+ const double * RESTRICT const N);
+
+double nucsyn_burn_helium(struct stardata_t * RESTRICT const stardata,
+ double * RESTRICT const Nin,
+ const double * RESTRICT const sigmav,
+ const double h)
+{
+ double err;
+ double **Jhelium;
+#define NSPECIES 5
+ if(stardata->tmpstore->network_jacobians[NUCSYN_NETWORK_helium] == NULL)
+ {
+ Allocate_network_jacobian(NUCSYN_NETWORK_helium,
+ NSPECIES);
+ }
+
+ Jhelium = stardata->tmpstore->network_jacobians[NUCSYN_NETWORK_helium];
+ double y_heburn[NSPECIES+1],
+ y2_heburn[NSPECIES+1],
+ f1_heburn[NSPECIES+1],
+ f2_heburn[NSPECIES+1],
+ f3_heburn[NSPECIES+1],
+ f4_heburn[NSPECIES+1];
+ int indx_heburn[NSPECIES+1];
+ y_heburn[1] = Nin[XHe4];
+ y_heburn[2] = Nin[XC12];
+ y_heburn[3] = Nin[XO16];
+ y_heburn[4] = Nin[XNe20];
+ y_heburn[5] = Nin[XMg24];
+
+ err=kaps_rentrop_GSL(stardata,
+ h,
+ y_heburn,
+ y2_heburn,
+ NSPECIES,
+ Jhelium,
+ sigmav,
+ Nin,
+ &helium_Jacobian,
+ &helium_derivatives,
+ f1_heburn,
+ f2_heburn,
+ f3_heburn,
+ f4_heburn,
+ stardata->preferences->nucsyn_network_error[NUCSYN_NETWORK_helium],
+ TRUE,
+ indx_heburn);
+ if(err<=1.0)
+ {
+ Nin[XHe4] = y_heburn[1];
+ Nin[XC12] = y_heburn[2];
+ Nin[XO16] = y_heburn[3];
+ Nin[XNe20] = y_heburn[4];
+ Nin[XMg24] = y_heburn[5];
+ }
+ return err;
+}
+
+static void helium_derivatives(const double * RESTRICT const y,
+ double * RESTRICT const dydt,
+ const double * RESTRICT const sigmav,
+ const double * RESTRICT const N MAYBE_UNUSED)
+{
+ const double a = sigmav[SIGMAV_3He4] * Pow3(y[1]);
+ const double b = sigmav[SIGMAV_C12_He4] * y[2] * y[1];
+ 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[3] = b - c;
+ dydt[4] = c - d;
+ dydt[5] = d;
+
+ /*
+ dydt[1] =
+ -3.0*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])
+ - sigmav[SIGMAV_C12_He4] * N[XC12] * N[XHe4];
+ dydt[3] =
+ + sigmav[SIGMAV_C12_He4] * N[XC12] * N[XHe4]
+ - sigmav[SIGMAV_O16_He4] * N[XO16] * N[XHe4];
+ dydt[4] =
+ + sigmav[SIGMAV_O16_He4] * N[XO16] * N[XHe4]
+ - sigmav[SIGMAV_Ne20_He4] * N[XNe20] * N[XHe4];
+ dydt[5] =
+ + sigmav[SIGMAV_Ne20_He4] * N[XNe20] * N[XHe4];
+ */
+}
+
+static void helium_Jacobian(double *const*const RESTRICT J,
+ const double * RESTRICT const N,
+ const double * RESTRICT const sigmav)
+{
+ const double y2 = 3.0 * Pow2(N[XHe4]) * sigmav[SIGMAV_3He4];
+ 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];
+ const double dummy4 = sigmav[SIGMAV_C12_He4] * N[XHe4];
+ const double dummy5 = sigmav[SIGMAV_O16_He4] * N[XHe4];
+ const double dummy6 = sigmav[SIGMAV_Ne20_He4] * N[XHe4];
+
+ J[1][1] =
+ -3.0 * y2
+ -dummy1
+ -dummy2
+ -dummy3;
+ J[1][2] = -dummy4;
+ J[1][3] = -dummy5;
+ J[1][4] = -dummy6;
+ J[1][5] = 0.0;
+
+ J[2][1] =
+ + y2
+ - dummy1;
+ J[2][2] = - dummy4;
+ J[2][3] = - dummy5;
+ J[2][4] = 0.0;
+ J[2][5] = 0.0;
+
+ J[3][1] =
+ + dummy1
+ - dummy2;
+ J[3][2] =
+ + dummy4;
+ J[3][3] =
+ - dummy5;
+ J[3][4] = 0.0;
+ J[3][5] = 0.0;
+
+ J[4][1] =
+ + dummy2
+ - dummy3;
+ J[4][2] = 0.0;
+ J[4][3] =
+ - dummy5;
+ J[4][4] =
+ + dummy6;
+ J[4][5] = 0.0;
+
+ J[5][1] = dummy3;
+ J[5][2] = 0.0;
+ J[5][3] = 0.0;
+ J[5][4] = dummy6;
+ J[5][5] = 0.0;
+}
+
+#endif//NUCSYN
diff --git a/src/nucsyn/nucsyn_burn_ppfast.c b/src/nucsyn/nucsyn_burn_ppfast.c
index 558d89ff27efe02494b529c53bb5d2c9f4407d26..224b3269296c70fce0c0a44c31139bdb70ffb86e 100644
--- a/src/nucsyn/nucsyn_burn_ppfast.c
+++ b/src/nucsyn/nucsyn_burn_ppfast.c
@@ -37,7 +37,7 @@
No_empty_translation_unit_warning;
-#if defined NUCSYN
+#ifdef NUCSYN
/*
* Network code : ppfast
* Automatically generated by make_reaction_network_code.pl
diff --git a/src/nucsyn/nucsyn_hbb.c b/src/nucsyn/nucsyn_hbb.c
index c7898f682c9c167cd7a167c93bcbb9d286731c52..f7972f4ed36975fc9c84480625e62a87ec84d59d 100644
--- a/src/nucsyn/nucsyn_hbb.c
+++ b/src/nucsyn/nucsyn_hbb.c
@@ -21,14 +21,50 @@ void nucsyn_clayton_table(double *mnuc);
void nucsyn_hbb(const double T, /* (log10) temperature */
const double rho, /* density */
const double dt, /* time we burn for (years) */
- Abundance * RESTRICT Xhbb, /* incoming mass fractions */
- struct star_t *star,
- struct stardata_t * stardata
+ Abundance * RESTRICT const Xhbb, /* incoming mass fractions */
+ struct star_t * const star,
+ struct stardata_t * const stardata
)
{
const Nuclear_mass * mnuc = stardata->store->mnuc;
const Nuclear_mass *imnuc = stardata->store->imnuc;
+
+#ifdef __TEST_HE_BURN
+ Abundance * X = New_clear_isotope_array;
+ Number_density * N = New_isotope_array;
+ const double testrho = 1e3;
+
+ /* start with a lot of helium */
+ X[XHe4] = 1.0;
+ X_to_N(imnuc,testrho,N,X,ISOTOPE_ARRAY_SIZE);
+
+ const double timemax = 1000;
+ double time = 0.0;
+ double dtime = timemax/1000.0;
+ double logtemp = log10(3e8);
+
+ while(time < timemax)
+ {
+ nucsyn_burn(logtemp,dtime,N,star,stardata);
+ N_to_X(mnuc,testrho,N,X,ISOTOPE_ARRAY_SIZE);
+ printf("HEBURN %g %g %g %g %g %g %g\n",
+ time,
+ X[XHe4],
+ X[XC12],
+ X[XO16],
+ X[XNe20],
+ X[XMg24],
+ XXsum(stardata,X)
+ );
+ time += dtime;
+ }
+ Safe_free(X);
+ Safe_free(N);
+
+ Flexit;
+#endif//__TEST_HE_BURN
+
#ifdef NUCSYN_HBB_RENORMALIZE_MASS_FRACTIONS
Abundance Xtot;
#endif
diff --git a/src/nucsyn/nucsyn_networks.def b/src/nucsyn/nucsyn_networks.def
index 3b708e1c38d2169ad5efc15ea0bcaee19eb7e4f1..afc80e37690f7696db0287a9c3c4fcb043599db8 100644
--- a/src/nucsyn/nucsyn_networks.def
+++ b/src/nucsyn/nucsyn_networks.def
@@ -25,9 +25,10 @@
X(COLDCNO, TRUE, "Cold CNO",&nucsyn_burn_kaps_rentrop_coldCNO,&nucsyn_burn_lsoda_coldCNO,_CVODE(&nucsyn_burn_cvode_coldCNO)) \
X(HOTCNO, FALSE, "Hot CNO",&nucsyn_burn_hotCNO,NULL,_CVODE(NULL)) \
X(NeNa, FALSE, "NeNa",NULL,NULL,_CVODE(NULL)) \
- X(NeNaMgAlnoleak, FALSE, "NeNaMgAl no leak",NULL,NULL,_CVODE(NULL)) \
- X(NeNaMgAl, TRUE, "NeNaMgAl",&nucsyn_burn_kaps_rentrop_NeNaMgAl,&nucsyn_burn_lsoda_coldCNO,_CVODE(&nucsyn_burn_cvode_NeNaMgAl)) \
- X(NUMBER, FALSE, "",NULL,NULL,NULL)
+ X(NeNaMgAlnoleak, FALSE, "NeNaMgAl no leak",NULL,NULL,_CVODE(NULL)) \
+ X(NeNaMgAl, TRUE, "NeNaMgAl",&nucsyn_burn_kaps_rentrop_NeNaMgAl,&nucsyn_burn_lsoda_coldCNO,_CVODE(&nucsyn_burn_cvode_NeNaMgAl)) \
+ X(helium,TRUE,"helium",&nucsyn_burn_helium,NULL,_CVODE(NULL)) \
+ X(NUMBER, FALSE, "",NULL,NULL,NULL)
diff --git a/src/nucsyn/nucsyn_prototypes.h b/src/nucsyn/nucsyn_prototypes.h
index 4ad79586a35c9172a999061f9b5d3fb0ce915b57..842295885a13cdca7c75a01f851ecbe5dd3766ad 100644
--- a/src/nucsyn/nucsyn_prototypes.h
+++ b/src/nucsyn/nucsyn_prototypes.h
@@ -51,9 +51,9 @@ void nucsyn_long_log(struct stardata_t * RESTRICT const stardata);
void nucsyn_hbb(const double temp,
const double rho,
const double thisdt,
- Abundance * RESTRICT Xhbb,
- struct star_t *star,
- struct stardata_t * stardata);
+ Abundance * RESTRICT const Xhbb,
+ struct star_t *const star,
+ struct stardata_t * const stardata);
void nucsyn_set_nuc_masses(Nuclear_mass * RESTRICT mnuc,
Nuclear_mass * RESTRICT imnuc,
Nuclear_mass * RESTRICT mnuc_anu,
@@ -499,6 +499,12 @@ double nucsyn_burn_lsoda_NeNaMgAl(struct stardata_t * RESTRICT const stardata,
const double * RESTRICT const sigmav,
const double h);
+double nucsyn_burn_helium(struct stardata_t * RESTRICT const stardata,
+ double * RESTRICT const Nin,
+ const double * RESTRICT const sigmav,
+ const double h);
+
+
double Constant_function nucsyn_CN_timescale(const double t9,
diff --git a/src/nucsyn/nucsyn_set_sigmav.c b/src/nucsyn/nucsyn_set_sigmav.c
index 8227bc3040432857f1e643b5584cf21280c1fee9..d179d6ba5e3eb7906e14a67fd01779bbb9112e39 100644
--- a/src/nucsyn/nucsyn_set_sigmav.c
+++ b/src/nucsyn/nucsyn_set_sigmav.c
@@ -1609,14 +1609,15 @@ void nucsyn_set_sigmav(struct stardata_t * const stardata,
* page 67
* formulae on page 151
*/
+ const double sigmav_C12_tmp = -32.123 * tm13;
const double sigmav_C12_E1 =
- 6.66e7 * Pow2(i9) * exp(-32.123 * tm13 - t92/Pow2(4.6)) *
+ 6.66e7 * Pow2(i9) * exp(sigmav_C12_tmp - t92/Pow2(4.6)) *
(1.0 + 2.54*t9 + 1.04*t92 - 0.226*t93)
+
1.39e3 * t15 * exp(-28.930*i9);
const double sigmav_C12_E2 =
- 6.56e7 * Pow2(i9) * exp(-32.123 * tm13 - t92/Pow2(1.3)) *
+ 6.56e7 * Pow2(i9) * exp(sigmav_C12_tmp - t92/Pow2(1.3)) *
(1.0 + 9.23*t9 - 13.7*t92 + 7.4*t93);
const double sigmav_C12_res =