diff --git a/src/MINT/MINT_prototypes.h b/src/MINT/MINT_prototypes.h
index 1e0d6b24ceff910db10f305db533cccb88319fc4..4e0b0a72d38af734a6b44dc21f549be64694a260 100644
--- a/src/MINT/MINT_prototypes.h
+++ b/src/MINT/MINT_prototypes.h
@@ -105,7 +105,9 @@ Stellar_type MINT_stellar_structure_NS(struct stardata_t * RESTRICT const starda
Stellar_type MINT_stellar_structure_BH(struct stardata_t * RESTRICT const stardata,
struct star_t * const prevstar,
struct star_t * const newstar);
-
+Stellar_type MINT_stellar_structure_massless_remnant(struct stardata_t * RESTRICT const stardata MAYBE_UNUSED,
+ struct star_t * const prevstar MAYBE_UNUSED,
+ struct star_t * const newstar);
diff --git a/src/MINT/MINT_stellar_structure.c b/src/MINT/MINT_stellar_structure.c
index 9c7f40fa6c00f57fa85375693289e06fd1aab936..11fc5f1743c58da2b165f3ca0897ff2e49007971 100644
--- a/src/MINT/MINT_stellar_structure.c
+++ b/src/MINT/MINT_stellar_structure.c
@@ -11,6 +11,16 @@ No_empty_translation_unit_warning;
*
* However, oldstar *can* be NULL, in which case
* we simply update the structure in newstar.
+ *
+ * While we are developing MINT, you can choose
+ * whether to call the MINT function before or
+ * after BSE, and whether to use the BSE function
+ * for each stellar type.
+ *
+ * Please see
+ * MINT_stellar_structure_functions.def
+ * for instructions.
+ *
*/
static int call_BSE_stellar_structure(struct stardata_t * const stardata,
struct star_t * const oldstar,
@@ -22,13 +32,16 @@ int MINT_stellar_structure(struct stardata_t * const stardata,
struct star_t * const newstar,
const Caller_id caller_id MAYBE_UNUSED)
{
- const Boolean vb = FALSE;
+ Boolean vb = FALSE;
/*
* Trim shells, e.g. because of mass loss
*/
MINT_trim_shells(stardata,
newstar);
-
+ printf("MINT model %d (%d %d)\n",
+ stardata->model.model_number,
+ oldstar->stellar_type,
+ newstar->stellar_type);
if(oldstar != NULL)
{
/*
@@ -64,8 +77,28 @@ int MINT_stellar_structure(struct stardata_t * const stardata,
newstar->tm);
}
}
+ else if(oldstar->stellar_type == CHeB)
+ {
+ if(oldstar->mint->XHec <= 0.0)
+ {
+ /*
+ * Finished CHeB
+ */
+ static int count = 0;
+ newstar->stellar_type = EAGB;
+ newstar->mint->XHec = 0.0;
+ printf("MINT to EAGB %d\n",count++);
+ if(count>0)
+ {
+ Backtrace;
+ //Flexit;
+ }
+ }
+ }
}
+ if(newstar->stellar_type == EAGB) vb = TRUE;
+
int ret = newstar->stellar_type;
#undef X
#define X(STELLAR_TYPE,FUNCTION,CALL_BSE) FUNCTION,
@@ -78,17 +111,23 @@ int MINT_stellar_structure(struct stardata_t * const stardata,
static const unsigned int call_bse_action[] =
{ MINT_STELLAR_STRUCTURE_FUNCTIONS_LIST };
#undef X
+#define X(ARG) #ARG,
+ static const char * const call_bse_action_strings[] =
+ { MINT_CALL_BSE_ACTIONS };
+#undef X
/*
* Call BSE before MINT?
*/
- if(vb)printf("MINT t=%g st = %d action = %u\n",
+ const unsigned int action = call_bse_action[newstar->stellar_type];
+ if(vb)printf("MINT t=%g st = %d action = %u %s\n",
stardata->model.time,
newstar->stellar_type,
- call_bse_action[newstar->stellar_type]);
+ call_bse_action[newstar->stellar_type],
+ call_bse_action_strings[action]);
- if(call_bse_action[newstar->stellar_type] == MINT_CALL_BSE_BEFORE ||
- call_bse_action[newstar->stellar_type] == MINT_CALL_BSE_BOTH)
+ if(action == MINT_CALL_BSE_BEFORE ||
+ action == MINT_CALL_BSE_BOTH)
{
ret = call_BSE_stellar_structure(stardata,
oldstar,
@@ -109,8 +148,8 @@ int MINT_stellar_structure(struct stardata_t * const stardata,
/*
* Call BSE after MINT?
*/
- if(call_bse_action[newstar->stellar_type] == MINT_CALL_BSE_AFTER ||
- call_bse_action[newstar->stellar_type] == MINT_CALL_BSE_BOTH)
+ if(action == MINT_CALL_BSE_AFTER ||
+ action == MINT_CALL_BSE_BOTH)
{
ret = call_BSE_stellar_structure(stardata,
oldstar,
@@ -118,6 +157,10 @@ int MINT_stellar_structure(struct stardata_t * const stardata,
caller_id);
}
+ if(vb)printf("MINT_stellar_structure return %d (stellar types %d %d)\n",
+ ret,
+ oldstar->stellar_type,
+ newstar->stellar_type);
return ret;
}
diff --git a/src/MINT/MINT_stellar_structure_CHeB.c b/src/MINT/MINT_stellar_structure_CHeB.c
index 81ad9755594dbf1efdcc8138d86d98afdf561db4..3ed4a491840859bf2fb04b0c91ce19e9b5e05e4c 100644
--- a/src/MINT/MINT_stellar_structure_CHeB.c
+++ b/src/MINT/MINT_stellar_structure_CHeB.c
@@ -336,15 +336,6 @@ Stellar_type MINT_stellar_structure_CHeB(struct stardata_t * RESTRICT const star
}
- if(newstar->mint->XHec <= 0.0)
- {
- /*
- * Finished CHeB
- */
- newstar->stellar_type = EAGB;
- newstar->mint->XHec = 0.0;
- }
-
if(vb)
{
printf("MINT CHeB star %d: stellar_type=%d (Z=%g, M=%g, XHec=%g) at t=%g = L=%g R=%g dYHec/dt=%g\n",
diff --git a/src/MINT/MINT_stellar_structure_functions.def b/src/MINT/MINT_stellar_structure_functions.def
index 955100e9fdf73d5b238fc5dba365ae2ce1638541..16789fbe8687f115b913047e3a0ba88511733b3f 100644
--- a/src/MINT/MINT_stellar_structure_functions.def
+++ b/src/MINT/MINT_stellar_structure_functions.def
@@ -1,12 +1,19 @@
#pragma once
#ifndef MINT_STELLAR_STRUCTURE_FUNCTIONS_DEF
#define MINT_STELLAR_STRUCTURE_FUNCTIONS_DEF
+
+#define MINT_CALL_BSE_ACTIONS \
+ X(FALSE) \
+ X(BEFORE) \
+ X(AFTER) \
+ X(BOTH)
+
+#undef X
+#define X(ARG) MINT_CALL_BSE_##ARG,
enum {
- MINT_CALL_BSE_FALSE,
- MINT_CALL_BSE_BEFORE,
- MINT_CALL_BSE_AFTER,
- MINT_CALL_BSE_BOTH
+ MINT_CALL_BSE_ACTIONS
};
+#undef X
/*
* Stellar structure functions called when
@@ -30,14 +37,15 @@ enum {
X(GIANT_BRANCH, &MINT_stellar_structure_GB, BEFORE) \
X(CORE_HELIUM_BURNING, &MINT_stellar_structure_CHeB, FALSE) \
X(EAGB, &MINT_stellar_structure_EAGB, BEFORE) \
- X(TPAGB, &MINT_stellar_structure_TPAGB, BEFORE) \
- X(HeMS, &MINT_stellar_structure_HeMS, BEFORE) \
- X(HeHG, &MINT_stellar_structure_HeHG, BEFORE) \
- X(HeGB, &MINT_stellar_structure_HeGB, BEFORE) \
- X(HeWD, &MINT_stellar_structure_HeWD, BEFORE) \
- X(COWD, &MINT_stellar_structure_COWD, BEFORE) \
- X(ONeWD, &MINT_stellar_structure_ONeWD, BEFORE) \
- X(NS, &MINT_stellar_structure_NS, BEFORE) \
- X(BH, &MINT_stellar_structure_BH, BEFORE)
+ X(TPAGB, &MINT_stellar_structure_TPAGB, BEFORE) \
+ X(HeMS, &MINT_stellar_structure_HeMS, BEFORE) \
+ X(HeHG, &MINT_stellar_structure_HeHG, BEFORE) \
+ X(HeGB, &MINT_stellar_structure_HeGB, BEFORE) \
+ X(HeWD, &MINT_stellar_structure_HeWD, BEFORE) \
+ X(COWD, &MINT_stellar_structure_COWD, BEFORE) \
+ X(ONeWD, &MINT_stellar_structure_ONeWD, BEFORE) \
+ X(NS, &MINT_stellar_structure_NS, BEFORE) \
+ X(BH, &MINT_stellar_structure_BH, BEFORE) \
+ X(MASSLESS_REMNANT, &MINT_stellar_structure_massless_remnant, BEFORE)
#endif // MINT_STELLAR_STRUCTURE_FUNCTIONS_DEF
diff --git a/src/MINT/MINT_stellar_structure_massless_remnant.c b/src/MINT/MINT_stellar_structure_massless_remnant.c
new file mode 100644
index 0000000000000000000000000000000000000000..fd5629dced2d33681e84d79fa1da627fbde79fc5
--- /dev/null
+++ b/src/MINT/MINT_stellar_structure_massless_remnant.c
@@ -0,0 +1,24 @@
+#include "../binary_c.h"
+No_empty_translation_unit_warning;
+
+
+#ifdef MINT
+
+/*
+ * MINT function to calculate stellar structure
+ * of a massless remnant.
+ */
+
+Stellar_type MINT_stellar_structure_massless_remnant(struct stardata_t * RESTRICT const stardata MAYBE_UNUSED,
+ struct star_t * const prevstar MAYBE_UNUSED,
+ struct star_t * const newstar)
+{
+ /*
+ * Do nothing
+ */
+
+
+ return newstar->stellar_type;
+}
+
+#endif // MINT
diff --git a/src/binary_c_debug.h b/src/binary_c_debug.h
index 0d288a0d4481bcb269c085a1448fde980a325d30..5cf5947ecdfd115ddc5d4e11dd1a7a347a1c5256 100644
--- a/src/binary_c_debug.h
+++ b/src/binary_c_debug.h
@@ -78,7 +78,7 @@
* the same, but stardata is defined everywhere as a global
* variable.
*/
-#define Debug_expression (1)
+#define Debug_expression (stardata->model.model_number>=349)
/*
* If you define Debug_stop_expression, and it is at any time TRUE,
@@ -106,7 +106,7 @@
* this is not shown.
*/
-#define Debug_show_expression " "
+#define Debug_show_expression " %d %d ",stardata->model.model_number,stardata->star[0].stellar_type
//#undef Debug_show_expression
/*
diff --git a/src/binary_c_structures.h b/src/binary_c_structures.h
index 6130d5bc995c59c3d6d427323d0eb6eaa2abf3cc..7a09ab6558b3afc982a3f5b26cea05cb87f87ce0 100644
--- a/src/binary_c_structures.h
+++ b/src/binary_c_structures.h
@@ -728,6 +728,8 @@ struct preferences_t {
double vwind_beta;
// wind gas/dust ratio (typically 200-500)
double wind_gas_to_dust_ratio;
+ /* Tout and Pringle 1992 wind multiplier */
+ double Tout_Pringle_1992_multiplier;
/* artificial accretion rates */
double artificial_mass_accretion_rate[NUMBER_OF_STARS];
@@ -775,7 +777,6 @@ struct preferences_t {
double PN_fast_wind_dm_AGB;
double PN_fast_wind_dm_GB;
-
/* Common envelope parameters */
double alpha_ce[MAX_NUMBER_OF_COMENVS];
double lambda_ce[MAX_NUMBER_OF_COMENVS];
diff --git a/src/binary_star_functions/angular_momentum_and_eccentricity_derivatives.c b/src/binary_star_functions/angular_momentum_and_eccentricity_derivatives.c
index 587d68c861ba263042a1bed739284bebba2a69ac..8aeccb009e58500a4d26a056a29b16605bec8103 100644
--- a/src/binary_star_functions/angular_momentum_and_eccentricity_derivatives.c
+++ b/src/binary_star_functions/angular_momentum_and_eccentricity_derivatives.c
@@ -248,8 +248,7 @@ void angular_momentum_and_eccentricity_derivatives(
stellar_wind_angmom(stardata,RLOF_boolean,star->starnum);
/* spin down by magnetic braking */
- Dprint("call mag braking %d\n",star->starnum);
- magnetic_braking(stardata,star);
+ magnetic_braking(stardata,star);
#ifdef OMEGA_CORE
/* couple core and envelope */
diff --git a/src/logging/log_every_timestep.c b/src/logging/log_every_timestep.c
index e4420343b66d11f5aaeae551e24cfc1c213c928a..d3f3a1961a47faec42ddab9250a4e1ed8fb7364d 100644
--- a/src/logging/log_every_timestep.c
+++ b/src/logging/log_every_timestep.c
@@ -3055,4 +3055,24 @@ void log_every_timestep(struct stardata_t * RESTRICT const stardata)
log_orbiting_objects(stardata);
#endif// ORBITING_OBJECTS
+ /*
+ if(stardata->star[0].stellar_type<10)
+ {
+ printf("XXX %d %d %g %g\n",
+ stardata->model.model_number,
+ stardata->star[0].stellar_type,
+ HRdiagram(0));
+ }
+ */
+/*
+ printf("YYY %g %g %g %g %g %g %g\n",
+ stardata->model.time,
+ stardata->star[0].mass,
+ stardata->star[0].radius,
+ stardata->star[0].luminosity,
+ stardata->star[0].omega / stardata->star[0].omega_crit,
+ stardata->star[0].v_eq / 2.0,
+ stardata->star[0].derivative[DERIVATIVE_STELLAR_MASS_WIND_LOSS]
+ );
+*/
}
diff --git a/src/logging/log_orbiting_objects.c b/src/logging/log_orbiting_objects.c
index 4e6511a4abd4794b28bfcb4a3f1f1f5dc5bd2e6e..40530ce6b11c81af5e10099fcd6bfd2f13c60888 100644
--- a/src/logging/log_orbiting_objects.c
+++ b/src/logging/log_orbiting_objects.c
@@ -60,7 +60,7 @@ void log_orbiting_objects(struct stardata_t * RESTRICT const stardata)
Printf("* * * * * * * ");
}
Printf("%d\n",
- o->in_disc == TRUE ? 1 : 0 // 12
+ o->in_disc == TRUE ? 1 : 0 // 14
);
}
}
diff --git a/src/nucsyn/nucsyn_burning_cycles.c b/src/nucsyn/nucsyn_burning_cycles.c
index d3edc5130a5a3ff64c8fa360e12775bf3d7f42ae..cb5bf4b813c7b48e1790081e9f221fdf9f137926 100644
--- a/src/nucsyn/nucsyn_burning_cycles.c
+++ b/src/nucsyn/nucsyn_burning_cycles.c
@@ -7,6 +7,11 @@ No_empty_translation_unit_warning;
//#define TESTING
+/*
+ * Wrapper to call the various nuclear-burning cycles
+ * in binary_c
+ */
+
void nucsyn_burning_cycles(Abundance * RESTRICT const N,
Reaction_rate * RESTRICT const sigmav,
const double dt,
@@ -16,11 +21,17 @@ void nucsyn_burning_cycles(Abundance * RESTRICT const N,
{
#include "nucsyn_solver_functions.h"
+ /*
+ * We have a vb parameter otherwise
+ * the debug output would be too much
+ */
+ const Boolean vb = FALSE;
+
#undef NUCSYN_NORMALIZE_NUCLEONS
- Dprint("nucsyn_burning_cycles...\n");
+ if(vb)Dprint("nucsyn_burning_cycles...\n");
/*
- * If dt = 0, why bother?
+ * Burn only if dt != 0
*/
if(!Is_zero(dt))
{
@@ -76,7 +87,7 @@ void nucsyn_burning_cycles(Abundance * RESTRICT const N,
density) :
TRUE;
- Dprint("Network %u \"%s\" : burn? %s && %s : for time %g\n",
+ if(vb)Dprint("Network %u \"%s\" : burn? %s && %s : for time %g\n",
i,
nuclear_network_strings[i],
Yesno(stardata->preferences->nucsyn_network[i]),
@@ -131,7 +142,8 @@ void nucsyn_burning_cycles(Abundance * RESTRICT const N,
NULL,
stardata->common.dtguess[i],
stardata,
- i
+ i,
+ vb
);
}
#ifdef TESTING
@@ -146,7 +158,8 @@ void nucsyn_burning_cycles(Abundance * RESTRICT const N,
NULL,
stardata->common.dtguess[i],
stardata,
- i
+ i,
+ vb
);
time += time_step;
diff --git a/src/nucsyn/nucsyn_functions.h b/src/nucsyn/nucsyn_functions.h
index 2162a1e5b2fdb8337cfeef74a089dfb478a477c6..6a37342685aa3fe2ca88db420eaf0e7b991f6f5e 100644
--- a/src/nucsyn/nucsyn_functions.h
+++ b/src/nucsyn/nucsyn_functions.h
@@ -28,6 +28,7 @@ typedef double(*nucsyn_solver_func)(double * const,
nucsyn_logfunc,
const double,
struct stardata_t * const,
- const Reaction_network);
+ const Reaction_network,
+ const Boolean);
#endif // NUCSYN_FUNCTIONS_H
diff --git a/src/nucsyn/nucsyn_network_burn.c b/src/nucsyn/nucsyn_network_burn.c
index aafab914a48375e4f8de3407eeb26732b4577655..e195d9747a6c4cda055b91bb50b2b2a944e7c136 100644
--- a/src/nucsyn/nucsyn_network_burn.c
+++ b/src/nucsyn/nucsyn_network_burn.c
@@ -54,7 +54,8 @@ double nucsyn_network_burn(double * RESTRICT const Nin,
nucsyn_logfunc logfunc,
const double dtguess,
struct stardata_t * const RESTRICT stardata,
- const Reaction_network network_id
+ const Reaction_network network_id,
+ const Boolean vb
)
{
/*
@@ -86,7 +87,7 @@ double nucsyn_network_burn(double * RESTRICT const Nin,
{
logfunc(stardata,Nin,t);
}
-// printf("dt in %g, burn for %g\n",dt,maxt-t);
+
for(;(t<maxt2 && dt>0);) // NB dt=small at input, update dealt with below
{
/*
@@ -98,7 +99,6 @@ double nucsyn_network_burn(double * RESTRICT const Nin,
Nin,
sigmav,
Min(dt,dtr));
-// printf("err %g\n",err);
if(err < 1.0)
{
@@ -106,9 +106,6 @@ double nucsyn_network_burn(double * RESTRICT const Nin,
success_count++;
t += dt;
-// printf("ok err = %g -> return_dt = Min(%g,%g) = %g\n",
- // err,return_dt,dt,Min(return_dt,dt));
-
/* choose a timestep for next time which is the minimum we require */
return_dt = Min(return_dt,dt);
@@ -121,13 +118,8 @@ double nucsyn_network_burn(double * RESTRICT const Nin,
}
/* adjust timestep */
-
- /* standard, linear raise? */
- //dt*=TIMESTEP_RAISE_MULTIPLIER;
-
- /* intelligent algorithm based on the error */
- dt *= Min(MAX_RAISE,INTELLIGENT_RAISE_MULTIPLIER * (-log10(err)));
-// printf("newdt %g\n",dt);
+ dt *= Min(MAX_RAISE,
+ INTELLIGENT_RAISE_MULTIPLIER * (-log10(err)));
}
else
{
@@ -137,9 +129,6 @@ double nucsyn_network_burn(double * RESTRICT const Nin,
failure_count++;
/* failure : reduce dt */
- //dt*=TIMESTEP_LOWER_MULTIPLIER;
-
- /* or make a more intelligient guess */
dt *= Max(MIN_LOWER,INTELLIGENT_LOWER_MULTIPLIER / err);
}
@@ -163,9 +152,7 @@ double nucsyn_network_burn(double * RESTRICT const Nin,
dtr = Min(dt,
maxt2-t);
-// printf("dt %g : %g\n",dt,Nin[XH1]);
-
- Dprint("dt=%12.12e -> err=%g\n",dt,err);
+ if(vb)Dprint("dt=%12.12e -> err=%g\n",dt,err);
}
@@ -183,7 +170,6 @@ double nucsyn_network_burn(double * RESTRICT const Nin,
(unsigned int)NUCSYN_NETWORK_BURN_MAX_FAILURE_COUNT);
}
-// printf("return dt %g\n",return_dt);
return return_dt;
}
diff --git a/src/nucsyn/nucsyn_network_burn_cvode.c b/src/nucsyn/nucsyn_network_burn_cvode.c
index 01cdac0c9cedfc37ca1ef7cb01d43fa2b21c1fbf..6010585fd28ac6851a5b3bd5ca4b00277c1af52c 100644
--- a/src/nucsyn/nucsyn_network_burn_cvode.c
+++ b/src/nucsyn/nucsyn_network_burn_cvode.c
@@ -11,7 +11,8 @@ double nucsyn_network_burn_cvode(double * RESTRICT const Nin MAYBE_UNUSED,
nucsyn_logfunc logfunc MAYBE_UNUSED,
const double dtguess MAYBE_UNUSED,
struct stardata_t * const RESTRICT stardata MAYBE_UNUSED,
- const Reaction_network network_id MAYBE_UNUSED
+ const Reaction_network network_id MAYBE_UNUSED,
+ const Boolean vb
)
{
diff --git a/src/nucsyn/nucsyn_network_burn_lsoda.c b/src/nucsyn/nucsyn_network_burn_lsoda.c
index 805ead50d3d1ca864fc5d1ca9b403d53991ef623..5a40490c2aac366113f2e937d3695a23d842f41a 100644
--- a/src/nucsyn/nucsyn_network_burn_lsoda.c
+++ b/src/nucsyn/nucsyn_network_burn_lsoda.c
@@ -9,7 +9,8 @@ double nucsyn_network_burn_lsoda(double * RESTRICT const Nin,
nucsyn_logfunc logfunc MAYBE_UNUSED,
const double dtguess MAYBE_UNUSED,
struct stardata_t * const RESTRICT stardata,
- const Reaction_network network_id MAYBE_UNUSED
+ const Reaction_network network_id MAYBE_UNUSED,
+ const Boolean vb MAYBE_UNUSED
)
{
diff --git a/src/nucsyn/nucsyn_prototypes.h b/src/nucsyn/nucsyn_prototypes.h
index bee1dbb1e931e96c8609f514e9cf23a78a09d3ff..5943ed4e49e64fbc2c8305d576cabdf1eeca8bf7 100644
--- a/src/nucsyn/nucsyn_prototypes.h
+++ b/src/nucsyn/nucsyn_prototypes.h
@@ -419,7 +419,8 @@ double nucsyn_network_burn(double * RESTRICT const Nin,
nucsyn_logfunc logfunc,
const double dtguess,
struct stardata_t * const RESTRICT stardata,
- const Reaction_network network_id);
+ const Reaction_network network_id,
+ const Boolean vb);
double nucsyn_network_burn_lsoda(double * RESTRICT const Nin,
@@ -429,7 +430,8 @@ double nucsyn_network_burn_lsoda(double * RESTRICT const Nin,
nucsyn_logfunc logfunc,
const double dtguess,
struct stardata_t * const RESTRICT stardata,
- const Reaction_network network_id);
+ const Reaction_network network_id,
+ const Boolean vb);
#ifdef __HAVE_LIBSUNDIALS_CVODE__
double nucsyn_network_burn_cvode(double * RESTRICT const Nin,
@@ -439,7 +441,8 @@ double nucsyn_network_burn_cvode(double * RESTRICT const Nin,
nucsyn_logfunc logfunc MAYBE_UNUSED,
const double dtguess MAYBE_UNUSED,
struct stardata_t * const RESTRICT stardata,
- const Reaction_network network_id MAYBE_UNUSED);
+ const Reaction_network network_id MAYBE_UNUSED,
+ const Boolean vb);
#endif // __HAVE_LIBSUNDIALS_CVODE__
/* nuclear burning functions */
diff --git a/src/setup/cmd_line_args_list.h b/src/setup/cmd_line_args_list.h
index 9260dfac954d945e6c8084a4fccaae86912cb6db..069b9d289ce4d68896b21f5eef3e8eb8ecd4012e 100644
--- a/src/setup/cmd_line_args_list.h
+++ b/src/setup/cmd_line_args_list.h
@@ -1282,6 +1282,15 @@ BATCHMODE_ARGS
Var(stardata->preferences->tpagb_reimers_eta),
1.0
CMD_LINE_ARG_T_OPTIONAL_ARGS },
+{
+ ARG_SECTION_STARS,
+ "Tout_Pringle_1992_multiplier",
+ "Multiplier for the Tout & Pringle (1992) magnetic wind. (0.0)",
+ ARG_FLOAT ,
+ WTTS_USE_DEFAULT,
+ Var(stardata->preferences->Tout_Pringle_1992_multiplier),
+ 1.0
+ CMD_LINE_ARG_T_OPTIONAL_ARGS },
{
ARG_SECTION_STARS,
diff --git a/src/stellar_structure/stellar_structure_AGB.c b/src/stellar_structure/stellar_structure_AGB.c
index 8fecc6ce1083ac8b3adb45d8993c41ace60a28ba..caf4e338e2865148cfca5a4086d97d8a9ad973f3 100644
--- a/src/stellar_structure/stellar_structure_AGB.c
+++ b/src/stellar_structure/stellar_structure_AGB.c
@@ -78,12 +78,16 @@ Stellar_type stellar_structure_AGB(struct star_t * const oldstar,
newstar->core_mass[CORE_He]
);
+ /*
+ * Obtain Chandrasekhar mass
+ */
+ const double MCh = chandrasekhar_mass(stardata,newstar);
/*
* Equivalent to Hurley et al. 2000's Eq. (75)
*/
const double mcSN = Max(0.7730 * mcbagb - 0.350,
- chandrasekhar_mass(stardata,newstar));
+ MCh);
/*
@@ -134,8 +138,8 @@ Stellar_type stellar_structure_AGB(struct star_t * const oldstar,
* luminosity, radius and age.
*/
if(newstar->age < newstar->timescales[T_TPAGB_FIRST_PULSE]
- &&
- newstar->stellar_type!=TPAGB)
+ &&
+ newstar->stellar_type!=TPAGB)
{
/* EAGB star : updates mc_CO */
convert_to_helium_star=
@@ -160,6 +164,7 @@ Stellar_type stellar_structure_AGB(struct star_t * const oldstar,
*/
else if((newstar->stellar_type == TPAGB ||
(Is_not_zero(newstar->core_mass[CORE_CO]) &&
+ newstar->core_mass[CORE_He] < MCh &&
newstar->core_mass[CORE_CO] < maximum_degenerate_CO_mcbagb)))
{
/* TPAGB star */
@@ -275,11 +280,11 @@ Stellar_type stellar_structure_AGB(struct star_t * const oldstar,
* exceeds either the total mass
* or the maximum allowed core mass.
*/
- Dprint("COMPARE st=%d mc_CO=%12.12e vs mcmax=%12.12e (Mch=%g) r=%12.12e\n",
+ Dprint("COMPARE st=%d mc_CO=%12.12e vs mcmax=%12.12e (MCh=%g) r=%12.12e\n",
newstar->stellar_type,
newstar->core_mass[CORE_CO],
mcmax,
- stardata->preferences->chandrasekhar_mass,
+ MCh,
newstar->radius);
Dprint("MC_CO=%g mcmax=%g mdot=%g dt=%g\n",
@@ -288,13 +293,13 @@ Stellar_type stellar_structure_AGB(struct star_t * const oldstar,
newstar->derivative[DERIVATIVE_STELLAR_MASS_WIND_LOSS],
stardata->model.dt);
- Dprint("Mc1TP0 st=%d mc_CO=%12.12e (McCO start AGB = %g) vs mcmax=%12.12e (Mch=%g McSN=%g) r=%12.12e\n",
+ Dprint("Mc1TP0 st=%d mc_CO=%12.12e (McCO start AGB = %g) vs mcmax=%12.12e (MCh=%g McSN=%g) r=%12.12e\n",
newstar->stellar_type,
newstar->core_mass[CORE_CO],
mc_CO_start_AGB,
mcmax,
mcSN,
- stardata->preferences->chandrasekhar_mass,
+ MCh,
newstar->radius);
const double envelope_mass = newstar->mass - newstar->core_mass[CORE_He];
@@ -508,7 +513,7 @@ Stellar_type stellar_structure_AGB(struct star_t * const oldstar,
/*
* Star is not massive enough to start C burning,
* which means mcbagb < 1.6 Msun,
- * yet the star reaches Mch with a hydrogen envelope.
+ * yet the star reaches MCh with a hydrogen envelope.
*
* The core is degenerate carbon and oxygen,
* so no remnant is left after the SN (type "IIa" SN)
diff --git a/src/stellar_structure/stellar_structure_BSE.c b/src/stellar_structure/stellar_structure_BSE.c
index 7227d0a419f2b23a75ff3b73746d0655390b2c4b..5573f5f13c3aaab96b6640dd0d35c8e572569e00 100644
--- a/src/stellar_structure/stellar_structure_BSE.c
+++ b/src/stellar_structure/stellar_structure_BSE.c
@@ -12,8 +12,7 @@
No_empty_translation_unit_warning;
#include "stellar_structure_debug.h"
-//#undef DEBUG
-//#define DEBUG 1
+
#ifdef BSE
int stellar_structure_BSE(struct stardata_t * const stardata,
const Caller_id caller_id,
diff --git a/src/stellar_structure/stellar_structure_EAGB.c b/src/stellar_structure/stellar_structure_EAGB.c
index 3e826257fa20f47f2f725247f919f9b5b533979d..55e2794c4383b363e5cc1dfb44126937fcdc05e8 100644
--- a/src/stellar_structure/stellar_structure_EAGB.c
+++ b/src/stellar_structure/stellar_structure_EAGB.c
@@ -11,10 +11,10 @@ No_empty_translation_unit_warning;
#define D (newstar->GB[GB_D])
#define A (newstar->GB[GB_A_HE])
-Boolean stellar_structure_EAGB(struct star_t * oldstar,
- struct star_t * newstar,
+Boolean stellar_structure_EAGB(struct star_t * const oldstar,
+ struct star_t * const newstar,
const double mcbagb,
- struct stardata_t * stardata,
+ struct stardata_t * const stardata,
const Caller_id caller_id MAYBE_UNUSED,
const Boolean force_to_explosion)
{
diff --git a/src/stellar_structure/stellar_structure_prototypes.h b/src/stellar_structure/stellar_structure_prototypes.h
index 8e7d5b6265bba17d691b567ff9f33f0e78a7af70..53167625ab4fe31dfccfcfdbb120133a12f69133 100644
--- a/src/stellar_structure/stellar_structure_prototypes.h
+++ b/src/stellar_structure/stellar_structure_prototypes.h
@@ -46,7 +46,7 @@ Stellar_type stellar_structure_AGB(struct star_t * const oldstar,
struct stardata_t * RESTRICT const stardata,
const Caller_id caller_id);
-Boolean stellar_structure_EAGB(struct star_t * oldstar,
+Boolean stellar_structure_EAGB(struct star_t * const oldstar,
struct star_t * RESTRICT const newstar,
const double mcbagb,
struct stardata_t * RESTRICT const stardata,
diff --git a/src/wind/wind_Tout_Pringle_1992.c b/src/wind/wind_Tout_Pringle_1992.c
new file mode 100644
index 0000000000000000000000000000000000000000..a8f5e895e9d28b653ee6889d0d87e8c8c65d525c
--- /dev/null
+++ b/src/wind/wind_Tout_Pringle_1992.c
@@ -0,0 +1,63 @@
+#include "../binary_c.h"
+No_empty_translation_unit_warning;
+
+#include "wind.h"
+
+/*
+ * Magnetic dynamo-driven wind from Tout and Pringle (1992)
+ */
+
+double Constant_function_if_no_debug wind_Tout_Pringle_1992(WIND_PROTOTYPE_ARGS)
+{
+ double mdot;
+ if(Is_not_zero(stardata->preferences->Tout_Pringle_1992_multiplier))
+ {
+ /*
+ * Tout's R* should be Renv, the convective envelope
+ */
+ const double Rstar = star->renv;
+
+ if(Is_not_zero(Rstar))
+ {
+ /*
+ * Convective turnover length (Hurley PhD Eq. B.21)
+ */
+ const double l = Rstar / 3.0;
+
+ /*
+ * eta = 3.0 * radius / convective turnover length
+ * (end section 3)
+ */
+ const double eta = 3.0 * (Rstar / l);
+
+ /*
+ * f is defined just after Eq. 5.6
+ */
+ const double f =
+ Is_zero(star->omega_crit) ? 0.0 :
+ (star->omega / star->omega_crit);
+
+ /*
+ * Eq. 5.12
+ */
+ mdot =
+ 1.4e-6 *
+ f *
+ pow(star->luminosity/4.0, 2.0/3.0) *
+ pow(eta/30.0, -2.0/3.0) *
+ pow(star->mass, -1.0/6.0) *
+ pow(Rstar/3.0, 1.0/6.0)
+ *
+ stardata->preferences->Tout_Pringle_1992_multiplier;
+ }
+ else
+ {
+ mdot = 0.0;
+ }
+ }
+ else
+ {
+ mdot = 0.0;
+ }
+ return mdot;
+}
diff --git a/src/wind/wind_binary_c_2020.c b/src/wind/wind_binary_c_2020.c
index 8eeef6b80c2b50c8a33746b32b6d93fd90a6d70a..2337a50644425d6a965e3f9e88aa55f09075aab8 100644
--- a/src/wind/wind_binary_c_2020.c
+++ b/src/wind/wind_binary_c_2020.c
@@ -16,8 +16,8 @@ double wind_binary_c_2020(WIND_PROTOTYPE_ARGS)
mdot = wind_WD(WIND_CALL_ARGS);
}
else if(stardata->preferences->wind_mass_loss == WIND_ALGORITHM_BINARY_C_2020 &&
- NAKED_HELIUM_STAR(star->stellar_type) &&
- log10(star->luminosity) > 4.85)
+ NAKED_HELIUM_STAR(star->stellar_type) &&
+ log10(star->luminosity) > 4.85)
{
/*
* Sander and Vink 2020
diff --git a/src/wind/wind_mass_loss_rate.c b/src/wind/wind_mass_loss_rate.c
index 4c8f3b2a4d51f38c6961d8e8706488389005f99b..d33acad3ea2510718dfd868d4cca3b3b71f413b1 100644
--- a/src/wind/wind_mass_loss_rate.c
+++ b/src/wind/wind_mass_loss_rate.c
@@ -58,20 +58,18 @@ double wind_mass_loss_rate(WIND_PROTOTYPE_ARGS)
);
}
+
+
+ /*
+ * Apply Tout and Pringle (1992) magnetic wind
+ */
+ const double mdot_TP92 = wind_Tout_Pringle_1992(WIND_CALL_ARGS);
+ mdot_wind = Max(mdot_wind,mdot_TP92);
+
+
/*
* modulate by wind_multiplier_%d where %d is the stellar type
*/
-/*
- if(star->starnum==0 && star->stellar_type==TPAGB)
- printf("model %d wind %g tscl %g : M %g dt %g st %d\n",
- stardata->model.model_number,
- mdot_wind,
- (star->mass - Outermost_core_mass(star)) / mdot_wind,
- star->mass,
- stardata->model.dt,
- star->stellar_type
- );
-*/
mdot_wind *= stardata->preferences->wind_multiplier[star->stellar_type];
return mdot_wind;
diff --git a/src/wind/wind_prototypes.h b/src/wind/wind_prototypes.h
index 101ac977b29154517bad3a112524e052f27da01b..c878a9c0dd835ad1d694419b8b2ab0a5fb4cd91e 100644
--- a/src/wind/wind_prototypes.h
+++ b/src/wind/wind_prototypes.h
@@ -34,5 +34,6 @@ double Pure_function wind_mixing_factor(const struct stardata_t * const stardata
double wind_velocity(struct stardata_t * const stardata,
struct star_t * const star);
double Mira_wind_velocity(const double period);
+double Constant_function_if_no_debug wind_Tout_Pringle_1992(WIND_PROTOTYPE_ARGS);
#endif // WIND_PROTOTYPES_H
diff --git a/src/zfuncs/mcgbtf.c b/src/zfuncs/mcgbtf.c
index 86a7d5eeb55d79850101997bbe1588dc30bda259..d54a2aedf8b90e17d235630c3f9ce358dadf13af 100644
--- a/src/zfuncs/mcgbtf.c
+++ b/src/zfuncs/mcgbtf.c
@@ -21,29 +21,11 @@ double Pure_function mcgbtf(const double t,
{
y=GB[GB_p]-1.0;
result = pow(y*A*GB[GB_D]*(tinf1 - t),-1.0/y);
- /*printf("mcx y1 %g : y=%g GB_D=%g tinf1=%g t=%g -> result %g\n",
- result,
- y,
- GB[GB_D],
- tinf1,
- t,
- result
- );
- */
}
else
{
y=GB[GB_q]-1.0;
result = pow(y*A*GB[GB_B]*(tinf2 - t),-1.0/y);
- /*printf("mcx y2 %g : y=%g GB_D=%g tinf2=%30.20g t=%30.20g -> result %g\n",
- result,
- y,
- GB[GB_B],
- tinf2,
- t,
- result
- );
- */
}
return result;
diff --git a/tbse b/tbse
index ee67f9ea2b7716b23f3ae916e8a7cbfd898a883a..7bf358db9e842bea7b55ce27507639f52d618b07 100755
--- a/tbse
+++ b/tbse
@@ -1017,6 +1017,9 @@ VWIND_MULTIPLIER=1.0
# 200 is typical in the Milky Way, 500 in LMC (see Goldman et al. 2017)
WIND_GAS_TO_DUST_RATIO=200.0
+# Tout and Pringle (1992) magnetic wind multiplier (0.0)
+TOUT_PRINGLE_1992_MULTIPLIER=1.0
+
# Wind during the first Giant Branch:
#
# GB_WIND_REIMERS = 0 : Reimers (original, choose eta~0.5),
@@ -2240,6 +2243,7 @@ DEFAULT_ARGS=" \
--Teff_postAGB_min $TEFF_POSTAGB_MIN \
--wind_gas_to_dust_ratio $WIND_GAS_TO_DUST_RATIO \
--gb_reimers_eta $GB_REIMERS_ETA \
+--Tout_Pringle_1992_multiplier $TOUT_PRINGLE_1992_MULTIPLIER \
--probability $PROBABILITY \
--accretion_limit_eddington_steady_multiplier $ACCRETION_LIMIT_EDDINGTON_STEADY_MULTIPLIER \
--accretion_limit_eddington_WD_to_remnant_multiplier $ACCRETION_LIMIT_EDDINGTON_WD_TO_REMNANT_MULTIPLIER \