From 83c75270d3759f0d2e55a68cae08fb8775757e3a Mon Sep 17 00:00:00 2001
From: Robert Izzard <r.izzard@surrey.ac.uk>
Date: Sat, 10 Jul 2021 14:38:45 +0100
Subject: [PATCH] add He4(He4He4,gamma), O16(He4,gamma), Ne20(He4,gamma) cross
sections
---
src/nucsyn/nucsyn_parameters.h | 22 +--
src/nucsyn/nucsyn_set_sigmav.c | 331 +++++++++++++++++++--------------
src/nucsyn/nucsyn_sigmav.def | 107 +++++------
src/nucsyn/nucsyn_sigmav.h | 13 +-
src/setup/version.c | 12 +-
src/tables/table_sigmav.c | 49 +++--
6 files changed, 296 insertions(+), 238 deletions(-)
diff --git a/src/nucsyn/nucsyn_parameters.h b/src/nucsyn/nucsyn_parameters.h
index c7c1d5592..5e247ead4 100644
--- a/src/nucsyn/nucsyn_parameters.h
+++ b/src/nucsyn/nucsyn_parameters.h
@@ -761,25 +761,23 @@
#define NUCSYN_SIGMAV_PRE_INTERPOLATE
/*
- * Resolution in t9, i.e. number of steps in t9 space: 3000 ~
- * maximum 5% error.
+ * Resolution in t9 for the sigmav table
*
- * For safety, i.e. not change in yields >1%, choose 1000.
- *
- * For speed with a ~1% error in yields, choose 100.
- *
- * To reduce the bumpiness in Al26 plots, choose 200.
- *
- * Update: you can go as low as 20 if you define
- * NUCSYN_SIGMAV_INTERPOLATE_LOGSIGMAV (see below)
+ * We have this many points between _MIN and _MID,
+ * then the same again between _MID and _MAX
*/
#define NUCSYN_SIGMAV_INTERPOLATION_RESOLUTION 200
/*
- * Range in t9, 0.001 for >MS, 0.13 is ~ max of T9 in STPAGB stars
+ * Range in t9,
+ * _MIN : 0.001 for >MS
+ * _MID : 0.13 is ~ max of T9 in STPAGB stars
+ * _MAX : 10.0 is ~NSE
*/
#define NUCSYN_SIGMAV_INTERPOLATION_MIN 0.001
-#define NUCSYN_SIGMAV_INTERPOLATION_MAX 0.13
+#define NUCSYN_SIGMAV_INTERPOLATION_MID 0.13
+#define NUCSYN_SIGMAV_INTERPOLATION_MAX 10.0
+
/*
* Define this to interpolate in log(t9) rather than t9
diff --git a/src/nucsyn/nucsyn_set_sigmav.c b/src/nucsyn/nucsyn_set_sigmav.c
index 34c52a2f3..8227bc304 100644
--- a/src/nucsyn/nucsyn_set_sigmav.c
+++ b/src/nucsyn/nucsyn_set_sigmav.c
@@ -9,11 +9,10 @@ No_empty_translation_unit_warning;
* to stderr
*/
//#define REACTION_RATES_TABLE
-#define SIGMAV_STRING_SIZE 100
void nucsyn_set_sigmav(struct stardata_t * const stardata,
struct store_t * const store,
- const double temp, /* log temperature */
+ const double temp, /* log temperature */
Reaction_rate * const sigmav, /* reaction rates */
const double * RESTRICT const multipliers /* reaction rate multipliers (ignored if NULL) */
#ifdef RATES_OF_AMANDA_NE22PG_FIX
@@ -21,31 +20,15 @@ void nucsyn_set_sigmav(struct stardata_t * const stardata,
#endif
)
{
- double t9;
-
#ifdef NUCSYN_SIGMAV_PRE_INTERPOLATE
/*
* if temp < 0 then we do not want to use the
* interpolation table, and we really want -1 * temp
*/
- int called_by_preset;
- if(temp<0.0)
- {
- called_by_preset=TRUE;
- t9=exp10(-temp-9.0);
- }
- else
- {
- called_by_preset=FALSE;
- t9=exp10(temp-9.0);
- }
+ const double t9 = exp10(fabs(temp) - 9.0) ;
+ Boolean called_by_preset = temp < 0.0 ? TRUE : FALSE;
#else
- t9=exp10(temp-9.0);
-#endif
-
-#ifdef REACTION_RATES_TABLE
- char reactions[SIGMAV_SIZE][SIGMAV_STRING_SIZE];
- memset(reactions,0,sizeof(char)*SIGMAV_SIZE*SIGMAV_STRING_SIZE);
+ const double t9 = exp10(temp-9.0);
#endif
/* default to zero, as some rates are not set */
@@ -55,7 +38,7 @@ void nucsyn_set_sigmav(struct stardata_t * const stardata,
if(called_by_preset==FALSE &&
In_range(t9,
NUCSYN_SIGMAV_INTERPOLATION_MIN,
- NUCSYN_SIGMAV_INTERPOLATION_MAX))
+ NUCSYN_SIGMAV_INTERPOLATION_MID))
{
nucsyn_set_sigmav_from_table(stardata,store,t9,sigmav);
}
@@ -87,46 +70,54 @@ void nucsyn_set_sigmav(struct stardata_t * const stardata,
/* speed advances (not pretty, but functional) */
- const double tm13=pow(t9,(-0.333333333333333));
- const double t23=Pow2(tm13);//pow(t9,(-2.00/3.00));
- const double t15=pow(t9,-1.50);
- const double t92=t9*t9;
- const double t93=t92*t9;
- const double t94=t93*t9;
- const double t95=t94*t9;
- const double i9=1.00/t9; // on the assumption that *i9 is quicker than /t9
-
+ 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
#ifdef RATES_OF_AREND_JAN
/* Sorry for the fortran variables */
- const double T9=t9;
- const double TP12=sqrt(T9);
- const double TP13=pow(T9,(1.0e0/3.0e0));
- const double TP14=sqrt(TP12);
- const double TP15=pow(T9,(1.0e0/5.0e0));
- const double TP23=Pow2(TP13);
- const double TP32=T9*TP12;
- const double TP43=T9*TP13;
- const double TP53=TP43*TP13;
- const double TP54=T9*TP14;
- const double TP65=T9*TP15;
- const double TM13=1.0e0/TP13;
- const double TM23=pow(TM13,2);
- const double TM32=1.0e0/TP32;
- const double TM9=1.0e0/T9;
- const double TM12=1.0e0/TP12;
- const double TM54=1.0e0/TP54;
- const double TM65=1.0e0/TP65;
- const double T913=TP13,T923=TP23,T943=TP43,T953=TP53,T915=TP15,T912=sqrt(T9);
- const double T932=pow(T9,3.0/2.0),T972=pow(T9,7.0/2.0);
- const double T935=pow(T9,3.0/5.0),T958=pow(T9,5.0/8.0);
- const double T914=pow(T9,1.0/4.0),T927=pow(T9,2.0/7.0);
+ const double T9 = t9;
+ const double TP12 = sqrt(T9);
+ const double TP13 = pow(T9,(1.0e0/3.0e0));
+ const double TP14 = sqrt(TP12);
+ const double TP15 = pow(T9,(1.0e0/5.0e0));
+ const double TP23 = Pow2(TP13);
+ const double TP32 = T9*TP12;
+ const double TP43 = T9*TP13;
+ const double TP53 = TP43*TP13;
+ const double TP54 = T9*TP14;
+ const double TP65 = T9*TP15;
+ const double TM13 = 1.0e0/TP13;
+ const double TM23 = pow(TM13,2);
+ const double TM32 = 1.0e0/TP32;
+ const double TM9 = 1.0e0/T9;
+ const double TM12 = 1.0e0/TP12;
+ const double TM54 = 1.0e0/TP54;
+ const double TM65 = 1.0e0/TP65;
+ const double T913 = TP13;
+ const double T923 = TP23;
+ const double T943 = TP43;
+ const double T953 = TP53;
+ const double T915 = TP15;
+ const double T912 = sqrt(T9);
+ const double T932 = pow(T9,3.0/2.0);
+ const double T972 = pow(T9,7.0/2.0);
+ const double T935 = pow(T9,3.0/5.0);
+ const double T958 = pow(T9,5.0/8.0);
+ const double T914 = pow(T9,1.0/4.0);
+ const double T927 = pow(T9,2.0/7.0);
#endif // RATES_OF_AREND_JAN
#ifdef RATES_OF_MARIA
- const double t913=pow(t9,0.333333333333333333333333);
- const double t953=pow(t9,1.66666666666666666666);
- const double lt9=log(t9);
+ 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
@@ -147,74 +138,6 @@ void nucsyn_set_sigmav(struct stardata_t * const stardata,
const double t914=pow(t9,0.25);
#endif // RATES_OF_AMANDA
-#ifdef REACTION_RATES_TABLE
- {
- /*
- * Set strings for each reaction
- */
-#define SETR(A,B) strlcpy(reactions[(A)],B,SIGMAV_STRING_SIZE);
-
- /* pp */
- SETR(SIGMAV_PP,"H + H -> D");
- SETR(SIGMAV_He3_He3,"He3 + He3 -> He4 + 2p");
- SETR(SIGMAV_He3_He4,"He3 + He4 -> Be7 + g");
-
- SETR(SIGMAV_Li6_p,"Li6 + p -> Be7");
- SETR(SIGMAV_Li6_p_alpha,"Li6 + p -> He4 + He3");
- SETR(SIGMAV_Li7_p,"Li7 + p -> Be8");
- SETR(SIGMAV_Li7_p_alpha,"Li7 + p -> He4 + He4");
- SETR(SIGMAV_Be7_p,"Be7 + p -> g + B8");
- SETR(SIGMAV_Be7_e,"Be7 + e -> Li7 + nu");
-
- /* CNO */
- SETR(SIGMAV_T12,"C12 + H1 -> N13 + g");
- SETR(SIGMAV_T15,"N15 + H1 -> C12 + He4");
- SETR(SIGMAV_T15_BRANCH,"N15 + H1 -> O16 + g");
- SETR(SIGMAV_TBETA13,"N13 -> C13 + e-");
- SETR(SIGMAV_TPN13,"N13 + p -> O14 + g");
- SETR(SIGMAV_T13,"C13 + H1 -> N14 + g");
- SETR(SIGMAV_T14,"N14 + H1 -> O15 + g");
- SETR(SIGMAV_TO14BETA,"O14 -> N14 + e-");
- SETR(SIGMAV_T17,"O17 + H1 -> N14 + He4");
- SETR(SIGMAV_TBETA15,"O15 -> N15 + e-");
- SETR(SIGMAV_T15P,"O15 + p -> F16");
- SETR(SIGMAV_T16,"O16 + H1 -> F17 + g");
- SETR(SIGMAV_O17G,"O17 + p -> F18 + gamma");
- SETR(SIGMAV_O18A,"O18 + p -> N15 + alpha");
- SETR(SIGMAV_O18G,"O18 + p -> F19 + gamma");
- SETR(SIGMAV_F18BETA,"F18 -> e+ + O18");
- SETR(SIGMAV_F19A,"F19 + p -> O16 + alpha");
- SETR(SIGMAV_F19n,"F19 + p -> Ne19 (maybe -> F18 (-> O18 + e+) +e+) + n");
- SETR(SIGMAV_TBETA17,"F17 -> O17 + e+ + nu");
- SETR(SIGMAV_TF17,"F17 + p -> Ne18 + g");
- SETR(SIGMAV_BRANCH_GAMMA,"ON branching ratio");
- SETR(SIGMAV_BRANCH_ALPHA,"CN branching ratio");
- SETR(SIGMAV_T19,"F19 + p -> Ne20 + g");
-
- /* NeNa */
- SETR(SIGMAV_T20,"Ne20 + p -> Na21 + g");
- SETR(SIGMAV_T21,"Ne21 + p -> Na22 + g");
- SETR(SIGMAV_T22,"Ne22 + p -> Na23 + g"); // reaction 18
- SETR(SIGMAV_T23,"Na23 + p -> Ne20 + alpha");
- SETR(SIGMAV_T23P,"Na23 + p -> Mg24 + gamma");
-
- /* MgAl */
- SETR(SIGMAV_T24,"Mg24 + p -> Al25 + gamma");
- SETR(SIGMAV_T25,"Mg25 + p -> Al26 + gamma");
- SETR(SIGMAV_T25m,"Mg25 + p -> Al26m + gamma");
- SETR(SIGMAV_TBETA26,"Al26 -> Mg26 + e-");
- SETR(SIGMAV_T26P,"Al26 + p -> Si27 + gamma");
- SETR(SIGMAV_T26,"Mg26 + p -> Al27 + gamma");
-
-
- SETR(SIGMAV_T27,"Al27 + p -> Si28 + gamma");
- SETR(SIGMAV_T27a,"Al27 + p -> Mg24 + alpha");
- SETR(SIGMAV_T27an,"Al27 + alpha -> P30 + n");
- SETR(SIGMAV_T28,"Si28 + p -> P29 + gamma");
- }
-#endif //REACTION_RATES_TABLE
-
-
/************************************************************/
/************************************************************/
@@ -488,8 +411,8 @@ void nucsyn_set_sigmav(struct stardata_t * const stardata,
* hence this extra factor
*/
#ifdef NUCSYN_THERMALIZED_CORRECTIONS
- sigmav[SIGMAV_T13]*=(1.00-2.0700*exp(-37.9380*i9));
-#endif
+ sigmav[SIGMAV_T13] *= 1.00-2.0700*exp(-37.9380*i9);
+#endif // NUCSYN_THERMALIZED_CORRECTIONS
#endif // RATE_CHOOSER
sigmav[SIGMAV_T13] /= AVOGADRO_CONSTANT;
@@ -1619,6 +1542,125 @@ void nucsyn_set_sigmav(struct stardata_t * const stardata,
sigmav[SIGMAV_T28]*=(1.0-23.813*exp(-22.872*i9-0.353*t9));
#endif
sigmav[SIGMAV_T28] /= AVOGADRO_CONSTANT;
+
+
+ /*
+ * triple-alpha rate : NACRE
+ */
+
+ /*
+ * He4+He4+n -> Be9 + n
+ * This is NA<sigmav>^{alphalpha}_{gs}
+ */
+ const double sigmav_alpha_alpha =
+ 2.43e9 * t23 * exp(-13.490 * tm13 - t92/Pow2(0.15)) * (1.0 + 74.5*t9)
+ +
+ 6.09e5 * t23 * exp(-1.054 * i9);
+
+ /*
+ * This is NA<sigmav>^{alphaalphan}_{gs}
+
+ const double sigmav_alpha_alpha_n =
+ sigmav_alpha_alpha *
+ t9 < 0.03
+ ?
+ (
+ 6.69e-12 * (1.0 - 192*t9 + 2.48e4*t92 - 1.50e6*t93 + 4.13e7*t94 - 3.90e8 * t95)
+ )
+ :
+ (
+ 2.42e-12 * (1.0 - 1.52*logt9 + 0.448*Pow2(logt9) + 0.435*Pow3(logt9))
+ );
+ */
+
+ /*
+ * He4+He4 -> Be8
+ * This is NA<sigmav>^{alphaBe8}_{gs}
+ */
+ const double sigmav_Be8 =
+ 2.76e7 * t23 * exp(-23.570*tm13 - t92/Pow2(0.4)) *
+ (1.0 + 5.47 * t9 + 326.0 * t92)
+ +
+ 130.7 * t15 * exp(-3.338 * i9)
+ +
+ 2.51e4 * t15 * exp(-20.307 * i9);
+
+ /*
+ * hence triple-alpha rate
+ * Table on p42-43
+ * Equations on p148
+ */
+ const double part3 =
+ t9 <= 0.03
+ ?
+ ( 3.07e-16 * (1.0 - 29.1*t9 + 1308.0*t92) )
+ :
+ ( 3.44e-16 * (1.0 + 0.0158 * pow(t9,-0.65)) );
+
+ sigmav[SIGMAV_3He4] =
+ sigmav_alpha_alpha *
+ sigmav_Be8 *
+ part3;
+
+ sigmav[SIGMAV_3He4] /= Pow2(AVOGADRO_CONSTANT);
+
+ /*
+ * C12 + He4 -> O16 + gamma
+ * page 67
+ * formulae on page 151
+ */
+ const double sigmav_C12_E1 =
+ 6.66e7 * Pow2(i9) * exp(-32.123 * tm13 - 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)) *
+ (1.0 + 9.23*t9 - 13.7*t92 + 7.4*t93);
+
+ const double sigmav_C12_res =
+ 19.2 * t92 * exp(-26.9*i9);
+
+ sigmav[SIGMAV_C12_He4] =
+ sigmav_C12_E1 + sigmav_C12_E2 + sigmav_C12_res;
+
+ sigmav[SIGMAV_C12_He4] /= AVOGADRO_CONSTANT;
+
+ /*
+ * O16 + He4 -> Ne20 + gamma
+ * page 87
+ * formulae on page 157
+ */
+ sigmav[SIGMAV_O16_He4] =
+ 2.68e10 * t23 * exp(-39.760*tm13 - t92/Pow2(1.6))
+ +
+ 51.1 * t15 * exp(-10.32 * i9)
+ +
+ 616.6 * t15 * exp(-12.200 * i9)
+ +
+ 0.41 * pow(t9,2.966) * exp(-11.900 * i9);
+
+ sigmav[SIGMAV_O16_He4] /= AVOGADRO_CONSTANT;
+
+ /*
+ * Ne20 + alpha -> Mg24 + gamma
+ * page 110
+ * formulae on page 160
+ */
+ sigmav[SIGMAV_Ne20_He4] =
+ t9 <= 1.0
+ ?
+ ( 8.72 * pow(t9,-0.532) * exp(-8.995 * i9) )
+ :
+ ( 3.74e2 * pow(t9,2.229) * exp(-12.681 * i9) );
+#ifdef NUCSYN_THERMALIZED_CORRECTIONS
+ sigmav[SIGMAV_Ne20_He4] *=
+ (1.0 - 7.787 * exp(-19.821 * i9 - 0.114 * t9));
+#endif // NUCSYN_THERMALIZED_CORRECTIONS
+
+ sigmav[SIGMAV_Ne20_He4] /= AVOGADRO_CONSTANT;
+
/************************************************************/
/************************************************************/
/************************************************************/
@@ -1642,30 +1684,31 @@ void nucsyn_set_sigmav(struct stardata_t * const stardata,
{
for(Reaction_rate_index i=1;i<SIGMAV_SIZE;i++)
{
- if(!(Fequal(multipliers[i],1.0)))
- {
- sigmav[i]*= multipliers[i];
- }
+ sigmav[i] *= multipliers[i];
}
}
#ifdef REACTION_RATES_TABLE
- FILE *fp = stderr;
+ FILE * fp = stderr;
fprintf(fp,"Reaction rates table : log10 T=%g, T=%g, t9=%g\n",temp,exp10(temp),t9);
- fprintf(fp," Reac# ...............................................................Reaction Rate Rate*6e23 \n");
- for(i=1;i<SIGMAV_SIZE;i++)
- {
- _printf("i = %d\n",i);
- fflush(stdout);
- fprintf(fp,"% 6d %70s %10.3e %10.3e\n",
- i,
- reactions[i],
- sigmav[i],
- sigmav[i]*AVOGADRO_CONSTANT
- );
- }
+ fprintf(fp," Reac# .......enum...............................................................Reaction Rate Rate*6e23 \n");
+
+#undef X
+ unsigned int i = 0;
+#define X(REACTION,STRING) \
+ fprintf(fp, \
+ "%3u %14s %70s %10.3e %10.3e\n", \
+ i, \
+ #REACTION, \
+ STRING, \
+ sigmav[i], \
+ sigmav[i]*AVOGADRO_CONSTANT); \
+ i++;
+ NUCSYN_SIGMAV_CODES;
+#undef X
fprintf(fp,"\n");
- Exit_binary_c_no_stardata(BINARY_C_NORMAL_EXIT,"");
+ fflush(NULL);
+ Exit_binary_c_no_stardata(BINARY_C_NORMAL_EXIT,"Sigmav table exit\n");
#endif /* REACTION_RATES_TABLE */
}
diff --git a/src/nucsyn/nucsyn_sigmav.def b/src/nucsyn/nucsyn_sigmav.def
index b9291c0ed..b1955644e 100644
--- a/src/nucsyn/nucsyn_sigmav.def
+++ b/src/nucsyn/nucsyn_sigmav.def
@@ -2,57 +2,60 @@
#ifndef NUCSYN_SIGMAV_CODES_DEF
#define NUCSYN_SIGMAV_CODES_DEF
-#define NUCSYN_SIGMAV_CODES \
- X( F18BETA ) \
- X( T12 ) \
- X( T15 ) \
- X( TBETA13 ) \
- X( TPN13 ) \
- X( T13 ) \
- X( T14 ) \
- X( T17 ) \
- X( TBETA15 ) \
- X( T15P ) \
- X( T16 ) \
- X( TBETA17 ) \
- X( TF17 ) \
- X( BRANCH_ALPHA ) \
- X( BRANCH_GAMMA ) \
- X( T19 ) \
- X( T20 ) \
- X( T21 ) \
- X( T22 ) \
- X( T23 ) \
- X( T23P ) \
- X( T24 ) \
- X( T25 ) \
- X( T26 ) \
- X( T26P ) \
- X( TBETA26 ) \
- X( O17G ) \
- X( O18A ) \
- X( O18G ) \
- X( F19A ) \
- X( F19n ) \
- X( T15_BRANCH ) \
- X( T27 ) \
- X( T27a ) \
- X( T27an ) \
- X( T28 ) \
- X( T25m ) \
- X( TO14BETA ) \
- X( F18pa ) \
- X( He3_He3 ) \
- X( PP ) \
- X( He3_He4 ) \
- X( Be7_p ) \
- X( Li7_p ) \
- X( DP ) \
- X( He3D ) \
- X( DD ) \
- X( Li6_p ) \
- X( Li6_p_alpha ) \
- X( Li7_p_alpha ) \
- X( Be7_e ) \
+#define NUCSYN_SIGMAV_CODES \
+ X( F18BETA,"F18 -> e+ + O18" ) \
+ X( T12,"C12 + H1 -> N13 + g" ) \
+ X( T15,"N15 + H1 -> C12 + He4" ) \
+ X( TBETA13,"N13 -> C13 + e-" ) \
+ X( TPN13,"N13 + p -> O14 + g" ) \
+ X( T13,"C13 + H1 -> N14 + g" ) \
+ X( T14,"N14 + H1 -> O15 + g" ) \
+ X( T17,"O17 + H1 -> N14 + He4" ) \
+ X( TBETA15,"O15 -> N15 + e-" ) \
+ X( T15P,"O15 + p -> F16" ) \
+ X( T16,"O16 + H1 -> F17 + g" ) \
+ X( TBETA17,"F17 -> O17 + e+ + nu" ) \
+ X( TF17,"F17 + p -> Ne18 + g" ) \
+ X( BRANCH_ALPHA,"CN branching ratio" ) \
+ X( BRANCH_GAMMA,"ON branching ratio" ) \
+ X( T19,"F19 + p -> Ne20 + g" ) \
+ X( T20,"Ne20 + p -> Na21 + g" ) \
+ X( T21,"Ne21 + p -> Na22 + g" ) \
+ X( T22,"Ne22 + p -> Na23 + g" ) \
+ X( T23,"Na23 + p -> Ne20 + alpha" ) \
+ X( T23P,"Na23 + p -> Mg24 + gamma" ) \
+ X( T24,"Mg24 + p -> Al25 + gamma" ) \
+ X( T25,"Mg25 + p -> Al26 + gamma" ) \
+ X( T26 ,"Mg26 + p -> Al27 + gamma") \
+ X( T26P,"Al26 + p -> Si27 + gamma" ) \
+ X( TBETA26 ,"Al26 -> Mg26 + e-") \
+ X( O17G,"O17 + p -> F18 + gamma" ) \
+ X( O18A,"O18 + p -> N15 + alpha" ) \
+ X( O18G,"O18 + p -> F19 + gamma" ) \
+ X( F19A,"F19 + p -> O16 + alpha" ) \
+ X( F19n,"F19 + p -> Ne19 (maybe -> F18 (-> O18 + e+) +e+) + n" ) \
+ X( T15_BRANCH,"N15 + H1 -> O16 + g" ) \
+ X( T27,"Al27 + p -> Si28 + gamma" ) \
+ X( T27a ,"Al27 + p -> Mg24 + alpha") \
+ X( T27an,"Al27 + alpha -> P30 + n" ) \
+ X( T28,"Si28 + p -> P29 + gamma" ) \
+ X( T25m,"Mg25 + p -> Al26m + gamma" ) \
+ X( TO14BETA,"O14 -> N14 + e-" ) \
+ X( F18pa, "F18 + p -> O16 + alpha" ) \
+ X( He3_He3,"He3 + He3 -> He4 + 2p" ) \
+ X( PP,"H + H -> D" ) \
+ X( He3_He4,"He3 + He4 -> Be7 + g" ) \
+ X( Be7_p,"Be7 + p -> g + B8" ) \
+ X( Li7_p,"Li7 + p -> Be8" ) \
+ X( DP, "H2 + H1 -> He3 + beta" ) \
+ X( DD, "H2 + H2 -> He4 + gamma" ) \
+ X( Li6_p,"Li6 + p -> Be7" ) \
+ X( Li6_p_alpha,"Li6 + p -> He4 + He3" ) \
+ X( Li7_p_alpha,"Li7 + p -> He4 + He4" ) \
+ X( Be7_e,"Be7 + e -> Li7 + nu" ) \
+ X( 3He4, "3 * He4 -> C12 + gamma") \
+ X( C12_He4, "C12 + He4 -> O16 + gamma") \
+ X( O16_He4, "O16 + He4 -> Ne20 + gamma") \
+ X( Ne20_He4, "Ne20 + He4 -> Mg24 + gamma")
#endif // NUCSYN_SIGMAV_CODES_DEF
diff --git a/src/nucsyn/nucsyn_sigmav.h b/src/nucsyn/nucsyn_sigmav.h
index 7c3f3a491..8c2370fcb 100644
--- a/src/nucsyn/nucsyn_sigmav.h
+++ b/src/nucsyn/nucsyn_sigmav.h
@@ -11,20 +11,23 @@
* Construct evolution loop codes
*/
#undef X
-#define X(CODE) SIGMAV_##CODE,
-enum { NUCSYN_SIGMAV_CODES } ;
+#define X(CODE,STRING) SIGMAV_##CODE,
+enum {
+ NUCSYN_SIGMAV_CODES
+ SIGMAV_SIZE
+} ;
#undef X
-#define X(CODE) #CODE,
+#define X(CODE,STRING) #CODE,
static char * NUCSYN_SIGMAV_code_macros[] MAYBE_UNUSED = { NUCSYN_SIGMAV_CODES };
#undef X
-#define X(CODE) #CODE,
+#define X(CODE,STRING) #CODE,
static char * NUCSYN_SIGMAV_code_strings[] MAYBE_UNUSED = { NUCSYN_SIGMAV_CODES };
#define NUCSYN_SIGMAV_string(N) Array_string(NUCSYN_SIGMAV_code_strings,(N))
-#define SIGMAV_SIZE Array_size(NUCSYN_SIGMAV_code_strings)
+//#define SIGMAV_SIZE Array_size(NUCSYN_SIGMAV_code_strings)
#undef X
diff --git a/src/setup/version.c b/src/setup/version.c
index 33649ba1d..2cabc7309 100644
--- a/src/setup/version.c
+++ b/src/setup/version.c
@@ -1761,10 +1761,10 @@ void version(struct stardata_t * RESTRICT const stardata)
Show_float_macro(NUCSYN_NETWORK_BURN_MAX_FAILURE_COUNT);
Macrotest(NUCSYN_TPAGB_HBB);
Show_float_macro(HBBTFAC_DEFAULT);
- Macrotest(NUCSYN_HBB_RENORMALIZE_MASS_FRACTIONS)
- Macrotest(NUCSYN_NORMALIZE_NUCLEONS)
- Macrotest2(RATES_OF_AMANDA,
- "Using Amanda's nuclear reaction rates");
+ Macrotest(NUCSYN_HBB_RENORMALIZE_MASS_FRACTIONS);
+ Macrotest(NUCSYN_NORMALIZE_NUCLEONS);
+ Macrotest2(RATES_OF_AMANDA,
+ "Using Amanda's nuclear reaction rates");
Macrotest2(RATES_OF_AMANDA_NE22PG_FIX,"with Ne22(pg) fix");
Macrotest2(RATES_OF_MARIA,"Using Maria's nuclear reaction rates\n");
Macrotest2(RATES_OF_AREND_JAN,"Using Arend Jan's nuclear reaction rates\n");
@@ -1772,12 +1772,14 @@ void version(struct stardata_t * RESTRICT const stardata)
#ifdef NUCSYN_SIGMAV_PRE_INTERPOLATE
Show_int_macro(NUCSYN_SIGMAV_INTERPOLATION_RESOLUTION);
Show_float_macro(NUCSYN_SIGMAV_INTERPOLATION_MIN);
- Show_float_macro(NUCSYN_SIGMAV_INTERPOLATION_MAX);
+ Show_float_macro(NUCSYN_SIGMAV_INTERPOLATION_MID);
Macrotestif2(NUCSYN_SIGMAV_INTERPOLATE_LOGT9,"interpolate log T9");
Macrotestifnot2(NUCSYN_SIGMAV_INTERPOLATE_LOGT9,"interpolate linear T9");
Macrotestif2(NUCSYN_SIGMAV_INTERPOLATE_LOGSIGMAV,"interpolate log sigmav");
Macrotestifnot2(NUCSYN_SIGMAV_INTERPOLATE_LOGSIGMAV,"interpolate linear sigmav");
#endif
+#include "../nucsyn/nucsyn_sigmav.h"
+ Show_enum(SIGMAV_SIZE);
Show_float_macro(SIGMAV_TINY);
Macrotest(NUCSYN_HOT_SIGMAV);
Macrotest(NUCSYN_THERMALIZED_CORRECTIONS);
diff --git a/src/tables/table_sigmav.c b/src/tables/table_sigmav.c
index 6aeac6a55..ea933331c 100644
--- a/src/tables/table_sigmav.c
+++ b/src/tables/table_sigmav.c
@@ -12,43 +12,42 @@ void table_sigmav(struct store_t * store)
* Pre-make a table of sigma-v values for interpolation, rather than
* calculate them each time
*/
- double t9,dt9; // temperature
Reaction_rate_index i;
- /* calculate sigmav_size, allow one extra element for the temperature */
- const size_t sigmav_size = (size_t)(SIGMAV_SIZE+1);
+ const size_t table_width = (size_t)(1 + SIGMAV_SIZE);
Reaction_rate * const data =
Malloc((size_t)
- sizeof(Reaction_rate)*(SIGMAV_SIZE+1)*NUCSYN_SIGMAV_INTERPOLATION_RESOLUTION);
+ sizeof(Reaction_rate)
+ *table_width
+ *2*NUCSYN_SIGMAV_INTERPOLATION_RESOLUTION);
/*
* Loop through t9 values and set the sigmav values
*/
- dt9=(NUCSYN_SIGMAV_INTERPOLATION_MAX-NUCSYN_SIGMAV_INTERPOLATION_MIN)/
+ double t9 = NUCSYN_SIGMAV_INTERPOLATION_MIN;
+ double dt9 =
+ (NUCSYN_SIGMAV_INTERPOLATION_MID - NUCSYN_SIGMAV_INTERPOLATION_MIN)/
(double)(NUCSYN_SIGMAV_INTERPOLATION_RESOLUTION-1);
+ int location = 0;
i=0;
-
- for(t9=NUCSYN_SIGMAV_INTERPOLATION_MIN;
- t9<=NUCSYN_SIGMAV_INTERPOLATION_MAX+TINY;
- t9+=dt9)
+ while(t9 <= NUCSYN_SIGMAV_INTERPOLATION_MAX*(1.0 + TINY))
{
- /* pointer to the sigmav array to be set */
- Reaction_rate *p = data + i * sigmav_size;
+ const double logT = log10(t9);
+ Reaction_rate * p = data + i * table_width;
i++;
- double logT=log10(t9);
/* set the first element to t9 */
#ifdef NUCSYN_SIGMAV_INTERPOLATE_LOGT9
- *p = logT;
+ *(p++) = logT;
#else
- *p = t9;
+ *(p++) = t9;
#endif
nucsyn_set_sigmav(NULL,
store,
- -(logT+9.0), // want log10(T)
- p+1, // fill at next element
+ -(logT+9.0), // log10(T)
+ p, // fill at next element
NULL // no multipliers!
#ifdef RATES_OF_AMANDA_NE22PG_FIX
/* this should be fixed! */
@@ -58,11 +57,22 @@ void table_sigmav(struct store_t * store)
#ifdef NUCSYN_SIGMAV_INTERPOLATE_LOGSIGMAV
/* convert sigmav to log10(sigmav) */
- for(size_t j=1;j<SIGMAV_SIZE+1;j++)
+ for(size_t j=0; j<SIGMAV_SIZE; j++)
{
- *(p+j)=log10(*(p+j));
+ *(p+j) = log10(*(p+j));
}
-#endif
+#endif // NUCSYN_SIGMAV_INTERPOLATE_LOGSIGMAV
+
+ if(location == 0 && t9 > NUCSYN_SIGMAV_INTERPOLATION_MID)
+ {
+ /*
+ * Change resolution above 0.13
+ */
+ location = 1;
+ dt9 = (double)(NUCSYN_SIGMAV_INTERPOLATION_MAX-NUCSYN_SIGMAV_INTERPOLATION_MID)/
+ (double)(NUCSYN_SIGMAV_INTERPOLATION_RESOLUTION-1);
+ }
+ t9 += dt9;
}
NewDataTable_from_Pointer(
@@ -72,7 +82,6 @@ void table_sigmav(struct store_t * store)
SIGMAV_SIZE,
NUCSYN_SIGMAV_INTERPOLATION_RESOLUTION
);
-
}
--
GitLab