From d561a0dfeeb8aa5132627fa420b575086790b04e Mon Sep 17 00:00:00 2001
From: Robert Izzard <r.izzard@surrey.ac.uk>
Date: Sun, 4 Mar 2018 12:22:39 +0000
Subject: [PATCH] added code to follow Lithium-7 according to George Angelou's
requirements
This involved cleaning up quite a lot of old code, esp. timestepping
and relating to loop logic, as well as disabling nucsyn routines'
access to lithium.
---
example_systems | 58 +++
src/binary_c_function_macros.h | 2 +
src/binary_c_stellar_types.h | 58 +--
src/binary_c_structures.h | 16 +-
src/binary_c_types.h | 1 +
src/evolution/check_reject_flags.c | 29 ++
src/evolution/evolution.c | 5 +-
src/evolution/evolution_loop_top.c | 2 +-
src/evolution/evolution_nucsyn.c | 4 -
src/evolution/evolution_prototypes.h | 2 +
src/evolution/evolve_detached_system.c | 10 +-
src/evolution/initialize_every_timestep.c | 44 ++-
src/evolution/stellar_evolution.c | 4 +-
src/logging/log_every_timestep.c | 366 ++++++++++--------
src/nucsyn/nucsyn_angelou.h | 15 +-
src/nucsyn/nucsyn_angelou_lithium.c | 260 +++++++++----
src/nucsyn/nucsyn_macros.h | 19 +-
.../nucsyn_observed_surface_abundances.c | 6 +-
src/nucsyn/nucsyn_parameters.h | 21 +-
src/nucsyn/nucsyn_prototypes.h | 3 -
src/nucsyn/nucsyn_set_1st_dup_abunds.c | 4 +-
src/nucsyn/nucsyn_set_third_dredgeup_abunds.c | 4 +-
src/setup/cmd_line_args.h | 111 ++++--
src/setup/cmd_line_macros.h | 16 +-
src/setup/parse_arguments.c | 38 +-
src/setup/version.c | 6 +-
.../stellar_structure_nucsyn.c | 5 +-
.../stellar_structure_prototypes.h | 3 +
src/timestep/timestep_limits.c | 72 +++-
tbse | 102 ++---
30 files changed, 835 insertions(+), 451 deletions(-)
create mode 100644 example_systems
create mode 100644 src/evolution/check_reject_flags.c
diff --git a/example_systems b/example_systems
new file mode 100644
index 000000000..e257bc00b
--- /dev/null
+++ b/example_systems
@@ -0,0 +1,58 @@
+#!/bin/bash
+############################################################
+#
+# Example binary_c systems
+#
+############################################################
+
+# RLOF checker
+#M1=1;M2=0.9;ORBITAL_PERIOD=0;ORBITAL_SEPARATION=7
+
+# early stripper -> sdO + NS
+#M1=5.2; M2=5.1; ORBITAL_PERIOD=30
+
+# double comenv route?
+#M1=7; M2=1.5; ORBITAL_PERIOD=1500
+
+# novae
+#M1=5.1; M2=5.0; ORBITAL_PERIOD=30
+
+#M1=10; M2=1; ORBITAL_PERIOD=10000
+
+# RGB test
+#M1=2; M2=1; ORBITAL_PERIOD=30
+
+# AGB test
+#M1=2; M2=0.2; ORBITAL_PERIOD=1000
+
+# CBDISC
+
+# RGB
+#M1=2; M2=1; ORBITAL_PERIOD=40; MAX_EVOLUTION_TIME=2000.0; METALLICITY=0.014
+
+# eagb
+#M1=2; M2=1; ORBITAL_PERIOD=1000; MAX_EVOLUTION_TIME=2000; METALLICITY=0.014
+
+# tpagb
+#M1=2; M2=1; ORBITAL_PERIOD=2000; MAX_EVOLUTION_TIME=2000; METALLICITY=0.014
+
+# Karsten's system in NGC 6791 with age 8.3Gyr, [Fe/H]=0.3
+#ORBITAL_SEPARATION=0; METALLICITY=0.03; MAX_EVOLUTION_TIME=8300
+
+# To model the future system
+#M1=1.77; M2=0.22; ORBITAL_PERIOD=1.4; MAX_EVOLUTION_TIME=20000
+
+# To model the system from the start
+# Pure Case A (MS/HG)
+#M1=1.2; M2=1.1; ORBITAL_PERIOD=1.1; MAX_EVOLUTION_TIME=20000
+# Case AB RLOF (HG)
+#M1=1.2; M2=1.1; ORBITAL_PERIOD=1.3; MAX_EVOLUTION_TIME=20000
+# Case B RLOF
+#M1=1.2; M2=1.1; ORBITAL_PERIOD=1.7; MAX_EVOLUTION_TIME=20000
+
+# or start later, with more magnetic braking
+#M1=1.2; M2=1.1; ORBITAL_PERIOD=4.0; MAX_EVOLUTION_TIME=20000
+
+# CEMP binary
+#M1=2; M2=0.75; ORBITAL_PERIOD=3500; MAX_EVOLUTION_TIME=13700; METALLICITY=0.0001
+
diff --git a/src/binary_c_function_macros.h b/src/binary_c_function_macros.h
index e4a8f7670..b7caea2ef 100644
--- a/src/binary_c_function_macros.h
+++ b/src/binary_c_function_macros.h
@@ -111,6 +111,8 @@
/* Convert a C (integer) expression to Boolean type */
#define Boolean_(EXPR) ((EXPR) ? (TRUE) : (FALSE))
+#define Orset(VAR,VAL) ((VAR) = (VAR) == TRUE ? TRUE : (VAL))
+
/*
* Set, unset and append log strings
* N = log number
diff --git a/src/binary_c_stellar_types.h b/src/binary_c_stellar_types.h
index b92735d72..99fd3ec30 100644
--- a/src/binary_c_stellar_types.h
+++ b/src/binary_c_stellar_types.h
@@ -1,27 +1,15 @@
-
-
#ifndef STELLAR_TYPES_H
-
#define STELLAR_TYPES_H
-/*
- * we need binary_c_parameters.h for the minimum and maximum convective masses
- */
#include "binary_c_parameters.h"
#include "binary_c_macros.h"
-/*
- * This file contains some useful macros for the stellar types... This is Rob's
- * attempt to make the code slightly more readable!
- */
-
-/** stellar types - in pretty array form **/
-#define STELLAR_TYPES "Deeply or fully convective low mass MS","Main Sequence","Hertzsprung Gap","First Giant Branch","Core Helium Burning","First Asymptotic Giant Branch","Second Asymptotic Giant Branch","Main Sequence Naked Helium","Hertzsprung Gap Naked Helium","Giant Branch Naked Helium","Helium White Dwarf","Carbon/Oxygen White Dwarf","Oxygen/Neon White Dwarf","Neutron","Black Hole","Massless"
/*
- * Stellar types taken from Hurley et al. "Evolution of Binary Stars and the
+ * Stellar types for binary_c,
+ * taken from Hurley et al. "Evolution of Binary Stars and the
* Effect of Tides on Binary Populations", MNRAS 2002
*/
@@ -33,7 +21,6 @@
#define HG 2
#define HERTZSPRUNG_GAP 2
-/* Best to use "FIRST_GIANT_BRANCH" to avoid confusion */
#define GIANT_BRANCH 3
#define FIRST_GIANT_BRANCH 3
#define CHeB 4
@@ -65,7 +52,28 @@
#define NUMBER_OF_STELLAR_TYPES 16
-/* Array macro */
+/*
+ * Strings for logging
+ */
+#define STELLAR_TYPES \
+ "Deeply or fully convective low mass MS", \
+ "Main Sequence", \
+ "Hertzsprung Gap", \
+ "First Giant Branch", \
+ "Core Helium Burning", \
+ "First Asymptotic Giant Branch", \
+ "Second Asymptotic Giant Branch", \
+ "Main Sequence Naked Helium", \
+ "Hertzsprung Gap Naked Helium", \
+ "Giant Branch Naked Helium", \
+ "Helium White Dwarf", \
+ "Carbon/Oxygen White Dwarf", \
+ "Oxygen/Neon White Dwarf", \
+ "Neutron", \
+ "Black Hole", \
+ "Massless"
+
+/* String macro */
#define Stellar_type_string(S) \
( \
(S)==LOW_MASS_MAIN_SEQUENCE ? "Convective low mass MS" : \
@@ -86,7 +94,7 @@
(S)==MASSLESS_REMNANT ? "Massless Remnant" : \
"Unknown" \
)
-/* Now what follows are macros defining certain types of stars... */
+
/* White dwarf 10,11,12 */
#define WHITE_DWARF(A) (((A)>HeGB)&&((A)<NEUTRON_STAR))
@@ -138,14 +146,9 @@
#define NUCLEAR_BURNING_BINARY ((NUCLEAR_BURNING(stardata->star[0].stellar_type)) || (NUCLEAR_BURNING(stardata->star[1].stellar_type<HeWD)))
-/* Giant filter for WR candidates */
-//#define GIANT_FILTER(A) ((ON_EITHER_GIANT_BRANCH(A)))
-#define GIANT_FILTER(A) (AGB(A))
-
/* Filter for convective envelopes */
-//#define CONVECTIVE_ENVELOPE(A) ((A==LOW_MASS_MS)||(ON_EITHER_GIANT_BRANCH(A)))
-//#define CONVECTIVE_ENVELOPE(A) ((A==LOW_MASS_MS)||(ON_EITHER_GIANT_BRANCH(A))||((A)==HG)||((A)==HeHG))
-#define CONVECTIVE_ENVELOPE(A) ((ON_EITHER_GIANT_BRANCH(A))||((A)==HG)||((A)==HeHG))
+#define CONVECTIVE_ENVELOPE(A) \
+ ((ON_EITHER_GIANT_BRANCH(A))||((A)==HG)||((A)==HeHG))
/**************************************************/
@@ -178,13 +181,14 @@
#define SINGLE_STAR_SYSTEM -1
/************************************************************/
-
-
#define EITHER_STAR_EVOLVING (EVOLVING(0) || EVOLVING(1))
#define EITHER_STAR_MASSLESS (MASSLESS(0) || MASSLESS(1))
#define NEITHER_STAR_MASSLESS (!(EITHER_STAR_MASSLESS))
#define NOTHING_LEFT (MASSLESS(0) && MASSLESS(1))
-#define OSTAR(A) ((spectral_type(&(stardata->star[(A)]))==SPECTRAL_TYPE_O)&&(stardata->star[(A)].stellar_type<TPAGB))
+#define OSTAR(A) \
+ ((spectral_type(&(stardata->star[(A)]))==SPECTRAL_TYPE_O) \
+ && \
+ (stardata->star[(A)].stellar_type<TPAGB))
#endif /* STELLAR_TYPES_H */
diff --git a/src/binary_c_structures.h b/src/binary_c_structures.h
index 8ec071b13..c87ed6585 100644
--- a/src/binary_c_structures.h
+++ b/src/binary_c_structures.h
@@ -256,9 +256,10 @@ struct preferences_t {
#endif // LITHIUM_TABLES
#ifdef NUCSYN_ANGELOU_LITHIUM
int angelou_lithium_decay_function;
-
+
/* start times */
double angelou_lithium_MS_time;
+ double angelou_lithium_LMMS_time;
double angelou_lithium_HG_time;
double angelou_lithium_GB_time;
double angelou_lithium_CHeB_time;
@@ -271,6 +272,7 @@ struct preferences_t {
/* decay times */
double angelou_lithium_MS_decay_time;
+ double angelou_lithium_LMMS_decay_time;
double angelou_lithium_HG_decay_time;
double angelou_lithium_GB_decay_time;
double angelou_lithium_CHeB_decay_time;
@@ -279,6 +281,7 @@ struct preferences_t {
/* mass fraction of lithium when it is created */
double angelou_lithium_MS_massfrac;
+ double angelou_lithium_LMMS_massfrac;
double angelou_lithium_HG_massfrac;
double angelou_lithium_GB_massfrac;
double angelou_lithium_CHeB_massfrac;
@@ -1115,7 +1118,9 @@ struct star_t {
double r_star_postcomenv_mc;
double r_star_comenv_time;
#endif // LOG_COMENV_RSTARS
-
+#ifdef NUCSYN_ANGELOU_LITHIUM
+ Boolean angelou_lithium_boost_this_timestep;
+#endif// NUCSYN_ANGELOU_LITHIUM
#endif /* NUCSYN */
#ifdef STELLAR_COLOURS
@@ -1189,9 +1194,10 @@ struct star_t {
#ifdef PRE_MAIN_SEQUENCE
Boolean PMS;
#endif
-#if defined(NUCSYN) && defined (CN_THICK_DISC)
+#if defined NUCSYN && defined CN_THICK_DISC
Boolean was_blue_straggler;
#endif
+ Boolean reject;
Stellar_type compact_core_type;
};
@@ -1290,7 +1296,8 @@ struct model_t {
int comenv_type; // non-zero if common envelope occurred.
/* Booleans */
-
+ Boolean reject;
+
/* start and end of RLOF */
Boolean RLOF_start,RLOF_end;
@@ -1317,7 +1324,6 @@ struct model_t {
#if (defined(NUCSYN_NOVAE)&&defined(FILE_LOG))
Boolean novalogged;
#endif
-
#ifdef EVOLUTION_SPLITTING
Boolean doingsplit;
Boolean split_last_time;
diff --git a/src/binary_c_types.h b/src/binary_c_types.h
index 5d6d2f4ab..f402abb25 100644
--- a/src/binary_c_types.h
+++ b/src/binary_c_types.h
@@ -48,5 +48,6 @@ typedef unsigned int Boolean;
#define Supernova_type int
#define Caller_id int
#define Stellar_structure_algorithm int
+#define Composition double
#endif
diff --git a/src/evolution/check_reject_flags.c b/src/evolution/check_reject_flags.c
new file mode 100644
index 000000000..cff0d55fd
--- /dev/null
+++ b/src/evolution/check_reject_flags.c
@@ -0,0 +1,29 @@
+#include "../binary_c.h"
+
+Boolean check_reject_flags(struct stardata_t * stardata)
+{
+ /*
+ * Check global and stellar timestep reject flags,
+ * return TRUE if any of them are TRUE.
+ */
+ Boolean reject = FALSE;
+ if(stardata->model.reject == TRUE)
+ {
+ reject = TRUE;
+ }
+ else
+ {
+ Star_number k;
+ STARLOOP(k)
+ {
+ SETstar(k);
+ if(star->reject == TRUE)
+ {
+ reject = TRUE;
+ break;
+ }
+ }
+ }
+
+ return reject;
+}
diff --git a/src/evolution/evolution.c b/src/evolution/evolution.c
index 414b804e3..38a9517b5 100644
--- a/src/evolution/evolution.c
+++ b/src/evolution/evolution.c
@@ -12,7 +12,8 @@
int evolution(struct stardata_t * RESTRICT stardata,
Boolean system_type)
{
- struct model_t * model = &stardata->model;
+ initialize_system_every_timestep(stardata);
+
/*
* Apply stellar mass and angular momentum changes,
* e.g. wind loss, wind accretion, RLOF
@@ -27,6 +28,7 @@ int evolution(struct stardata_t * RESTRICT stardata,
* He-star fail suggests we should, but
* other systems suggest we should not.
*/
+ struct model_t * model = &stardata->model;
if(system_type==DETACHED_SYSTEM_CALL)
{
/*
@@ -68,7 +70,6 @@ int evolution(struct stardata_t * RESTRICT stardata,
NANCHECK(model->dt);
NANCHECK(model->dtm);
-
/*
* Do stellar evolution over the timestep
*/
diff --git a/src/evolution/evolution_loop_top.c b/src/evolution/evolution_loop_top.c
index 45bdbbf5f..9500a6444 100644
--- a/src/evolution/evolution_loop_top.c
+++ b/src/evolution/evolution_loop_top.c
@@ -22,7 +22,7 @@ int evolution_loop_top(struct stardata_t * RESTRICT stardata)
/* do some nucleosynthesis */
#ifdef NUCSYN
- evolution_nucsyn(stardata);
+ //evolution_nucsyn(stardata);
#endif // NUCSYN
diff --git a/src/evolution/evolution_nucsyn.c b/src/evolution/evolution_nucsyn.c
index 611e4f6c9..0140b4fa1 100644
--- a/src/evolution/evolution_nucsyn.c
+++ b/src/evolution/evolution_nucsyn.c
@@ -8,9 +8,5 @@ void evolution_nucsyn(struct stardata_t * RESTRICT stardata)
/* convenient place to do radioactive decay */
nucsyn_radioactive_decay(stardata);
#endif //NUCSYN_RADIOACTIVE_DECAY
-
-#ifdef NUCSYN_ANGELOU_LITHIUM
- nucsyn_angelou_lithium(stardata);
-#endif // NUCSYN_ANGELOU_LITHIUM
#endif //NUCSYN
}
diff --git a/src/evolution/evolution_prototypes.h b/src/evolution/evolution_prototypes.h
index b087f50f6..0ed42cbef 100644
--- a/src/evolution/evolution_prototypes.h
+++ b/src/evolution/evolution_prototypes.h
@@ -60,6 +60,8 @@ void update_the_time(struct stardata_t * stardata);
void update_phase_start_times(struct stardata_t * stardata);
void catch_events(struct stardata_t * stardata);
void erase_events(struct stardata_t * stardata);
+Boolean check_reject_flags(struct stardata_t * stardata);
+void initialize_system_every_timestep(struct stardata_t * RESTRICT stardata);
#endif /* EVOLUTION_PROTOTYPES_H */
diff --git a/src/evolution/evolve_detached_system.c b/src/evolution/evolve_detached_system.c
index 773491b97..434dc8de9 100644
--- a/src/evolution/evolve_detached_system.c
+++ b/src/evolution/evolve_detached_system.c
@@ -52,8 +52,16 @@ int evolve_detached_system(struct stardata_t * RESTRICT stardata,
{
stop = LOOP_NOW;
Dprint("SN true -> stop = LOOP_NOW");
- }
+ }
+ /*
+ * Check for manual rejection
+ */
+ if(check_reject_flags(stardata) == TRUE)
+ {
+ stop = LOOP_REJECT;
+ }
+
/*
* Do Roche lobe and mass checks, determine mass ratios etc.
*/
diff --git a/src/evolution/initialize_every_timestep.c b/src/evolution/initialize_every_timestep.c
index 16a7c02fd..43c1a88b0 100644
--- a/src/evolution/initialize_every_timestep.c
+++ b/src/evolution/initialize_every_timestep.c
@@ -8,18 +8,13 @@
static void initialize_detached_system_every_timestep_star(struct star_t * RESTRICT star);
-void initialize_detached_system_every_timestep(struct stardata_t * RESTRICT stardata)
-
+void initialize_system_every_timestep(struct stardata_t * RESTRICT stardata)
{
- Dprint("Initialize every timestep : masses %g %g\n",
+ Dprint("Initialize system every timestep : masses %g %g\n",
stardata->star[0].mass,
stardata->star[1].mass
);
struct model_t * model=&(stardata->model);
- model->iter = 0;
- model->intpol = 0;
- model->inttry = FALSE;
-
model->prec = FALSE;
model->supernova = FALSE;
model->com = FALSE;
@@ -27,24 +22,45 @@ void initialize_detached_system_every_timestep(struct stardata_t * RESTRICT star
model->coel = FALSE;
model->novae = FALSE;
model->tphys0 = model->time;
+#ifdef NUCSYN
#ifdef NUCSYN_GCE
/* force evolution to resolve at least one yields_dt time */
stardata->preferences->yields_dt = MIN(stardata->preferences->yields_dt,
stardata->model.max_evolution_time);
-#endif
- model->ndonor=0;
- model->naccretor=1;
+#endif//NUCSYN_GCE
+#endif//NUCSYN
+ model->ndonor = 0;
+ model->naccretor = 1;
zero_stellar_derivatives(stardata);
+#ifdef DISCS
+ disc_initialize_every_timestep(stardata);
+#endif
+ Star_number k;
+ STARLOOP(k)
+ {
+ SETstar(k);
+ star->reject = FALSE;
+#if defined NUCSYN && defined NUCSYN_ANGELOU_LITHIUM
+ star->angelou_lithium_boost_this_timestep = FALSE;
+#endif // NUCSYN && NUCSYN_LITHIUM_ANGELOU
+ }
+}
+
+void initialize_detached_system_every_timestep(struct stardata_t * RESTRICT stardata)
+
+{
+ Dprint("Initialize detached system every timestep : masses %g %g\n",
+ stardata->star[0].mass,
+ stardata->star[1].mass
+ );
+ initialize_system_every_timestep(stardata);
+
Star_number k;
STARLOOP(k)
{
SETstar(k);
initialize_detached_system_every_timestep_star(star);
-
}
-#ifdef DISCS
- disc_initialize_every_timestep(stardata);
-#endif
}
static void initialize_detached_system_every_timestep_star(struct star_t * RESTRICT star)
diff --git a/src/evolution/stellar_evolution.c b/src/evolution/stellar_evolution.c
index fa7d0c8c9..881868084 100644
--- a/src/evolution/stellar_evolution.c
+++ b/src/evolution/stellar_evolution.c
@@ -233,7 +233,6 @@ int stellar_evolution(struct stardata_t * RESTRICT stardata,
*/
set_derivatives(star,newstar,model->dtm * 1e6);
-
/*
* Todo : replace the use of rdot with derivative[]
*/
@@ -282,7 +281,7 @@ int stellar_evolution(struct stardata_t * RESTRICT stardata,
Dprint("Check for stellar type change k=%d, kw=%d, kstar=%d, sgl=%d\n",k,stellar_type,star->stellar_type,model->sgl);
- if(stellar_type!=star->stellar_type)
+ if(stellar_type != star->stellar_type)
{
newstar->mass = mt;
newstar->phase_start_core_mass = newstar->core_mass;
@@ -316,6 +315,7 @@ int stellar_evolution(struct stardata_t * RESTRICT stardata,
/*
* Determine stellar evolution timescales
+ * and the timestep required
*/
stellar_timestep(stellar_type,
newstar->age,
diff --git a/src/logging/log_every_timestep.c b/src/logging/log_every_timestep.c
index dd1abe015..15ffe7889 100644
--- a/src/logging/log_every_timestep.c
+++ b/src/logging/log_every_timestep.c
@@ -46,36 +46,37 @@ void log_every_timestep(struct stardata_t * RESTRICT stardata)
Boolean binary MAYBE_UNUSED = NOT(single);
if(0)PRINTF("TESTLOG %30.12e %g %g %d %g %g %g\n",
- stardata->model.time, // 1
- stardata->star[0].mass,//2
- stardata->star[0].core_mass,
- stardata->star[0].stellar_type,//4
- stardata->star[1].mass,
- stardata->common.separation ,//6
- stardata->common.eccentricity//7
+ stardata->model.time, // 1
+ stardata->star[0].mass,//2
+ stardata->star[0].core_mass,
+ stardata->star[0].stellar_type,//4
+ stardata->star[1].mass,
+ stardata->common.separation ,//6
+ stardata->common.eccentricity//7
);
if(0)PRINTF("TESTLOG t=%30.12e M1=%g Mc1=%g k1=%d M2=%g a=%g e=%g R=%g RL=%g dt=%g\n",
- stardata->model.time, // 1
- stardata->star[0].mass,//2
- stardata->star[0].core_mass,
- stardata->star[0].stellar_type,//4
- stardata->star[1].mass,
- stardata->common.separation ,//6
-
- stardata->common.eccentricity//7
- ,stardata->star[0].radius//8
- ,stardata->star[0].roche_radius //9
- ,stardata->model.dt
+ stardata->model.time, // 1
+ stardata->star[0].mass,//2
+ stardata->star[0].core_mass,
+ stardata->star[0].stellar_type,//4
+ stardata->star[1].mass,
+ stardata->common.separation ,//6
+
+ stardata->common.eccentricity//7
+ ,stardata->star[0].radius//8
+ ,stardata->star[0].roche_radius //9
+ ,stardata->model.dt
);
//if(stardata->star[0].stellar_type>=6)
if(0)printf("HRDHRD %d %g %g\n",
- stardata->star[0].stellar_type,
- Teff(0),
- stardata->star[0].luminosity);
+ stardata->star[0].stellar_type,
+ Teff(0),
+ stardata->star[0].luminosity);
-#if defined GAIA && defined STELLAR_COLOURS
+#if defined GAIA && \
+ defined STELLAR_COLOURS
/*
* Colours for Gaia
*/
@@ -339,7 +340,8 @@ void log_every_timestep(struct stardata_t * RESTRICT stardata)
wtts_log(stardata,dt);
#endif // WTTS_LOG
-#if (defined(NUCSYN) && defined(YONG_PROJECT))
+#if defined NUCSYN && \
+ defined YONG_PROJECT
{
/*
* Project to plot the magnesium ratios
@@ -350,8 +352,8 @@ void log_every_timestep(struct stardata_t * RESTRICT stardata)
{ \
PRINTF("YONG%c %g %g\n", \
single?'s':'b', \
- stardata->star[(K)].Xenv[XMg26]/ \
- stardata->star[(K)].Xenv[XMg24], \
+ stardata->star[(K)].Xenv[XMg26]/ \
+ stardata->star[(K)].Xenv[XMg24], \
dtp); \
}
@@ -843,7 +845,8 @@ void log_every_timestep(struct stardata_t * RESTRICT stardata)
nucsyn_network_test(stardata);
#endif
-#if defined ADAPTIVE_RLOF && defined ADAPTIVE_RLOF_LOG
+#if defined ADAPTIVE_RLOF && \
+ defined ADAPTIVE_RLOF_LOG
/* new RLOF log */
if(//(stardata->preferences->RLOF_method!=0)&&
((stardata->star[0].stellar_type<10)||
@@ -1016,7 +1019,8 @@ void log_every_timestep(struct stardata_t * RESTRICT stardata)
-#if DEBUG==1 && defined NUCSYN
+#if DEBUG==1 && \
+ defined NUCSYN
{
Abundance * Xenv=nucsyn_observed_surface_abundances(&(stardata->star[1]));
if(stardata->star[1].stellar_type<10)
@@ -1197,59 +1201,59 @@ void log_every_timestep(struct stardata_t * RESTRICT stardata)
#ifdef STELLAR_COLOURS
if(0){
- double colours[NUMBER_OF_STARS][9]; /* colour data for each star */
- double l[NUMBER_OF_STARS][5]; /* luminosities */
+ double colours[NUMBER_OF_STARS][9]; /* colour data for each star */
+ double l[NUMBER_OF_STARS][5]; /* luminosities */
// if(stardata->star[0].stellar_type<=1)
- if(stardata->model.time<TINY)
- {
- /* Obtain colours for star k */
- Star_number k;
- STARLOOP(k)
+ if(stardata->model.time<TINY)
{
- binary_colours(stardata,
- stardata->store,
- k,
- colours[k],
- colours[k]+1,
- colours[k]+2,
- colours[k]+3,
- colours[k]+4,
- colours[k]+5,
- colours[k]+6,
- colours[k]+7,
- l[k],
- l[k]+1,
- l[k]+2,
- l[k]+3);
+ /* Obtain colours for star k */
+ Star_number k;
+ STARLOOP(k)
+ {
+ binary_colours(stardata,
+ stardata->store,
+ k,
+ colours[k],
+ colours[k]+1,
+ colours[k]+2,
+ colours[k]+3,
+ colours[k]+4,
+ colours[k]+5,
+ colours[k]+6,
+ colours[k]+7,
+ l[k],
+ l[k]+1,
+ l[k]+2,
+ l[k]+3);
- PRINTF("STELLAR_COLOURS# Mu Mb Mv Mi U-B B-V V-I Mbol (star %d)\n",k);
- PRINTF("STELLAR_COLOURS %g %g %g %g %g %g %g %g\n",
- colours[k][5],colours[k][4],colours[k][0],colours[k][6],
- colours[k][1],colours[k][2],colours[k][3],colours[k][7]);
+ PRINTF("STELLAR_COLOURS# Mu Mb Mv Mi U-B B-V V-I Mbol (star %d)\n",k);
+ PRINTF("STELLAR_COLOURS %g %g %g %g %g %g %g %g\n",
+ colours[k][5],colours[k][4],colours[k][0],colours[k][6],
+ colours[k][1],colours[k][2],colours[k][3],colours[k][7]);
- PRINTF("\n\n");
+ PRINTF("\n\n");
- }
+ }
- /* calculate M_x for the binary */
- if(stardata->star[1].mass<=stardata->star[0].mass)
- {
- PRINTF("BINARY_BINARY_COLOURS %g %g %g %g %g %g %g\n",
- 6.34482 - 2.5 * log10(l[1][0]+l[2][0]),
- 6.07549 - 2.5 * log10(l[1][1]+l[2][1]),
- 5.35129 - 2.5 * log10(l[1][2]+l[2][2]),
- 4.60777 - 2.5 * log10(l[1][3]+l[2][3]),
- stardata->star[0].mass,
- stardata->star[1].mass,
- stardata->common.separation
- );
- }
- /* skip the rest of the code */
- stardata->model.time=16000;
+ /* calculate M_x for the binary */
+ if(stardata->star[1].mass<=stardata->star[0].mass)
+ {
+ PRINTF("BINARY_BINARY_COLOURS %g %g %g %g %g %g %g\n",
+ 6.34482 - 2.5 * log10(l[1][0]+l[2][0]),
+ 6.07549 - 2.5 * log10(l[1][1]+l[2][1]),
+ 5.35129 - 2.5 * log10(l[1][2]+l[2][2]),
+ 4.60777 - 2.5 * log10(l[1][3]+l[2][3]),
+ stardata->star[0].mass,
+ stardata->star[1].mass,
+ stardata->common.separation
+ );
+ }
+ /* skip the rest of the code */
+ stardata->model.time=16000;
- }
+ }
}
#endif
@@ -1296,7 +1300,8 @@ void log_every_timestep(struct stardata_t * RESTRICT stardata)
#endif
-#if defined NUCSYN && defined NUCSYN_AND_HRD
+#if defined NUCSYN && \
+ defined NUCSYN_AND_HRD
if(stardata->star[0].stellar_type<10)
{
struct star_t * star = &(stardata->star[0]);
@@ -1431,26 +1436,29 @@ void log_every_timestep(struct stardata_t * RESTRICT stardata)
#endif
-#if defined CGI_EXTENSIONS && defined CGI_EXTENSIONS_HRD
+#if defined CGI_EXTENSIONS && \
+ defined CGI_EXTENSIONS_HRD
/*
* HRD for cgi (web) output
*/
if((stardata->model.log_fp!=NULL)&&
(stardata->model.time>1e-5))
{
- fprintf(stardata->model.log_fp,
- "CGIHRD %d %g %g %d %g %g\n",
- stardata->star[0].stellar_type,
- Teff(0),
- stardata->star[0].luminosity,
- stardata->star[1].stellar_type,
- Teff(1),
- stardata->star[1].luminosity);
+ fprintf(stardata->model.log_fp,
+ "CGIHRD %d %g %g %d %g %g\n",
+ stardata->star[0].stellar_type,
+ Teff(0),
+ stardata->star[0].luminosity,
+ stardata->star[1].stellar_type,
+ Teff(1),
+ stardata->star[1].luminosity);
}
#endif
-#if defined NUCSYN && defined CGI_EXTENSIONS \
- && defined CGI_EXTENSIONS_ABUNDS && defined FILE_LOG
+#if defined NUCSYN && \
+ defined CGI_EXTENSIONS && \
+ defined CGI_EXTENSIONS_ABUNDS && \
+ defined FILE_LOG
/* Extra stuff to the logfile for cgi (web) output */
if((stardata->model.log_fp!=NULL)&&
@@ -1921,7 +1929,8 @@ void log_every_timestep(struct stardata_t * RESTRICT stardata)
#endif //MASSERON
-#if defined NUCSYN && defined CN_THICK_DISC
+#if defined NUCSYN && \
+ defined CN_THICK_DISC
{
/*
* CN logging in thick-disc giants
@@ -1935,9 +1944,9 @@ void log_every_timestep(struct stardata_t * RESTRICT stardata)
* if Z=Z_SMC=0.004
*/
-#define THICK_DISC_STAR (\
- star->stellar_type>=GIANT_BRANCH && \
- star->stellar_type<HeMS && \
+#define THICK_DISC_STAR ( \
+ star->stellar_type>=GIANT_BRANCH && \
+ star->stellar_type<HeMS && \
stardata->model.time > stardata->preferences->thick_disc_end_age && \
stardata->model.time < stardata->preferences->thick_disc_start_age && \
IN_RANGE(loggv[k], \
@@ -1954,15 +1963,15 @@ void log_every_timestep(struct stardata_t * RESTRICT stardata)
if(THICK_DISC_STAR) out = TRUE;
/*
printf("STAR %d, type %d : age = %g cf %g .. %g, logg = %g vs %g .. %g\n",
- k,
- star->stellar_type,
- stardata->model.time,
- stardata->preferences->thick_disc_end_age,
- stardata->preferences->thick_disc_start_age,
- loggv[k],
- stardata->preferences->thick_disc_logg_min,
- stardata->preferences->thick_disc_logg_max
- );
+ k,
+ star->stellar_type,
+ stardata->model.time,
+ stardata->preferences->thick_disc_end_age,
+ stardata->preferences->thick_disc_start_age,
+ loggv[k],
+ stardata->preferences->thick_disc_logg_min,
+ stardata->preferences->thick_disc_logg_max
+ );
*/
if(star->blue_straggler) star->was_blue_straggler=TRUE;
}
@@ -2035,26 +2044,26 @@ void log_every_timestep(struct stardata_t * RESTRICT stardata)
#endif //CN_THICK_DISC
if(0)
- printf("HRDHRD %g %g %g %g %g %g\n",
- stardata->model.time,
- HRdiagram(0),
- HRdiagram(1),
- MAX(stardata->star[0].derivative[DERIVATIVE_STELLAR_MASS_RLOF_GAIN],
- stardata->star[1].derivative[DERIVATIVE_STELLAR_MASS_RLOF_GAIN]));
-/*
- if(stardata->star[0].radius > stardata->star[0].roche_radius)
- {
- struct star_t * star = &(stardata->star[0]);
- printf("DTT %g %g %g %g %g %g\n",
+ printf("HRDHRD %g %g %g %g %g %g\n",
stardata->model.time,
- stardata->common.separation,
- star->radius,
- star->roche_radius,
- star->mass,
- -star->derivative[DERIVATIVE_STELLAR_MASS_RLOF_LOSS]
- );
- // ~0.002/Myr
- }
+ HRdiagram(0),
+ HRdiagram(1),
+ MAX(stardata->star[0].derivative[DERIVATIVE_STELLAR_MASS_RLOF_GAIN],
+ stardata->star[1].derivative[DERIVATIVE_STELLAR_MASS_RLOF_GAIN]));
+/*
+ if(stardata->star[0].radius > stardata->star[0].roche_radius)
+ {
+ struct star_t * star = &(stardata->star[0]);
+ printf("DTT %g %g %g %g %g %g\n",
+ stardata->model.time,
+ stardata->common.separation,
+ star->radius,
+ star->roche_radius,
+ star->mass,
+ -star->derivative[DERIVATIVE_STELLAR_MASS_RLOF_LOSS]
+ );
+ // ~0.002/Myr
+ }
*/
@@ -2091,7 +2100,8 @@ void log_every_timestep(struct stardata_t * RESTRICT stardata)
#endif
-#if defined COMENV_POLYTROPES && defined COMENV_POLYTROPE_LOGGING
+#if defined COMENV_POLYTROPES && \
+ defined COMENV_POLYTROPE_LOGGING
{
struct star_t * star = & stardata->star[0];
if(star->stellar_type>=HERTZSPRUNG_GAP &&
@@ -2205,12 +2215,12 @@ void log_every_timestep(struct stardata_t * RESTRICT stardata)
disc_logging(stardata);
#endif // DISCS
/*
- printf("ACCLOG %g %g %g %g\n",
- stardata->model.time,
- stardata->star[0].mass,
- stardata->star[0].stellar_type==HeMS ? stardata->star[0].mass : stardata->star[0].core_mass,
- stardata->star[0].phase_start_mass
- );
+ printf("ACCLOG %g %g %g %g\n",
+ stardata->model.time,
+ stardata->star[0].mass,
+ stardata->star[0].stellar_type==HeMS ? stardata->star[0].mass : stardata->star[0].core_mass,
+ stardata->star[0].phase_start_mass
+ );
*/
/*
@@ -2235,33 +2245,33 @@ void log_every_timestep(struct stardata_t * RESTRICT stardata)
* Count post-AGB stars for Hans
*/
/*
-#define POST_AGB_MIN_TEFF 6000
-#define POST_AGB_MAX_TEFF 20000
-#define POST_AGB_MAX_MENV 0.1
-#define POST_AGB_STELLAR_TYPE(ST) (ON_GIANT_BRANCH(ST))
+ #define POST_AGB_MIN_TEFF 6000
+ #define POST_AGB_MAX_TEFF 20000
+ #define POST_AGB_MAX_MENV 0.1
+ #define POST_AGB_STELLAR_TYPE(ST) (ON_GIANT_BRANCH(ST))
- Star_number k;
- STARLOOP(k)
- {
- SETstar(k);
- double teff = Teff(k);
- if(POST_AGB_STELLAR_TYPE(star->stellar_type) &&
- star->mass - star->core_mass < POST_AGB_MAX_MENV &&
- IN_RANGE(teff,
- POST_AGB_MIN_TEFF,
- POST_AGB_MAX_TEFF))
- {
- Boolean has_disc = Boolean_(stardata->common.ndiscs > 0);
- PRINTF("POSTAGB %g %g %d %d %g %g %g\n",
- stardata->model.time, // Myr
- dtp,
- star->stellar_type,
- has_disc,
- teff,
- stardata->common.orbital_period, // y
- stardata->common.eccentricity);
- }
- }
+ Star_number k;
+ STARLOOP(k)
+ {
+ SETstar(k);
+ double teff = Teff(k);
+ if(POST_AGB_STELLAR_TYPE(star->stellar_type) &&
+ star->mass - star->core_mass < POST_AGB_MAX_MENV &&
+ IN_RANGE(teff,
+ POST_AGB_MIN_TEFF,
+ POST_AGB_MAX_TEFF))
+ {
+ Boolean has_disc = Boolean_(stardata->common.ndiscs > 0);
+ PRINTF("POSTAGB %g %g %d %d %g %g %g\n",
+ stardata->model.time, // Myr
+ dtp,
+ star->stellar_type,
+ has_disc,
+ teff,
+ stardata->common.orbital_period, // y
+ stardata->common.eccentricity);
+ }
+ }
*/
#ifdef NUCSYN
@@ -2288,12 +2298,58 @@ void log_every_timestep(struct stardata_t * RESTRICT stardata)
#endif // NUCSYN
/*
- PRINTF("ROTATION %g %g %g %g %g\n",
- stardata->model.time,
- stardata->star[0].v_eq,
- stardata->star[0].v_crit_eq,
- stardata->star[1].v_eq,
- stardata->star[1].v_crit_eq);
+ PRINTF("ROTATION %g %g %g %g %g\n",
+ stardata->model.time,
+ stardata->star[0].v_eq,
+ stardata->star[0].v_crit_eq,
+ stardata->star[1].v_eq,
+ stardata->star[1].v_crit_eq);
*/
-}
+
+#if defined NUCSYN && \
+ defined NUCSYN_ANGELOU_LITHIUM && \
+ defined NUCSYN_ANGELOU_LITHIUM_LOGGING
+ /*
+ * Log stellar lithium for George Angelou
+ *
+ * NB only log hydrogen-rich stars
+ */
+ {
+ Star_number k;
+ STARLOOP(k)
+ {
+ SETstar(k);
+ if(star->stellar_type <= TPAGB)
+ {
+ // surface lithium abundance (mass fraction)
+ Abundance XLi = Observed_abundance(star,XLi7);
+ Number_density N[ISOTOPE_ARRAY_SIZE];
+ X_to_N(stardata->store->imnuc,
+ 1.0,
+ N,
+ Observed_surface(star),
+ ISOTOPE_ARRAY_SIZE);
+ Abundance ALi = MAX(1e-30,log10(N[XLi7] / MAX(1e-100,N[XH1])) + 12.0);
+
+ PRINTF("ANGELOU%d %d %g %g %g %g %g %g %g %g %g %g %g\n",
+ star->starnum, // star number
+ star->stellar_type, // 1 stellar type (see stellar_types.h)
+ stardata->model.time, // 2 time
+ stardata->model.dt * stardata->model.probability, // 3 weighting
+ XLi, // 4 mass fraction
+ ALi, // 5 abundance in log(X/H)+12 notation
+ star->mass, // 6 mass/Msun
+ star->luminosity, // 7 L/Lsun
+ star->radius, // 8 R/Rsun
+ logg(star), // 9 logg (cgs)
+ star->v_eq, // 10 equatorial rotation rate v_eq (km/s)
+ star->v_crit_eq, // 11 critical equatorial rotation rate v_eq_crit (km/s)
+ star->v_eq_ratio // 12 v_eq / v_eq_crit
+ );
+ }
+ }
+ }
+#endif // NUCSYN, NUCSYN_ANGELOU_LITHIUM, NUCSYN_ANGELOU_LITHIUM_LOGGING
+
+}
diff --git a/src/nucsyn/nucsyn_angelou.h b/src/nucsyn/nucsyn_angelou.h
index d4034bc89..14c769eb0 100644
--- a/src/nucsyn/nucsyn_angelou.h
+++ b/src/nucsyn/nucsyn_angelou.h
@@ -10,7 +10,8 @@
*/
#define angelou_factor(TYPE,ST) \
( \
- (ST)<=MAIN_SEQUENCE ? stardata->preferences->angelou_lithium_MS_##TYPE : \
+ (ST)==LOW_MASS_MS ? stardata->preferences->angelou_lithium_LMMS_##TYPE : \
+ (ST)==MAIN_SEQUENCE ? stardata->preferences->angelou_lithium_MS_##TYPE : \
(ST)==HERTZSPRUNG_GAP ? stardata->preferences->angelou_lithium_HG_##TYPE : \
(ST)==GIANT_BRANCH ? stardata->preferences->angelou_lithium_GB_##TYPE : \
(ST)==CHeB ? stardata->preferences->angelou_lithium_CHeB_##TYPE : \
@@ -20,4 +21,16 @@
)
+#define NUCSYN_ANGELOU_BOOST_NONE 0
+#define NUCSYN_ANGELOU_BOOST_IN_TIME_RANGE 1
+#define NUCSYN_ANGELOU_BOOST_VROT 2
+#define NUCSYN_ANGELOU_BOOST_VROTFRAC 3
+
+#define NUCSYN_ANGELOU_BOOST_String(N) ( \
+ (N)==NUCSYN_ANGELOU_BOOST_NONE ? "None" : \
+ (N)==NUCSYN_ANGELOU_BOOST_IN_TIME_RANGE ? "Time range" : \
+ (N)==NUCSYN_ANGELOU_BOOST_VROT ? "v_eq > trigger" : \
+ (N)==NUCSYN_ANGELOU_BOOST_VROTFRAC ? "v_eq / v_rot > trigger" : \
+ "Unknown")
+
#endif // NUCSYN_ANGELOU_H
diff --git a/src/nucsyn/nucsyn_angelou_lithium.c b/src/nucsyn/nucsyn_angelou_lithium.c
index c69634583..a88045593 100644
--- a/src/nucsyn/nucsyn_angelou_lithium.c
+++ b/src/nucsyn/nucsyn_angelou_lithium.c
@@ -12,12 +12,16 @@ static void nucsyn_angelou_decay_lithium(
struct star_t * star
);
+static void react_lithium(struct stardata_t * stardata,
+ const double fac,
+ double * X);
+
#include "nucsyn_angelou.h"
/*
* Set debug to 1 for output
*/
-#define ANGELOU_DEBUG 1
+#define ANGELOU_DEBUG 0
/* debug statement */
#define Aprint(...) \
@@ -33,79 +37,112 @@ static void nucsyn_angelou_decay_lithium(
}
-void nucsyn_angelou_lithium(struct stardata_t * stardata)
+void nucsyn_angelou_lithium(struct stardata_t * stardata,
+ struct star_t * star)
{
/*
* Lithium is spiked when boost is TRUE.
*/
- Star_number k;
- STARLOOP(k)
+ int boost = NUCSYN_ANGELOU_BOOST_NONE;
+
+ /*
+ * Boost lithium at a lithium_time relative to the
+ * phase start time (star->stellar_type_tstart).
+ */
+ double lithium_time = angelou_factor(time,
+ star->stellar_type);
+ double phase_time = stardata->model.time - star->stellar_type_tstart;
+
+ if(IS_NOT_ZERO(lithium_time) &&
+ IN_RANGE(lithium_time - phase_time,
+ 0.0,
+ stardata->model.dtm))
{
- SETstar(k);
- Boolean boost = FALSE;
- //Aprint("lithium check star %d\n",k);
-
/*
- * Boost lithium at a given time
+ * We want lithium to spike in this timestep
*/
- double lithium_time = angelou_factor(time,star->stellar_type);
- if(IS_NOT_ZERO(lithium_time) &&
- IN_RANGE(lithium_time - (stardata->model.time - star->stellar_type_tstart),
- 0.0,
- stardata->model.dt))
- {
- /*
- * We want lithium to spike in this timestep
- */
- boost = TRUE;
- Aprint("Boost lithium by time algorithm : t = %g\n",
- stardata->model.time);
- }
+ boost = NUCSYN_ANGELOU_BOOST_IN_TIME_RANGE;
+ Aprint("Boost lithium by time algorithm : t = %g\n",
+ stardata->model.time);
+ }
- /*
- * Boost lithium when the star's rotation rate exceeds some limit
- */
- if(IS_NOT_ZERO(stardata->preferences->angelou_lithium_vrot_trigger) &&
- star->v_eq > stardata->preferences->angelou_lithium_vrot_trigger)
- {
- boost = TRUE;
- Aprint("Boost lithium because v_eq (=%g) > %g\n",
- star->v_eq,
- stardata->preferences->angelou_lithium_vrot_trigger);
- }
+ /*
+ * Boost lithium when the star's rotation rate exceeds some limit
+ */
+
+ if(IS_NOT_ZERO(stardata->preferences->angelou_lithium_vrot_trigger) &&
+ star->v_eq > stardata->preferences->angelou_lithium_vrot_trigger)
+ {
+ boost = NUCSYN_ANGELOU_BOOST_VROT;
+ Aprint("Boost lithium because v_eq (=%g) > %g\n",
+ star->v_eq,
+ stardata->preferences->angelou_lithium_vrot_trigger);
+ }
- /*
- * Boost lithium when the star's rotation rate relative to critical
- * exceeds some limit
- */
- if(IS_NOT_ZERO(stardata->preferences->angelou_lithium_vrotfrac_trigger) &&
- star->v_eq_ratio > stardata->preferences->angelou_lithium_vrotfrac_trigger)
- {
- boost = TRUE;
- Aprint("Boost lithium because v_eq/v_kep (=%g) > %g\n",
- star->v_eq_ratio,
- stardata->preferences->angelou_lithium_vrotfrac_trigger);
- }
+ /*
+ * Boost lithium when the star's rotation rate relative to critical
+ * exceeds some limit
+ */
+ if(IS_NOT_ZERO(stardata->preferences->angelou_lithium_vrotfrac_trigger) &&
+ star->v_eq_ratio > stardata->preferences->angelou_lithium_vrotfrac_trigger)
+ {
+ boost = NUCSYN_ANGELOU_BOOST_VROTFRAC;
+ Aprint("Boost lithium because v_eq/v_kep (=%g) > %g\n",
+ star->v_eq_ratio,
+ stardata->preferences->angelou_lithium_vrotfrac_trigger);
+ }
+ /*
+ * Apply lithium boost
+ */
+ if(boost != NUCSYN_ANGELOU_BOOST_NONE)
+ {
+ star->Xenv[XLi7] = star->Xacc[XLi7] =
+ angelou_factor(massfrac,star->stellar_type);
+
+ Aprint("Boost lithium to env=%g acc=%g (stellar type %d)\n",
+ star->Xenv[XLi7],
+ star->Xacc[XLi7],
+ star->stellar_type);
+ }
+ else if(star->Xenv[XLi7] > NUCSYN_ANGELOU_MINIMUM_Li7)
+ {
/*
- * Apply boost
+ * Otherwise decay...
*/
- if(boost==TRUE)
- {
- star->Xenv[XLi7] = star->Xacc[XLi7] =
- angelou_factor(massfrac,star->stellar_type);
- Aprint("Boost lithium to %g (stellar type %d)\n",
- star->Xenv[XLi7],
- star->stellar_type);
- }
- else if(star->Xenv[XLi7] > NUCSYN_ANGELOU_MINIMUM_Li7)
- {
- /*
- * Otherwise decay...
- */
- nucsyn_angelou_decay_lithium(stardata,star);
- }
+ nucsyn_angelou_decay_lithium(stardata,star);
}
+
+ /*
+ * Set the lithium boost flag which is used to calculate
+ * the timestep. We don't want a short timestep when the boost
+ * happens, only after it stops (or perhaps just before it).
+ */
+ if(boost != NUCSYN_ANGELOU_BOOST_NONE)
+ {
+ star->angelou_lithium_boost_this_timestep = TRUE;
+ if(0)fprintf(stderr,"BOOST star %d at %g reason %s at %g Li7 now %g\n",
+ star->starnum,
+ stardata->model.time,
+ NUCSYN_ANGELOU_BOOST_String(boost),
+ stardata->model.time,
+ star->Xenv[XLi7]
+ );
+ }
+ else
+ {
+ if(0)fprintf(stderr,"NOBOOST star %d at %g Li7 now %g\n",
+ star->starnum,
+ stardata->model.time,
+ star->Xenv[XLi7]
+ );
+ }
+
+ if(0&&star->starnum==0)
+ printf("TIME %g nucsyn done boost? %d prev %d\n",
+ stardata->model.time,
+ star->angelou_lithium_boost_this_timestep,
+ stardata->previous_stardata->star[star->starnum].angelou_lithium_boost_this_timestep);
}
@@ -114,32 +151,91 @@ static void nucsyn_angelou_decay_lithium(
struct star_t * star
)
{
- if(stardata->preferences->angelou_lithium_decay_function ==
- ANGELOU_LITHIUM_DECAY_FUNCTION_EXPONENTIAL)
+
+ /*
+ * Reject long timesteps during the decay phase:
+ *
+ */
+ if(
+ star->Xenv[XLi7] > (1.0-TINY) * NUCSYN_ANGELOU_MINIMUM_Li7 &&
+ stardata->previous_stardata->star[star->starnum].angelou_lithium_boost_this_timestep == TRUE &&
+ stardata->model.dtm > (1.0+TINY) * NUCSYN_ANGELOU_MAX_TIMESTEP_FACTOR * angelou_factor(decay_time,star->stellar_type)
+ )
{
- double decay_time = angelou_factor(decay_time,star->stellar_type);
- if(IS_NOT_ZERO(decay_time))
+ if(0)fprintf(stderr,"RJCT %g %g\n",
+ stardata->model.time,
+ star->Xenv[XLi7]);
+ star->reject = TRUE;
+ }
+ else
+ {
+ if(stardata->preferences->angelou_lithium_decay_function ==
+ ANGELOU_LITHIUM_DECAY_FUNCTION_EXPONENTIAL)
{
- double fac = exp(-stardata->model.dtm/decay_time);
-
- if(!FEQUAL(fac,1.0))
+ double decay_time = angelou_factor(decay_time,
+ star->stellar_type);
+ if(IS_NOT_ZERO(decay_time))
{
- star->Xenv[XLi7] *= fac;
- star->Xacc[XLi7] *= fac;
- Aprint("Decay lithium by %g : is now Xenv=%g Xacc=%g\n",
- fac,
- star->Xenv[XLi7],
- star->Xacc[XLi7]);
+ /*
+ * fac is the "decay" factor: this is integrated
+ * over the timestep (and, for short timesteps, is
+ * approximately exp(-timestep/decay_time) )
+ */
+ double fac = pow(exp(-stardata->model.dtm),1.0/decay_time);
+ if(!FEQUAL(fac,1.0))
+ {
+ react_lithium(stardata,fac,star->Xenv);
+
+ if(star->dmacc > 0.0)
+ react_lithium(stardata,fac,star->Xacc);
+
+ Aprint("Decay lithium by %g : is now Xenv=%g Xacc=%g\n",
+ fac,
+ star->Xenv[XLi7],
+ star->Xacc[XLi7]);
+ }
}
}
+ else
+ {
+ Exit_binary_c(
+ ALGORITHM_OUT_OF_RANGE,
+ "Unknown decay function %d for George Angelou's lithium",
+ stardata->preferences->angelou_lithium_decay_function
+ );
+ }
}
- else
- {
- Exit_binary_c(
- ALGORITHM_OUT_OF_RANGE,
- "Unknown decay function %d for George Angelou's lithium",
- stardata->preferences->angelou_lithium_decay_function
- );
- }
}
+
+
+static void react_lithium(struct stardata_t * stardata,
+ const double fac,
+ double * X)
+{
+ /*
+ * Li7 + p -> He4 + He4
+ *
+ * Convert mass fractions passed in (X)
+ * to number densities (N) at an arbitrary
+ * density (1.0). Then "decay" by factor (fac),
+ * and convert back to mass fractions.
+ */
+ Number_density N[ISOTOPE_ARRAY_SIZE];
+ X_to_N(stardata->store->imnuc,
+ 1.0,
+ N,
+ X,
+ ISOTOPE_ARRAY_SIZE);
+ Number_density was = N[XLi7];
+ N[XLi7] *= fac;
+ Number_density dN = N[XLi7] - was;
+ N[XH1] -= dN;
+ N[XHe4] += 2.0 * dN;
+ N_to_X(stardata->store->mnuc,
+ 1.0,
+ N,
+ X,
+ ISOTOPE_ARRAY_SIZE);
+}
+
#endif // NUCSYN_ANGELOU_LITHIUM
diff --git a/src/nucsyn/nucsyn_macros.h b/src/nucsyn/nucsyn_macros.h
index cab859d74..3660b7659 100644
--- a/src/nucsyn/nucsyn_macros.h
+++ b/src/nucsyn/nucsyn_macros.h
@@ -2,13 +2,28 @@
#ifdef NUCSYN
+/*
+ * Macros for nucleosynthesis in binary_c
+ */
+
#ifndef NUCSYN_MACROS_H
#define NUCSYN_MACROS_H
#include "nucsyn_sources.h"
-/* Compositions should be double precision */
-#define Composition double
+/*
+ * Observed surface abundance
+ */
+#define Observed_surface(S) \
+ ((S)->dmacc>TINY ? &((S)->Xacc[0]) : &((S)->Xenv[0]))
+
+/*
+ * Observed surface abundance of isotope I
+ * in star S
+ */
+#define Observed_abundance(S,I) \
+ (*(Observed_surface(S)+(int)I))
+
/* define conditions for burning of the dredge-up material */
#define DUP_BURN_CONDITIONS(A,B) ((A>=T_MIN_HBB)&&(B>=RHO_MIN_HBB)&&((m0>nucsyn_mhbbmin(stardata->common.metallicity))||(stardata->common.metallicity<=0.004)))
diff --git a/src/nucsyn/nucsyn_observed_surface_abundances.c b/src/nucsyn/nucsyn_observed_surface_abundances.c
index 31bef9b1c..01e494fc9 100644
--- a/src/nucsyn/nucsyn_observed_surface_abundances.c
+++ b/src/nucsyn/nucsyn_observed_surface_abundances.c
@@ -6,9 +6,9 @@ double * nucsyn_observed_surface_abundances(struct star_t *star)
{
/*
* This function selects one of Xenv or Xacc depending on whether
- * there is an accretion layer
+ * there is an accretion layer. The work is now done by the
+ * Observed_surface macro.
*/
- return star->dmacc>TINY ? &(star->Xacc[0]) : &(star->Xenv[0]);
-
+ return Observed_surface(star);
}
#endif
diff --git a/src/nucsyn/nucsyn_parameters.h b/src/nucsyn/nucsyn_parameters.h
index 31f9cb922..1e9aea98a 100644
--- a/src/nucsyn/nucsyn_parameters.h
+++ b/src/nucsyn/nucsyn_parameters.h
@@ -261,9 +261,24 @@ undances with s-process
*/
#define NUCSYN_ANGELOU_LITHIUM
-/* minimum lithium to be considered non-zero */
-#define NUCSYN_ANGELOU_MINIMUM_Li7 1e-50
+/*
+ * minimum lithium to be considered non-zero:
+ * algorithms are disabled below this
+ */
+#define NUCSYN_ANGELOU_MINIMUM_Li7 1e-30
+
+/*
+ * maxmimum timestep factor: use this multiplied
+ * by the decay time to set the maximum timestep
+ * when the lithium abundance exceeds
+ * NUCSYN_ANGELOU_MINIMUM_Li7
+ */
+#define NUCSYN_ANGELOU_MAX_TIMESTEP_FACTOR 0.1
+
+/* turn on logging */
+#define NUCSYN_ANGELOU_LITHIUM_LOGGING
+/* failed combinations */
#ifdef NUCSYN_ANGELOU_LITHIUM
#define REQUIRE_PREVIOUS_STARDATA
#ifdef LITHIUM_TABLES
@@ -271,6 +286,8 @@ You cannot define LITHIUM_TABLES and NUCSYN_ANGELOU_LITHIUM together
#endif
#endif // NUCSYN_ANGELOU_LITHIUM
+
+
// Enable this to check if we should get a late thermal pulse: this
// is not currently implemented so ignore it.
//#define LATE_THERMAL_PULSES
diff --git a/src/nucsyn/nucsyn_prototypes.h b/src/nucsyn/nucsyn_prototypes.h
index 54607e47b..6d20d6b89 100644
--- a/src/nucsyn/nucsyn_prototypes.h
+++ b/src/nucsyn/nucsyn_prototypes.h
@@ -618,9 +618,6 @@ void nucsyn_renormalize_abundance(struct stardata_t * stardata,
Abundance * RESTRICT X);
-#ifdef NUCSYN_ANGELOU_LITHIUM
-void nucsyn_angelou_lithium(struct stardata_t * stardata);
-#endif//NUCSYN_ANGELOU_LITHIUM
#endif /* NUCSYN */
diff --git a/src/nucsyn/nucsyn_set_1st_dup_abunds.c b/src/nucsyn/nucsyn_set_1st_dup_abunds.c
index dbf627068..36b7855a4 100644
--- a/src/nucsyn/nucsyn_set_1st_dup_abunds.c
+++ b/src/nucsyn/nucsyn_set_1st_dup_abunds.c
@@ -385,7 +385,7 @@ void nucsyn_first_dup_apply(double m,double z,
X[XH2]+=r[0];X[XHe3]+=r[1];
-#ifndef LITHIUM_TABLES
+#if !defined LITHIUM_TABLES && !defined NUCSYN_ANGELOU_LITHIUM
X[XLi7]+=r[2];
#endif
@@ -412,7 +412,7 @@ void nucsyn_first_dup_apply(double m,double z,
r[0],r[3],r[7],r[8],r[11],r[15]);
X[XH1]+=r[0]; X[XH2]+=r[1]; X[XHe3]+=r[2]; X[XHe4]+=r[3];
-#ifndef LITHIUM_TABLES
+#if !defined LITHIUM_TABLES && !defined NUCSYN_ANGELOU_LITHIUM
X[XLi7]+=r[4]; X[XLi7]=MAX(0.0,X[XLi7]);
#endif
X[XBe7]+=r[5]; X[XBe7]=MAX(0.0,X[XBe7]); X[XC12]+=r[7]; X[XC13]+=r[8]; X[XN13]+=r[10]; X[XN14]+=r[11]; X[XN15]+=r[12]; X[XO15]+=r[14]; X[XO16]+=r[15]; X[XO17]+=r[16]; X[XO18]+=r[17]; X[XF17]+=r[19]; X[XF18]+=r[20]; X[XF19]+=r[21]; X[XNe20]+=r[24]; X[XNe21]+=r[25]; X[XNe22]+=r[26]; X[XNa22]+=r[28]; X[XNa23]+=r[29];
diff --git a/src/nucsyn/nucsyn_set_third_dredgeup_abunds.c b/src/nucsyn/nucsyn_set_third_dredgeup_abunds.c
index ef8e6b691..dfe06199e 100644
--- a/src/nucsyn/nucsyn_set_third_dredgeup_abunds.c
+++ b/src/nucsyn/nucsyn_set_third_dredgeup_abunds.c
@@ -216,8 +216,10 @@ void nucsyn_set_third_dredgeup_abunds(Abundance * RESTRICT dup_material,
dup_material[XHe3]=intershell[3];
dup_material[XHe4]=intershell[4];
dup_material[XHe4]+=intershell[7]; // CAB: this is B8, decays in He4, half life=0.77s
-#ifdef LITHIUM_TABLES
+#if defined LITHIUM_TABLES
dup_material[XLi7]=0.0;
+#elif defined NUCSYN_LITHIUM_ANGELOU
+ /* do nothing */
#else
dup_material[XLi7]=intershell[5];
#endif
diff --git a/src/setup/cmd_line_args.h b/src/setup/cmd_line_args.h
index 3311fcb2a..dc606295b 100644
--- a/src/setup/cmd_line_args.h
+++ b/src/setup/cmd_line_args.h
@@ -2650,11 +2650,20 @@ struct cmd_line_arg_t
ANGELOU_LITHIUM_DECAY_FUNCTION_VAR, \
1.0 \
}, \
+{ \
+ ARG_SECTION_NUCSYN, \
+ "angelou_lithium_LMMS_time", \
+ "Time at which lithium manufacture is triggered in a low-mass (convective) main sequence (Myr). Requires NUCSYN and NUCSYN_ANGELOU_LITHIUM. Ignored if 0 (for the start, use 1e-6).\n", \
+ ARG_FLOAT, \
+ "0", \
+ ANGELOU_LITHIUM_LMMS_TIME_VAR, \
+ 1.0 \
+ }, \
{ \
ARG_SECTION_NUCSYN, \
"angelou_lithium_MS_time", \
"Time at which lithium manufacture is triggered on the main sequence (Myr). Requires NUCSYN and NUCSYN_ANGELOU_LITHIUM. Ignored if 0 (for the start, use 1e-6).\n", \
- ARG_FLOAT, \
+ ARG_FLOAT, \
"0", \
ANGELOU_LITHIUM_MS_TIME_VAR, \
1.0 \
@@ -2663,7 +2672,7 @@ struct cmd_line_arg_t
ARG_SECTION_NUCSYN, \
"angelou_lithium_HG_time", \
"Time at which lithium manufacture is triggered on the Hertzsprung gap (Myr). Requires NUCSYN and NUCSYN_ANGELOU_LITHIUM. Ignored if 0 (for the start, use 1e-6).\n", \
- ARG_FLOAT, \
+ ARG_FLOAT, \
"0", \
ANGELOU_LITHIUM_HG_TIME_VAR, \
1.0 \
@@ -2672,7 +2681,7 @@ struct cmd_line_arg_t
ARG_SECTION_NUCSYN, \
"angelou_lithium_GB_time", \
"Time at which lithium manufacture is triggered on the giant branch (Myr). Requires NUCSYN and NUCSYN_ANGELOU_LITHIUM. Ignored if 0 (for the start, use 1e-6).\n", \
- ARG_FLOAT, \
+ ARG_FLOAT, \
"0", \
ANGELOU_LITHIUM_GB_TIME_VAR, \
1.0 \
@@ -2681,7 +2690,7 @@ struct cmd_line_arg_t
ARG_SECTION_NUCSYN, \
"angelou_lithium_CHeB_time", \
"Time at which lithium manufacture is triggered during core helium burning (Myr). Requires NUCSYN and NUCSYN_ANGELOU_LITHIUM. Ignored if 0 (for the start, use 1e-6).\n", \
- ARG_FLOAT, \
+ ARG_FLOAT, \
"0", \
ANGELOU_LITHIUM_CHeB_TIME_VAR, \
1.0 \
@@ -2690,7 +2699,7 @@ struct cmd_line_arg_t
ARG_SECTION_NUCSYN, \
"angelou_lithium_EAGB_time", \
"Time at which lithium manufacture is triggered on the early AGB (Myr). Requires NUCSYN and NUCSYN_ANGELOU_LITHIUM. Ignored if 0 (for the start, use 1e-6).\n", \
- ARG_FLOAT, \
+ ARG_FLOAT, \
"0", \
ANGELOU_LITHIUM_EAGB_TIME_VAR, \
1.0 \
@@ -2699,117 +2708,135 @@ struct cmd_line_arg_t
ARG_SECTION_NUCSYN, \
"angelou_lithium_TPAGB_time", \
"Time at which lithium manufacture is triggered on the thermally pulsing AGB (Myr). Requires NUCSYN and NUCSYN_ANGELOU_LITHIUM. Ignored if 0 (for the start, use 1e-6).\n", \
- ARG_FLOAT, \
+ ARG_FLOAT, \
"0", \
ANGELOU_LITHIUM_TPAGB_TIME_VAR, \
1.0 \
}, \
{ \
ARG_SECTION_NUCSYN, \
- "angelou_lithium_MS_decay_time", \
+ "angelou_lithium_LMMS_decay_time", \
+ "Decay time for surface lithium abundance during the low-mass (convective) main sequence (Myr). Requires NUCSYN and NUCSYN_ANGELOU_LITHIUM. Ignored if 0 (for the start, use 1e-6).\n", \
+ ARG_FLOAT, \
+ "0", \
+ ANGELOU_LITHIUM_LMMS_DECAY_TIME_VAR, \
+ 1.0 \
+ }, \
+{ \
+ ARG_SECTION_NUCSYN, \
+ "angelou_lithium_MS_decay_time", \
"Decay time for surface lithium abundance on the main sequence (Myr). Requires NUCSYN and NUCSYN_ANGELOU_LITHIUM. Ignored if 0 (for the start, use 1e-6).\n", \
- ARG_FLOAT, \
+ ARG_FLOAT, \
"0", \
- ANGELOU_LITHIUM_MS_DECAY_TIME_VAR, \
+ ANGELOU_LITHIUM_MS_DECAY_TIME_VAR, \
1.0 \
}, \
{ \
ARG_SECTION_NUCSYN, \
- "angelou_lithium_HG_decay_time", \
+ "angelou_lithium_HG_decay_time", \
"Decay time for surface lithium abundance on the Hertzsprung gap (Myr). Requires NUCSYN and NUCSYN_ANGELOU_LITHIUM. Ignored if 0 (for the start, use 1e-6).\n", \
- ARG_FLOAT, \
+ ARG_FLOAT, \
"0", \
- ANGELOU_LITHIUM_HG_DECAY_TIME_VAR, \
+ ANGELOU_LITHIUM_HG_DECAY_TIME_VAR, \
1.0 \
}, \
{ \
ARG_SECTION_NUCSYN, \
- "angelou_lithium_GB_decay_time", \
+ "angelou_lithium_GB_decay_time", \
"Decay time for surface lithium abundance on the giant branch (Myr). Requires NUCSYN and NUCSYN_ANGELOU_LITHIUM. Ignored if 0 (for the start, use 1e-6).\n", \
- ARG_FLOAT, \
+ ARG_FLOAT, \
"0", \
- ANGELOU_LITHIUM_GB_DECAY_TIME_VAR, \
+ ANGELOU_LITHIUM_GB_DECAY_TIME_VAR, \
1.0 \
}, \
{ \
ARG_SECTION_NUCSYN, \
- "angelou_lithium_CHeB_decay_time", \
+ "angelou_lithium_CHeB_decay_time", \
"Decay time for surface lithium abundance during core helium burning (Myr). Requires NUCSYN and NUCSYN_ANGELOU_LITHIUM. Ignored if 0 (for the start, use 1e-6).\n", \
- ARG_FLOAT, \
+ ARG_FLOAT, \
"0", \
- ANGELOU_LITHIUM_CHeB_DECAY_TIME_VAR, \
+ ANGELOU_LITHIUM_CHeB_DECAY_TIME_VAR, \
1.0 \
}, \
{ \
ARG_SECTION_NUCSYN, \
- "angelou_lithium_EAGB_decay_time", \
+ "angelou_lithium_EAGB_decay_time", \
"Decay time for surface lithium abundance on the early AGB (Myr). Requires NUCSYN and NUCSYN_ANGELOU_LITHIUM. Ignored if 0 (for the start, use 1e-6).\n", \
- ARG_FLOAT, \
+ ARG_FLOAT, \
"0", \
- ANGELOU_LITHIUM_EAGB_DECAY_TIME_VAR, \
+ ANGELOU_LITHIUM_EAGB_DECAY_TIME_VAR, \
1.0 \
}, \
{ \
ARG_SECTION_NUCSYN, \
- "angelou_lithium_TPAGB_decay_time", \
+ "angelou_lithium_TPAGB_decay_time", \
"Decay time for surface lithium abundance on the thermally pulsing AGB (Myr). Requires NUCSYN and NUCSYN_ANGELOU_LITHIUM. Ignored if 0 (for the start, use 1e-6).\n", \
- ARG_FLOAT, \
+ ARG_FLOAT, \
+ "0", \
+ ANGELOU_LITHIUM_TPAGB_DECAY_TIME_VAR, \
+ 1.0 \
+ }, \
+{ \
+ ARG_SECTION_NUCSYN, \
+ "angelou_lithium_LMMS_massfrac", \
+ "Lithium mass fraction when its manufacture is triggered during the low-mass (convective) main sequence (Myr). Requires NUCSYN and NUCSYN_ANGELOU_LITHIUM. Ignored if 0 (for the start, use 1e-6).\n", \
+ ARG_FLOAT, \
"0", \
- ANGELOU_LITHIUM_TPAGB_DECAY_TIME_VAR, \
+ ANGELOU_LITHIUM_LMMS_MASSFRAC_VAR, \
1.0 \
}, \
{ \
ARG_SECTION_NUCSYN, \
- "angelou_lithium_MS_massfrac", \
+ "angelou_lithium_MS_massfrac", \
"Lithium mass fraction when its manufacture is triggered on the main sequence (Myr). Requires NUCSYN and NUCSYN_ANGELOU_LITHIUM. Ignored if 0 (for the start, use 1e-6).\n", \
- ARG_FLOAT, \
+ ARG_FLOAT, \
"0", \
- ANGELOU_LITHIUM_MS_MASSFRAC_VAR, \
+ ANGELOU_LITHIUM_MS_MASSFRAC_VAR, \
1.0 \
}, \
{ \
ARG_SECTION_NUCSYN, \
- "angelou_lithium_HG_massfrac", \
+ "angelou_lithium_HG_massfrac", \
"Lithium mass fraction when its manufacture is triggered on the Hertzsprung gap (Myr). Requires NUCSYN and NUCSYN_ANGELOU_LITHIUM. Ignored if 0 (for the start, use 1e-6).\n", \
- ARG_FLOAT, \
+ ARG_FLOAT, \
"0", \
- ANGELOU_LITHIUM_HG_MASSFRAC_VAR, \
+ ANGELOU_LITHIUM_HG_MASSFRAC_VAR, \
1.0 \
}, \
{ \
ARG_SECTION_NUCSYN, \
- "angelou_lithium_GB_massfrac", \
+ "angelou_lithium_GB_massfrac", \
"Lithium mass fraction when its manufacture is triggered on the giant branch (Myr). Requires NUCSYN and NUCSYN_ANGELOU_LITHIUM. Ignored if 0 (for the start, use 1e-6).\n", \
- ARG_FLOAT, \
+ ARG_FLOAT, \
"0", \
- ANGELOU_LITHIUM_GB_MASSFRAC_VAR, \
+ ANGELOU_LITHIUM_GB_MASSFRAC_VAR, \
1.0 \
}, \
{ \
ARG_SECTION_NUCSYN, \
- "angelou_lithium_CHeB_massfrac", \
+ "angelou_lithium_CHeB_massfrac", \
"Lithium mass fraction when its manufacture is triggered during core helium burning (Myr). Requires NUCSYN and NUCSYN_ANGELOU_LITHIUM. Ignored if 0 (for the start, use 1e-6).\n", \
- ARG_FLOAT, \
+ ARG_FLOAT, \
"0", \
- ANGELOU_LITHIUM_CHeB_MASSFRAC_VAR, \
+ ANGELOU_LITHIUM_CHeB_MASSFRAC_VAR, \
1.0 \
}, \
{ \
ARG_SECTION_NUCSYN, \
- "angelou_lithium_EAGB_massfrac", \
+ "angelou_lithium_EAGB_massfrac", \
"Lithium mass fraction when its manufacture is triggered on the early AGB (Myr). Requires NUCSYN and NUCSYN_ANGELOU_LITHIUM. Ignored if 0 (for the start, use 1e-6).\n", \
- ARG_FLOAT, \
+ ARG_FLOAT, \
"0", \
- ANGELOU_LITHIUM_EAGB_MASSFRAC_VAR, \
+ ANGELOU_LITHIUM_EAGB_MASSFRAC_VAR, \
1.0 \
}, \
{ \
ARG_SECTION_NUCSYN, \
- "angelou_lithium_TPAGB_massfrac", \
+ "angelou_lithium_TPAGB_massfrac", \
"Lithium mass fraction when its manufacture is triggered on the thermally pulsing AGB (Myr). Requires NUCSYN and NUCSYN_ANGELOU_LITHIUM. Ignored if 0 (for the start, use 1e-6).\n", \
- ARG_FLOAT, \
+ ARG_FLOAT, \
"0", \
- ANGELOU_LITHIUM_TPAGB_MASSFRAC_VAR, \
+ ANGELOU_LITHIUM_TPAGB_MASSFRAC_VAR, \
1.0 \
}, \
{ \
diff --git a/src/setup/cmd_line_macros.h b/src/setup/cmd_line_macros.h
index c139be1fd..a003da714 100644
--- a/src/setup/cmd_line_macros.h
+++ b/src/setup/cmd_line_macros.h
@@ -363,6 +363,7 @@
#define ANGELOU_LITHIUM_VROTFRAC_TRIGGER_VAR &(stardata->preferences->angelou_lithium_vrotfrac_trigger)
#define ANGELOU_LITHIUM_MS_TIME_VAR &(stardata->preferences->angelou_lithium_MS_time)
+#define ANGELOU_LITHIUM_LMMS_TIME_VAR &(stardata->preferences->angelou_lithium_LMMS_time)
#define ANGELOU_LITHIUM_HG_TIME_VAR &(stardata->preferences->angelou_lithium_HG_time)
#define ANGELOU_LITHIUM_GB_TIME_VAR &(stardata->preferences->angelou_lithium_GB_time)
#define ANGELOU_LITHIUM_CHeB_TIME_VAR &(stardata->preferences->angelou_lithium_CHeB_time)
@@ -370,13 +371,15 @@
#define ANGELOU_LITHIUM_TPAGB_TIME_VAR &(stardata->preferences->angelou_lithium_TPAGB_time)
#define ANGELOU_LITHIUM_MS_DECAY_TIME_VAR &(stardata->preferences->angelou_lithium_MS_decay_time)
-#define ANGELOU_LITHIUM_HG_DECAY_TIME_VAR &(stardata->preferences->angelou_lithium_MS_decay_time)
-#define ANGELOU_LITHIUM_GB_DECAY_TIME_VAR &(stardata->preferences->angelou_lithium_MS_decay_time)
-#define ANGELOU_LITHIUM_CHeB_DECAY_TIME_VAR &(stardata->preferences->angelou_lithium_MS_decay_time)
-#define ANGELOU_LITHIUM_EAGB_DECAY_TIME_VAR &(stardata->preferences->angelou_lithium_MS_decay_time)
-#define ANGELOU_LITHIUM_TPAGB_DECAY_TIME_VAR &(stardata->preferences->angelou_lithium_MS_decay_time)
+#define ANGELOU_LITHIUM_LMMS_DECAY_TIME_VAR &(stardata->preferences->angelou_lithium_LMMS_decay_time)
+#define ANGELOU_LITHIUM_HG_DECAY_TIME_VAR &(stardata->preferences->angelou_lithium_HG_decay_time)
+#define ANGELOU_LITHIUM_GB_DECAY_TIME_VAR &(stardata->preferences->angelou_lithium_GB_decay_time)
+#define ANGELOU_LITHIUM_CHeB_DECAY_TIME_VAR &(stardata->preferences->angelou_lithium_CHeB_decay_time)
+#define ANGELOU_LITHIUM_EAGB_DECAY_TIME_VAR &(stardata->preferences->angelou_lithium_EAGB_decay_time)
+#define ANGELOU_LITHIUM_TPAGB_DECAY_TIME_VAR &(stardata->preferences->angelou_lithium_TPAGB_decay_time)
#define ANGELOU_LITHIUM_MS_MASSFRAC_VAR &(stardata->preferences->angelou_lithium_MS_massfrac)
+#define ANGELOU_LITHIUM_LMMS_MASSFRAC_VAR &(stardata->preferences->angelou_lithium_LMMS_massfrac)
#define ANGELOU_LITHIUM_HG_MASSFRAC_VAR &(stardata->preferences->angelou_lithium_HG_massfrac)
#define ANGELOU_LITHIUM_GB_MASSFRAC_VAR &(stardata->preferences->angelou_lithium_GB_massfrac)
#define ANGELOU_LITHIUM_CHeB_MASSFRAC_VAR &(stardata->preferences->angelou_lithium_CHeB_massfrac)
@@ -385,6 +388,7 @@
#else
#define ANGELOU_LITHIUM_DECAY_FUNCTION_VAR NULL
+#define ANGELOU_LITHIUM_LMMS_TIME_VAR NULL
#define ANGELOU_LITHIUM_MS_TIME_VAR NULL
#define ANGELOU_LITHIUM_HG_TIME_VAR NULL
#define ANGELOU_LITHIUM_GB_TIME_VAR NULL
@@ -392,11 +396,13 @@
#define ANGELOU_LITHIUM_EAGB_TIME_VAR NULL
#define ANGELOU_LITHIUM_TPAGB_TIME_VAR NULL
#define ANGELOU_LITHIUM_MS_MASSFRAC_VAR NULL
+#define ANGELOU_LITHIUM_LMMS_MASSFRAC_VAR NULL
#define ANGELOU_LITHIUM_HG_MASSFRAC_VAR NULL
#define ANGELOU_LITHIUM_GB_MASSFRAC_VAR NULL
#define ANGELOU_LITHIUM_CHeB_MASSFRAC_VAR NULL
#define ANGELOU_LITHIUM_EAGB_MASSFRAC_VAR NULL
#define ANGELOU_LITHIUM_TPAGB_MASSFRAC_VAR NULL
+#define ANGELOU_LITHIUM_LMMS_DECAY_TIME_VAR NULL
#define ANGELOU_LITHIUM_MS_DECAY_TIME_VAR NULL
#define ANGELOU_LITHIUM_HG_DECAY_TIME_VAR NULL
#define ANGELOU_LITHIUM_GB_DECAY_TIME_VAR NULL
diff --git a/src/setup/parse_arguments.c b/src/setup/parse_arguments.c
index 8404cfbba..a88b241b5 100644
--- a/src/setup/parse_arguments.c
+++ b/src/setup/parse_arguments.c
@@ -210,21 +210,31 @@ void parse_arguments(int start,
/*
* Special cases
*/
- if(IS_ZERO(stardata->common.orbital_period) &&
- IS_NOT_ZERO(stardata->star[0].mass) &&
- IS_NOT_ZERO(stardata->star[1].mass) &&
- IS_NOT_ZERO(stardata->common.separation))
- {
- stardata->common.orbital_period =
- calculate_orbital_period(stardata);
+ if(IS_NOT_ZERO(stardata->star[0].mass) &&
+ IS_NOT_ZERO(stardata->star[1].mass))
+ {
+ if(IS_ZERO(stardata->common.orbital_period) &&
+ IS_NOT_ZERO(stardata->common.separation))
+ {
+ /*
+ * Period is set to zero, but separation is given.
+ * Hence calculate the initial period.
+ */
+ stardata->common.orbital_period =
+ calculate_orbital_period(stardata);
+ }
+ else if(IS_NOT_ZERO(stardata->common.orbital_period))
+ {
+ /*
+ * Update the separation here so that the "INITIAL" log is
+ * correct. (This is reset later anyway.)
+ */
+ stardata->common.orbital_period /= YEAR_LENGTH_IN_DAYS;
+ stardata->common.separation =
+ calculate_orbital_separation(stardata);
+ stardata->common.orbital_period *= YEAR_LENGTH_IN_DAYS;
+ }
}
-
-
- /*
- * Hence the separation in R_SUN : now done later
- */
- // stardata->common.separation =
- // calculate_orbital_separation(stardata);
}
}
diff --git a/src/setup/version.c b/src/setup/version.c
index 4534b4de6..7147b7016 100644
--- a/src/setup/version.c
+++ b/src/setup/version.c
@@ -742,6 +742,11 @@ void version(struct stardata_t * RESTRICT stardata)
Macrotest(NUCSYN_THERMALIZED_CORRECTIONS);
Macrotest(HELIUM_BURNING_REACTIONS);
Macrotest(NUCSYN_S_PROCESS_LOG);
+
+ Macrotest(NUCSYN_LITHIUM);
+ Macrotest(NUCSYN_LITHIUM_TABLES);
+ Macrotest(NUCSYN_ANGELOU_LITHIUM);
+
Macrotest(NUCSYN_R_PROCESS);
Macrotest(NUCSYN_R_PROCESS_ARLANDINI1999);
Macrotest(NUCSYN_R_PROCESS_SIMMERER2004);
@@ -803,7 +808,6 @@ void version(struct stardata_t * RESTRICT stardata)
Show_int_macro(GCE_ISOTOPES);
#endif
Macrotest(NUCSYN_GCE_ELSE);
-
Macrotest(NUCSYN_GCE_OUTFLOW_CHECKS);
#ifdef NUCSYN_GCE_OUTFLOW_CHECKS
Show_float_macro(NUCSYN_GCE_OUTFLOW_FACTOR);
diff --git a/src/stellar_structure/stellar_structure_nucsyn.c b/src/stellar_structure/stellar_structure_nucsyn.c
index cd5c4388a..caefaf873 100644
--- a/src/stellar_structure/stellar_structure_nucsyn.c
+++ b/src/stellar_structure/stellar_structure_nucsyn.c
@@ -131,13 +131,16 @@ void stellar_structure_nucsyn(Stellar_type stellar_type,
}
#endif
+#ifdef NUCSYN_ANGELOU_LITHIUM
+ nucsyn_angelou_lithium(stardata,star);
+#endif // NUCSYN_ANGELOU_LITHIUM
+
#ifdef NUCSYN_STRIP_AND_MIX
/*
* First dredge up and/or envelope stripping
* using TAMS abundances
*/
nucsyn_strip_and_mix(stardata,star);
-
#endif // NUCSYN_STRIP_AND_MIX
diff --git a/src/stellar_structure/stellar_structure_prototypes.h b/src/stellar_structure/stellar_structure_prototypes.h
index d0ad8613d..6ed4a78b1 100644
--- a/src/stellar_structure/stellar_structure_prototypes.h
+++ b/src/stellar_structure/stellar_structure_prototypes.h
@@ -187,4 +187,7 @@ void stellar_structure_modified_BSE(
Supernova_type * new_SN_type
);
+void nucsyn_angelou_lithium(struct stardata_t * stardata,
+ struct star_t * star);
+
#endif /*STELLAR_STRUCTURE_PROTOTYPES_H*/
diff --git a/src/timestep/timestep_limits.c b/src/timestep/timestep_limits.c
index 4d0ebc90b..70ee71e11 100644
--- a/src/timestep/timestep_limits.c
+++ b/src/timestep/timestep_limits.c
@@ -146,7 +146,7 @@ void timestep_limits(Timestep_prototype_args)
}
#endif
-#if (defined RESOLVE_PN_FAST_WIND && defined PN_FAST_WIND)
+#if defined RESOLVE_PN_FAST_WIND && defined PN_FAST_WIND
/*
* Reduce timestep to better resolve the fast wind at the end of the AGB
* or in post-common envelope.
@@ -348,15 +348,71 @@ void timestep_limits(Timestep_prototype_args)
#if defined NUCSYN && defined NUCSYN_ANGELOU_LITHIUM
- if(star->Xenv[XLi7] > NUCSYN_ANGELOU_MINIMUM_Li7)
+
{
- double decay_time = angelou_factor(decay_time,star->stellar_type);
- if(IS_NOT_ZERO(decay_time))
+ double lithium_time = angelou_factor(time,star->stellar_type);
+ struct star_t * previous = &stardata->previous_stardata->star[star->starnum];
+
+ if(IS_NOT_ZERO(lithium_time))
{
- Limit_timestep(*dt,
- 0.1 * decay_time,
- star,
- DT_LIMIT_NUCSYN_ANGELOU_LITHIUM);
+ double dtli = fabs(lithium_time - stardata->model.time + star->stellar_type_tstart);
+ if(IS_NOT_ZERO(dtli))
+ {
+ Limit_timestep(*dt,
+ MAX(1e-3,0.1 * dtli),
+ star,
+ DT_LIMIT_NUCSYN_ANGELOU_LITHIUM);
+ }
+ }
+
+ if(0)
+ fprintf(stderr,
+ ">>>TIMESTEP star %d at time %g this %d prev %d XLi7 %g\n",
+ star->starnum,
+ stardata->model.time,
+ star->angelou_lithium_boost_this_timestep,
+ previous->angelou_lithium_boost_this_timestep,
+ Observed_abundance(star,XLi7));
+
+ if(
+ star->angelou_lithium_boost_this_timestep == FALSE
+ &&
+ (
+ previous->angelou_lithium_boost_this_timestep == TRUE
+ ||
+ Observed_abundance(star,XLi7) > (1.0+TINY)*NUCSYN_ANGELOU_MINIMUM_Li7)
+ )
+ {
+ /*
+ * You might think we should only apply this timestep
+ * limiter if the Li7 abundance spikes, but it *may*
+ * spike in the next timestep, we just don't know.
+ */
+ double decay_time =
+ angelou_factor(decay_time,star->stellar_type);
+
+ if(IS_NOT_ZERO(decay_time))
+ {
+ Limit_timestep(*dt,
+ (1.0-TINY) * NUCSYN_ANGELOU_MAX_TIMESTEP_FACTOR * decay_time,
+ star,
+ DT_LIMIT_NUCSYN_ANGELOU_LITHIUM);
+ }
+
+ if(0)
+ fprintf(stderr,
+ "USE DECAY TIMESTEP star %d at time %g use decay time %g -> %g (this %d prev %d >lowerli7 %d [Li7=%g min=%g] reject %d)\n",
+ star->starnum,
+ stardata->model.time,
+ NUCSYN_ANGELOU_MAX_TIMESTEP_FACTOR * decay_time,
+ *dt,
+ star->angelou_lithium_boost_this_timestep,
+ previous->angelou_lithium_boost_this_timestep,
+ (int)(Observed_abundance(star,XLi7) > (1.0+TINY)*NUCSYN_ANGELOU_MINIMUM_Li7),
+ Observed_abundance(star,XLi7),
+ (1.0+TINY)*NUCSYN_ANGELOU_MINIMUM_Li7,
+ star->reject
+ );
}
}
#endif
diff --git a/tbse b/tbse
index 24abcd26c..96bda2595 100755
--- a/tbse
+++ b/tbse
@@ -75,6 +75,7 @@ M2=1.0
# (this requires that the masses are already set)
ORBITAL_PERIOD=0
ORBITAL_SEPARATION=2000
+# TODO fix the initial separation in the log : it is wrong!
# Orbital eccentricity should be in the range 0.0-1.0.
ECCENTRICITY=0.0
@@ -84,7 +85,7 @@ MAX_EVOLUTION_TIME=15000
# metallicity for stellar evolution, must be in the range 1e-4<=Z<=0.03
# (mass fraction)
-METALLICITY=0.0001
+METALLICITY=0.02
# Stellar structure algorithm
# 0 = modified BSE (default)
@@ -93,63 +94,10 @@ METALLICITY=0.0001
# 3 = binary_c (not yet implemented)
STELLAR_STRUCTURE_ALGORITHM=0
-############################################################
-#
-# Example systems
-#
-############################################################
-
-# RLOF checker
-#M1=1;M2=0.9;ORBITAL_PERIOD=0;ORBITAL_SEPARATION=7
-
-# early stripper -> sdO + NS
-#M1=5.2; M2=5.1; ORBITAL_PERIOD=30
-
-# double comenv route?
-#M1=7; M2=1.5; ORBITAL_PERIOD=1500
-
-# novae
-#M1=5.1; M2=5.0; ORBITAL_PERIOD=30
-
-#M1=10; M2=1; ORBITAL_PERIOD=10000
-
-# RGB test
-M1=2; M2=1; ORBITAL_PERIOD=30
-
-# AGB test
-#M1=2; M2=0.2; ORBITAL_PERIOD=1000
-
-# CBDISC
-
-# RGB
-#M1=2; M2=1; ORBITAL_PERIOD=40; MAX_EVOLUTION_TIME=2000.0; METALLICITY=0.014
-
-# eagb
-#M1=2; M2=1; ORBITAL_PERIOD=1000; MAX_EVOLUTION_TIME=2000; METALLICITY=0.014
-
-# tpagb
-M1=2; M2=1; ORBITAL_PERIOD=2000; MAX_EVOLUTION_TIME=2000; METALLICITY=0.014
-
-# Karsten's system in NGC 6791 with age 8.3Gyr, [Fe/H]=0.3
-ORBITAL_SEPARATION=0; METALLICITY=0.03; MAX_EVOLUTION_TIME=8300
-
-# To model the future system
-#M1=1.77; M2=0.22; ORBITAL_PERIOD=1.4; MAX_EVOLUTION_TIME=20000
-
-# To model the system from the start
-# Pure Case A (MS/HG)
-#M1=1.2; M2=1.1; ORBITAL_PERIOD=1.1; MAX_EVOLUTION_TIME=20000
-# Case AB RLOF (HG)
-M1=1.2; M2=1.1; ORBITAL_PERIOD=1.3; MAX_EVOLUTION_TIME=20000
-# Case B RLOF
-#M1=1.2; M2=1.1; ORBITAL_PERIOD=1.7; MAX_EVOLUTION_TIME=20000
-
-# or start later, with more magnetic braking
-M1=1.2; M2=1.1; ORBITAL_PERIOD=4.0; MAX_EVOLUTION_TIME=20000
-
-# CEMP binary
-#M1=2; M2=0.75; ORBITAL_PERIOD=3500; MAX_EVOLUTION_TIME=13700; METALLICITY=0.0001
-
+# possible example systems
+if [ -f "example_systems" ]; then
+ source "example_systems"
+fi
############################################################
# ZAMS rotation rates (km/s)
@@ -157,8 +105,8 @@ M1=1.2; M2=1.1; ORBITAL_PERIOD=4.0; MAX_EVOLUTION_TIME=20000
# 0 = H02 formula for vrot (assuming the star is spherical)
# -1 to start at breakup (which assumes the star is equatorially deformed)
# -2 to sync the angular velocity with the orbit
-VROT1=-1
-VROT2=-1
+VROT1=0
+VROT2=0
# random number seed (ignored if commented out)
#IDUM=-30291
@@ -990,7 +938,8 @@ angelou_lithium_decay_function=0
# perhaps trigger lithium at a certain time, relative to the
# start of the stellar evolutionary phase?
# start times : 0 = ignored, to set at the start use 1e-6
-angelou_lithium_MS_time=500.0
+angelou_lithium_LMMS_time=0.0
+angelou_lithium_MS_time=0.0
angelou_lithium_HG_time=0.0
angelou_lithium_GB_time=0.0
angelou_lithium_CHeB_time=0.0
@@ -1005,20 +954,24 @@ angelou_lithium_vrotfrac_trigger=0.0
# decay times as a function of stellar phase (Myr)
# ignored if 0
-angelou_lithium_MS_decay_time=0.1
-angelou_lithium_HG_decay_time=0.1
-angelou_lithium_GB_decay_time=0.1
-angelou_lithium_CHeB_decay_time=0.1
-angelou_lithium_EAGB_decay_time=0.1
-angelou_lithium_TPAGB_decay_time=0.1
+angelou_decay_time=0.1
+angelou_lithium_LMMS_decay_time=$angelou_decay_time
+angelou_lithium_MS_decay_time=$angelou_decay_time
+angelou_lithium_HG_decay_time=$angelou_decay_time
+angelou_lithium_GB_decay_time=$angelou_decay_time
+angelou_lithium_CHeB_decay_time=$angelou_decay_time
+angelou_lithium_EAGB_decay_time=$angelou_decay_time
+angelou_lithium_TPAGB_decay_time=$angelou_decay_time
# mass fraction of lithium when it is created
-angelou_lithium_MS_massfrac=1e-10
-angelou_lithium_HG_massfrac=0.0
-angelou_lithium_GB_massfrac=0.0
-angelou_lithium_CHeB_massfrac=0.0
-angelou_lithium_EAGB_massfrac=0.0
-angelou_lithium_TPAGB_massfrac=0.0
+angelou_massfrac=1e-8
+angelou_lithium_LMMS_massfrac=$angelou_massfrac
+angelou_lithium_MS_massfrac=$angelou_massfrac
+angelou_lithium_HG_massfrac=$angelou_massfrac
+angelou_lithium_GB_massfrac=$angelou_massfrac
+angelou_lithium_CHeB_massfrac=$angelou_massfrac
+angelou_lithium_EAGB_massfrac=$angelou_massfrac
+angelou_lithium_TPAGB_massfrac=$angelou_massfrac
############################################################
@@ -1597,6 +1550,7 @@ $HE_IG_EXTRA \
--maximum_nuclear_burning_timestep $MAXIMUM_NUCLEAR_BURNING_DT \
--stellar_structure_algorithm $STELLAR_STRUCTURE_ALGORITHM \
--angelou_lithium_decay_function $angelou_lithium_decay_function \
+--angelou_lithium_LMMS_time $angelou_lithium_LMMS_time \
--angelou_lithium_MS_time $angelou_lithium_MS_time \
--angelou_lithium_HG_time $angelou_lithium_HG_time \
--angelou_lithium_GB_time $angelou_lithium_GB_time \
@@ -1605,12 +1559,14 @@ $HE_IG_EXTRA \
--angelou_lithium_TPAGB_time $angelou_lithium_TPAGB_time \
--angelou_lithium_vrot_trigger $angelou_lithium_vrot_trigger \
--angelou_lithium_vrotfrac_trigger $angelou_lithium_vrotfrac_trigger \
+--angelou_lithium_LMMS_decay_time $angelou_lithium_LMMS_decay_time \
--angelou_lithium_MS_decay_time $angelou_lithium_MS_decay_time \
--angelou_lithium_HG_decay_time $angelou_lithium_HG_decay_time \
--angelou_lithium_GB_decay_time $angelou_lithium_GB_decay_time \
--angelou_lithium_CHeB_decay_time $angelou_lithium_CHeB_decay_time \
--angelou_lithium_EAGB_decay_time $angelou_lithium_EAGB_decay_time \
--angelou_lithium_TPAGB_decay_time $angelou_lithium_TPAGB_decay_time \
+--angelou_lithium_LMMS_massfrac $angelou_lithium_LMMS_massfrac \
--angelou_lithium_MS_massfrac $angelou_lithium_MS_massfrac \
--angelou_lithium_HG_massfrac $angelou_lithium_HG_massfrac \
--angelou_lithium_GB_massfrac $angelou_lithium_GB_massfrac \
--
GitLab