add helium burning

src/nucsyn/nucsyn_burn_helium.c

+#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

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

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

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)
 
 
 

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,

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 =