From bf6bad79af9b457cd7a906e6277d5d87be07c3b7 Mon Sep 17 00:00:00 2001
From: Robert Izzard <r.izzard@surrey.ac.uk>
Date: Tue, 13 Jul 2021 18:20:41 +0100
Subject: [PATCH] add post-AGB wind options
add function to detect when we're on the post-AGB
---
src/binary_c_structures.h | 8 +-
src/disc/disc_logging.c | 19 +++-
.../minimum_stable_orbit_separation.c | 10 +-
src/setup/cmd_line_args_list.h | 29 ++++++
src/setup/set_cmd_line_macro_pairs.c | 3 +
src/setup/set_default_preferences.c | 5 +
src/setup/version.c | 8 +-
.../single_star_functions_prototypes.h | 2 +
src/single_star_functions/star_is_postAGB.c | 30 ++++++
src/stellar_structure/stellar_structure_BSE.c | 1 +
src/wind/wind_AGB.c | 9 +-
src/wind/wind_GB.c | 18 +++-
src/wind/wind_Hurley2002.c | 24 ++++-
src/wind/wind_Schneider.c | 14 +++
src/wind/wind_WD.c | 15 +++
src/wind/wind_binary_c_2020.c | 8 +-
src/wind/wind_mass_loss_rate.c | 22 ++---
src/wind/wind_postAGB.c | 93 +++++++++++++++++++
src/wind/wind_prescriptions.def | 5 +
src/wind/wind_prescriptions.h | 6 ++
src/wind/wind_prototypes.h | 3 +
tbse | 11 +++
22 files changed, 314 insertions(+), 29 deletions(-)
create mode 100644 src/single_star_functions/star_is_postAGB.c
create mode 100644 src/wind/wind_WD.c
create mode 100644 src/wind/wind_postAGB.c
diff --git a/src/binary_c_structures.h b/src/binary_c_structures.h
index 03947bb09..7955add9b 100644
--- a/src/binary_c_structures.h
+++ b/src/binary_c_structures.h
@@ -749,6 +749,12 @@ struct preferences_t {
double wind_LBV_luminosity_lower_limit;
double wind_Nieuwenhuijzen_luminosity_lower_limit;
+ /*
+ * Temperature range of post-(A)GB stars
+ */
+ double Teff_postAGB_min;
+ double Teff_postAGB_max;
+
/*
* post-common envelope transition time in years (1e2 - 1e4)
*/
@@ -1075,7 +1081,7 @@ struct preferences_t {
#endif//EQUATION_OF_STATE_ALGORITHMS
/* different wind choices */
- int gbwind,eagbwind,tpagbwind,wr_wind;
+ int gbwind,eagbwind,tpagbwind,wr_wind,postagbwind;
/* algorithm to handle overspin */
int overspin_algorithm;
diff --git a/src/disc/disc_logging.c b/src/disc/disc_logging.c
index 0fcaefd7a..92b08a672 100644
--- a/src/disc/disc_logging.c
+++ b/src/disc/disc_logging.c
@@ -503,7 +503,7 @@ void disc_logging(struct stardata_t * const stardata,
stardata->tmpstore->disc_logfilecount++;
binary_c_fprintf(
Disc_log_filepointer,
- "\"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\"\n",
+ "\"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\" \"%s\"\n",
"t",
"st1",
"st2",
@@ -603,7 +603,11 @@ void disc_logging(struct stardata_t * const stardata,
"","",
#endif
"Mc1dot",
- "Mc2dot"
+ "Mc2dot",
+ "M1dotwind",
+ "M2dotwind",
+ "M1dotdisc",
+ "M2dotdisc"
);
binary_c_fprintf(Disc_log_filepointer,
@@ -685,7 +689,7 @@ void disc_logging(struct stardata_t * const stardata,
binary_c_filter_fprintf(
Disc_log_filepointer,
FILE_FILTER_STRIP_ARGUMENTS,
- "t=%.20e st1=%d st2=%d nzones=%u zone0=%d zone1=%d zone2=%d Rin=%20.12g Rout=%20.12g RMhalf=%20.12g RJhalf=%20.12g a=%20.12g a1=%20.12g a2=%20.12g aL1=%20.12g aL2=%20.12g aL3=%20.12g Revapin=%20.12g Revapout=%20.12g RISM=%20.12g R0=%20.12g R1=%20.12g RL0=%20.12g RL1=%20.12g RZ0=%20.12g RZ1=%20.12g torqueF=%20.12g Ldisc=%20.12g Idisc=%20.12g Fdisc=%20.12g Jdisc=%20.12g Jbinary=%20.12g hdisc=%20.12g Mdisc=%20.12g Mbinary=%20.12g Mc0=%20.12g Mc1=%20.12g Menv1=%20.12g Menv2=%20.12g Mdotvisc=%20.12g Mdoscglobal=%20.12g MdotL2=%20.12g MdotISM=%20.12g MdotFUV=%20.12g MdotX=%20.12g Mdotin=%20.12g Mdotout=%20.12g Jdotvisc=%20.12g Jdotglobal=%20.12g JdotL2=%20.12g JdotISM=%20.12g JdotFUV=%20.12g JdotX=%20.12g Jdotin=%20.12g Jdotout=%20.12g Jdottorque=%20.12g hwind=%20.12g e=%20.12g edot=%20.12g Tin=%20.12g Tout=%20.12g densityin=%20.12g densityout=%20.12g sigmain=%20.12g sigmaout=%20.12g Pgravin=%20.12g Pgravout=%20.12g PISM=%20.12g HRin=%20.12g HRout=%20.12g Teffmaxstar=%20.12g Tratio=%20.12g Lstar=%20.12g LstarX=%20.12g Ldisc=%20.12g tviscin=%20.12g tviscout=%20.12g PEdisc=%20.12g KEdisc=%20.12g sigma0=%20.12g Tvisc0=%20.12g epsM=%20.12g epsJ=%20.12g epsF=%20.12g epstorqueF=%20.12g solver=%d guess=%d fRing=%20.12g ToomreQ=%20.12g Terrorpc=%20.12g Mdoterrorpc=%20.12g Mc0dot=%20.12g Mc1dot=%20.12g\n",
+ "t=%.20e st1=%d st2=%d nzones=%u zone0=%d zone1=%d zone2=%d Rin=%20.12g Rout=%20.12g RMhalf=%20.12g RJhalf=%20.12g a=%20.12g a1=%20.12g a2=%20.12g aL1=%20.12g aL2=%20.12g aL3=%20.12g Revapin=%20.12g Revapout=%20.12g RISM=%20.12g R0=%20.12g R1=%20.12g RL0=%20.12g RL1=%20.12g RZ0=%20.12g RZ1=%20.12g torqueF=%20.12g Ldisc=%20.12g Idisc=%20.12g Fdisc=%20.12g Jdisc=%20.12g Jbinary=%20.12g hdisc=%20.12g Mdisc=%20.12g Mbinary=%20.12g Mc0=%20.12g Mc1=%20.12g Menv1=%20.12g Menv2=%20.12g Mdotvisc=%20.12g Mdoscglobal=%20.12g MdotL2=%20.12g MdotISM=%20.12g MdotFUV=%20.12g MdotX=%20.12g Mdotin=%20.12g Mdotout=%20.12g Jdotvisc=%20.12g Jdotglobal=%20.12g JdotL2=%20.12g JdotISM=%20.12g JdotFUV=%20.12g JdotX=%20.12g Jdotin=%20.12g Jdotout=%20.12g Jdottorque=%20.12g hwind=%20.12g e=%20.12g edot=%20.12g Tin=%20.12g Tout=%20.12g densityin=%20.12g densityout=%20.12g sigmain=%20.12g sigmaout=%20.12g Pgravin=%20.12g Pgravout=%20.12g PISM=%20.12g HRin=%20.12g HRout=%20.12g Teffmaxstar=%20.12g Tratio=%20.12g Lstar=%20.12g LstarX=%20.12g Ldisc=%20.12g tviscin=%20.12g tviscout=%20.12g PEdisc=%20.12g KEdisc=%20.12g sigma0=%20.12g Tvisc0=%20.12g epsM=%20.12g epsJ=%20.12g epsF=%20.12g epstorqueF=%20.12g solver=%d guess=%d fRing=%20.12g ToomreQ=%20.12g Terrorpc=%20.12g Mdoterrorpc=%20.12g Mc0dot=%20.12g Mc1dot=%20.12g M1dotwind=%20.12g M2dotwind=%20.12g M1dotdisc=%20.12g M2dotdisc=%20.12g\n",
disc->lifetime/YEAR_LENGTH_IN_SECONDS,
// stellar types
@@ -830,7 +834,14 @@ void disc_logging(struct stardata_t * const stardata,
/* core growth rates */
stardata->star[0].derivative[DERIVATIVE_STELLAR_CO_CORE_MASS],
- stardata->star[1].derivative[DERIVATIVE_STELLAR_CO_CORE_MASS]
+ stardata->star[1].derivative[DERIVATIVE_STELLAR_CO_CORE_MASS],
+ /* stellar wind rates (fabs: >0) */
+ fabs(stardata->star[0].derivative[DERIVATIVE_STELLAR_MASS_WIND_LOSS]),
+ fabs(stardata->star[1].derivative[DERIVATIVE_STELLAR_MASS_WIND_LOSS]),
+ /* accretion rates */
+ disc->lifetime > 0.0 ? stardata->star[0].derivative[DERIVATIVE_STELLAR_MASS_CBDISC_GAIN] : 0.0,
+ stardata->star[1].derivative[DERIVATIVE_STELLAR_MASS_CBDISC_GAIN]
+
);
binary_c_fflush(Disc_log_filepointer);
Safe_free(RZcross);
diff --git a/src/orbiting_objects/minimum_stable_orbit_separation.c b/src/orbiting_objects/minimum_stable_orbit_separation.c
index 4c622c313..bc8f13ebb 100644
--- a/src/orbiting_objects/minimum_stable_orbit_separation.c
+++ b/src/orbiting_objects/minimum_stable_orbit_separation.c
@@ -9,15 +9,21 @@ double minimum_stable_orbit_separation(struct stardata_t * RESTRICT const starda
const double semi_major_axis)
{
/*
- * Aarseth and Mardling minimum stable semi-major axis
+ * Mardling and Aarseth (2001)
+ * MNRAS 321, 398
+ *
+ * Minimum stable semi-major axis
+ * Their Eq. 90
*/
return
semi_major_axis *
(
2.8/(1.0 - o->orbit.eccentricity) *
+
(1.0 - 0.3 * o->orbit.inclination / PI) *
+
pow( (1.0 + o->mass / central_object_mass) *
- (1.0 + o->orbit.inclination) /
+ (1.0 + o->orbit.eccentricity) /
sqrt(1.0 - o->orbit.eccentricity), 0.4)
);
}
diff --git a/src/setup/cmd_line_args_list.h b/src/setup/cmd_line_args_list.h
index 563659d69..f83811df8 100644
--- a/src/setup/cmd_line_args_list.h
+++ b/src/setup/cmd_line_args_list.h
@@ -1079,6 +1079,35 @@ BATCHMODE_ARGS
Var(stardata->preferences->gbwind),
1.0
CMD_LINE_ARG_T_OPTIONAL_ARGS },
+{
+ ARG_SECTION_STARS,
+ "postagbwind",
+ "Apply special post-(A)GB prescription. Default is POSTAGB_WIND_USE_GIANT which means we just use whatever is prescribed on the giant branch. Other options include: POSTAGB_WIND_NONE = 1 (no wind on the post (A)GB), POSTAGB_WIND_KRTICKA2020 = 2 which uses Krticka, Kubát and Krticková (2020, A&A 635, A173).",
+ ARG_INTEGER ,
+ WTTS_USE_DEFAULT,
+ Var(stardata->preferences->postagbwind),
+ 1.0
+ CMD_LINE_ARG_T_OPTIONAL_ARGS },
+{
+ ARG_SECTION_STARS,
+ "Teff_postAGB_min",
+ "The minimum temperature for which we apply post-(A)GB winds. See also Teff_postAGB_max. (6000 K)",
+ ARG_FLOAT ,
+ WTTS_USE_DEFAULT,
+ Var(stardata->preferences->Teff_postAGB_min),
+ 1.0
+ CMD_LINE_ARG_T_OPTIONAL_ARGS },
+{
+ ARG_SECTION_STARS,
+ "Teff_postAGB_max",
+ "The maximum temperature for which we apply post-(A)GB winds. See also Teff_postAGB_min. (120000 K)",
+ ARG_FLOAT ,
+ WTTS_USE_DEFAULT,
+ Var(stardata->preferences->Teff_postAGB_max),
+ 1.0
+ CMD_LINE_ARG_T_OPTIONAL_ARGS },
+
+
{
ARG_SECTION_STARS,
"mattsson_Orich_tpagbwind",
diff --git a/src/setup/set_cmd_line_macro_pairs.c b/src/setup/set_cmd_line_macro_pairs.c
index 7a06d717c..31ba0404b 100644
--- a/src/setup/set_cmd_line_macro_pairs.c
+++ b/src/setup/set_cmd_line_macro_pairs.c
@@ -181,6 +181,9 @@ void set_cmd_line_macro_pairs(struct stardata_t * RESTRICT const stardata,
#define X(CODE) Add_value_macro_pair("tpagbwind",TPAGB_WIND_##CODE);
TPAGB_WIND_PRESCRIPTIONS_LIST;
#undef X
+#define X(CODE) Add_value_macro_pair("postagbwind",POSTAGB_WIND_##CODE);
+ POSTAGB_WIND_PRESCRIPTIONS_LIST;
+#undef X
#define X(CODE) Add_value_macro_pair("wr_wind",WR_WIND_##CODE);
WR_WIND_PRESCRIPTIONS_LIST;
#undef X
diff --git a/src/setup/set_default_preferences.c b/src/setup/set_default_preferences.c
index 0d2396cde..54b36dad8 100644
--- a/src/setup/set_default_preferences.c
+++ b/src/setup/set_default_preferences.c
@@ -622,4 +622,9 @@ preferences->WD_accretion_rate_novae_upper_limit_other_donor = ACCRETION_RATE_NO
preferences->orbiting_objects_wind_accretion_multiplier = BONDI_HOYLE_ACCRETION_FACTOR_DEFAULT;
preferences->orbiting_objects_close_pc_threshold = 1.0;
#endif
+
+ preferences->postagbwind;
+ preferences->Teff_postAGB_min = 6000.0;
+ preferences->Teff_postAGB_max = 120000.0;
+
}
diff --git a/src/setup/version.c b/src/setup/version.c
index fce970d52..8bce3226c 100644
--- a/src/setup/version.c
+++ b/src/setup/version.c
@@ -2067,7 +2067,13 @@ void version(struct stardata_t * RESTRICT const stardata)
Show_float_macro(POST_CE_ENVELOPE_DM_EAGB);
Show_float_macro(POST_CE_ENVELOPE_DM_TPAGB);
-
+#undef X
+#define X(TYPE) printf("Wind MDOT_%s = %d\n",#TYPE,_i++);
+ {
+ int _i = 0;
+ WIND_MDOTS_LIST;
+ }
+#undef X
Macrotest2(AGB_ROTATIONAL_WIND_ENHANCEMENT_DS06,"(Experiental rotational wind enhancement of Dijkstra & Speck 2006)");
Macrotest2(AGB_ROTATIONAL_WIND_ENHANCEMENT_MOE,"(Experiental rotational wind enhancement of Moe)");
Macrotest2(WRLOF_MASS_TRANSFER,"(Carlo Abate's Wind-RLOF)");
diff --git a/src/single_star_functions/single_star_functions_prototypes.h b/src/single_star_functions/single_star_functions_prototypes.h
index c4dae2472..aa759510b 100644
--- a/src/single_star_functions/single_star_functions_prototypes.h
+++ b/src/single_star_functions/single_star_functions_prototypes.h
@@ -89,6 +89,8 @@ Boolean ONe_core(const struct stardata_t * const stardata,
Boolean star_is_LBV(const struct stardata_t * const stardata,
const struct star_t * const star);
Boolean star_is_WR(const struct star_t * const star);
+Boolean star_is_postAGB(struct stardata_t * const stardata,
+ struct star_t * const star);
double convective_turnover_time(const double radius,
struct star_t * star,
diff --git a/src/single_star_functions/star_is_postAGB.c b/src/single_star_functions/star_is_postAGB.c
new file mode 100644
index 000000000..b360c43df
--- /dev/null
+++ b/src/single_star_functions/star_is_postAGB.c
@@ -0,0 +1,30 @@
+#include "../binary_c.h"
+No_empty_translation_unit_warning;
+
+/*
+ * Identify a star as a post-(A)GB star
+ *
+ * We use Teff_postAGB_min to detect them,
+ * but note we ignore Teff_postAGB_max because
+ * often the Hurley prescription seems to make
+ * very hot stars that would not be classed
+ * as "post-AGB". These should be!
+ */
+
+Boolean star_is_postAGB(struct stardata_t * const stardata,
+ struct star_t * const star)
+{
+ const double teff =
+ Teff_from_star_struct(star);
+ const Boolean post_AGB =
+ (Boolean)
+ (
+ (star->stellar_type == GIANT_BRANCH ||
+ star->stellar_type == EAGB ||
+ star->stellar_type == TPAGB ||
+ WHITE_DWARF(star->stellar_type))
+ &&
+ teff > stardata->preferences->Teff_postAGB_min
+ );
+ return post_AGB;
+}
diff --git a/src/stellar_structure/stellar_structure_BSE.c b/src/stellar_structure/stellar_structure_BSE.c
index c629d0afa..8a8bb9caa 100644
--- a/src/stellar_structure/stellar_structure_BSE.c
+++ b/src/stellar_structure/stellar_structure_BSE.c
@@ -210,6 +210,7 @@ int stellar_structure_BSE_with_newstar(struct stardata_t * const stardata,
oldstar->radius,
stardata->model.intermediate_step
);
+
return 0;
}
#endif//BSE
diff --git a/src/wind/wind_AGB.c b/src/wind/wind_AGB.c
index 2d95c3d13..acc155823 100644
--- a/src/wind/wind_AGB.c
+++ b/src/wind/wind_AGB.c
@@ -13,8 +13,15 @@ static double wind_TPAGB(WIND_PROTOTYPE_ARGS);
double wind_AGB(WIND_PROTOTYPE_ARGS)
{
double mdot_AGB;
+ const Boolean post_AGB =
+ star_is_postAGB(stardata,star);
- if(stellar_type==EAGB)
+ if(post_AGB == TRUE &&
+ stardata->preferences->postagbwind != POSTAGB_WIND_GIANT)
+ {
+ mdot_AGB = wind_postAGB(WIND_CALL_ARGS);
+ }
+ else if(stellar_type==EAGB)
{
mdot_AGB = wind_EAGB(WIND_CALL_ARGS);
}
diff --git a/src/wind/wind_GB.c b/src/wind/wind_GB.c
index 1e09568cc..aa7b13df3 100644
--- a/src/wind/wind_GB.c
+++ b/src/wind/wind_GB.c
@@ -11,11 +11,21 @@ double Pure_function wind_GB(WIND_PROTOTYPE_ARGS)
* "first giant branch" winds: AGB winds are dealt with in the AGB
* function: all other potentially giant stars have this wind applied
*/
+ const Boolean post_AGB =
+ star_is_postAGB(stardata,star);
- if((stellar_type==HERTZSPRUNG_GAP)||
- (stellar_type==GIANT_BRANCH)||
- (stellar_type==CHeB)||
- (stellar_type>=HeMS))
+ if(post_AGB == TRUE &&
+ stardata->preferences->postagbwind != POSTAGB_WIND_GIANT)
+ {
+ /*
+ * Use post-AGB mass loss rate
+ */
+ mdot_GB = wind_postAGB(WIND_CALL_ARGS);
+ }
+ else if((stellar_type==HERTZSPRUNG_GAP)||
+ (stellar_type==GIANT_BRANCH)||
+ (stellar_type==CHeB)||
+ (stellar_type>=HeMS))
{
if(stardata->preferences->PN_fast_wind &&
diff --git a/src/wind/wind_Hurley2002.c b/src/wind/wind_Hurley2002.c
index b65387f07..b895ff89b 100644
--- a/src/wind/wind_Hurley2002.c
+++ b/src/wind/wind_Hurley2002.c
@@ -70,9 +70,29 @@ double wind_Hurley2002(WIND_PROTOTYPE_ARGS)
else
{
/*
- * no mass loss for remnants
+ * no mass loss for remnants, except post-AGB stars
*/
- mdot_wind = 0.0;
+ if(WHITE_DWARF(stellar_type))
+ {
+ mdot_wind = wind_WD(WIND_CALL_ARGS);
+ }
+ else
+ {
+ mdot_wind = 0.0;
+ }
+ }
+
+ if(star->starnum==0 && stardata->common.ndiscs > 0)
+ {
+ printf("DISC %30.20e %g %g %g %g %g %g\n",
+ stardata->model.time,
+ mdot_wind,
+ mdot[MDOT_GB],
+ mdot[MDOT_AGB],
+ mdot[MDOT_WR],
+ mdot[MDOT_OTHER],
+ mdot[MDOT_LBV]
+ );
}
return mdot_wind;
}
diff --git a/src/wind/wind_Schneider.c b/src/wind/wind_Schneider.c
index ab18cef81..6f41d831c 100644
--- a/src/wind/wind_Schneider.c
+++ b/src/wind/wind_Schneider.c
@@ -92,6 +92,20 @@ double wind_Schneider(WIND_PROTOTYPE_ARGS)
mdot_dyn,
mdot_tenpc);
}
+ else
+ {
+ /*
+ * no mass loss for remnants, except post-AGB stars
+ */
+ if(WHITE_DWARF(stellar_type))
+ {
+ mdot_Schneider = wind_WD(WIND_CALL_ARGS);
+ }
+ else
+ {
+ mdot_Schneider = 0.0;
+ }
+ }
return mdot_Schneider;
}
diff --git a/src/wind/wind_WD.c b/src/wind/wind_WD.c
new file mode 100644
index 000000000..fcbcb8149
--- /dev/null
+++ b/src/wind/wind_WD.c
@@ -0,0 +1,15 @@
+#include "../binary_c.h"
+No_empty_translation_unit_warning;
+#include "wind.h"
+
+/*
+ * Mass-loss rate for white dwarfs
+ */
+
+double wind_WD(WIND_PROTOTYPE_ARGS)
+{
+ return Max(
+ stardata->preferences->postagbwind == POSTAGB_WIND_GIANT ? 0.0 : wind_postAGB(WIND_CALL_ARGS),
+ 0.0
+ );
+}
diff --git a/src/wind/wind_binary_c_2020.c b/src/wind/wind_binary_c_2020.c
index 54b06cbdd..8eeef6b80 100644
--- a/src/wind/wind_binary_c_2020.c
+++ b/src/wind/wind_binary_c_2020.c
@@ -11,9 +11,11 @@ No_empty_translation_unit_warning;
double wind_binary_c_2020(WIND_PROTOTYPE_ARGS)
{
double mdot;
-
-
- if(stardata->preferences->wind_mass_loss == WIND_ALGORITHM_BINARY_C_2020 &&
+ if(WHITE_DWARF(stellar_type))
+ {
+ 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)
{
diff --git a/src/wind/wind_mass_loss_rate.c b/src/wind/wind_mass_loss_rate.c
index 70fe8910d..4c8f3b2a4 100644
--- a/src/wind/wind_mass_loss_rate.c
+++ b/src/wind/wind_mass_loss_rate.c
@@ -62,17 +62,17 @@ double wind_mass_loss_rate(WIND_PROTOTYPE_ARGS)
* 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
- );
+ 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];
+ mdot_wind *= stardata->preferences->wind_multiplier[star->stellar_type];
- return mdot_wind;
+ return mdot_wind;
}
diff --git a/src/wind/wind_postAGB.c b/src/wind/wind_postAGB.c
new file mode 100644
index 000000000..57b223e41
--- /dev/null
+++ b/src/wind/wind_postAGB.c
@@ -0,0 +1,93 @@
+#include "../binary_c.h"
+No_empty_translation_unit_warning;
+#include "wind.h"
+
+
+double wind_postAGB(WIND_PROTOTYPE_ARGS)
+{
+ /*
+ * Post-(A)GB mass loss rate
+ */
+ double mdot=0.0;
+ const double teff = Teff_from_star_struct(star);
+
+ if(star_is_postAGB(stardata, star) == TRUE)
+ {
+ if(stardata->preferences->postagbwind == POSTAGB_WIND_GIANT)
+ {
+ Exit_binary_c(BINARY_C_ALGORITHM_BRANCH_FAILURE,
+ "Despite wanting to use the giant-branch wind, for some reason wind_postAGB has been called. This is a bug.");
+ }
+ else if(stardata->preferences->postagbwind == POSTAGB_WIND_NONE)
+ {
+ /*
+ * No wind
+ */
+ mdot = 0.0;
+ }
+ else if(stardata->preferences->postagbwind == POSTAGB_WIND_KRTICKA2020)
+ {
+ /*
+ * Krticka, Kubát and Krticková (2020, A&A 635, A173)
+ *
+ * Eq. 1. Note that I have introduced factors of 1e-9
+ * and converted temperatures to K, and their index
+ * runs from 1 to 4, while ours runs from 0 to 3,
+ * both inclusive.
+ *
+ * As they say, this is valid for 10-117kK only, so I
+ * limit application to Teff>10kK, but at the high
+ * temperature end I limit the Teff used in the calculation
+ * to 117kK.
+ *
+ * Also, white dwarfs require T>105kK to launch a wind
+ * (see their section 7.3) so I also put in that limit.
+ */
+ if(
+ (
+ stellar_type < HeWD &&
+ teff > 10e3)
+ ||
+ (
+ WHITE_DWARF(stellar_type) &&
+ teff > 105e3
+ )
+ )
+ {
+ const double teff_use = Min(teff,117e3);
+ const double alpha = 1.63;
+ const double zeta = 0.53;
+ const double Zsun = 0.0134;
+ static const double mi[] = {2.51e-9, 0.72e-9, 0.48e-9, 0.33e-9};
+ static const double Ti[] = {18.3e3, 46.4e3, 84.6e3, 116.7e3};
+ static const double deltaTi[] = {3.63e3, 17.7e3, 20.4e3, 1.7e3};
+ mdot = 0.0;
+ for(int i=0; i<4; i++)
+ {
+ mdot += mi[i] * exp(-Pow2((teff_use - Ti[i])/deltaTi[i]));
+ }
+ mdot *=
+ pow(star->luminosity*1e-3,alpha)
+ *
+ pow(stardata->common.effective_metallicity/Zsun, zeta);
+ }
+ else
+ {
+ mdot = 0.0;
+ }
+ }
+ else
+ {
+ Exit_binary_c(BINARY_C_ALGORITHM_BRANCH_FAILURE,
+ "Unknown post-AGB wind %d attempted.",
+ stardata->preferences->postagbwind);
+ }
+ }
+
+ printf("KRITCKA %30.20e %g %g %g\n",
+ stardata->model.time,
+ teff,
+ star->luminosity,
+ mdot);
+ return mdot;
+}
diff --git a/src/wind/wind_prescriptions.def b/src/wind/wind_prescriptions.def
index 41d029b11..5c8dbc5f7 100644
--- a/src/wind/wind_prescriptions.def
+++ b/src/wind/wind_prescriptions.def
@@ -56,6 +56,11 @@
X(GOLDMAN_ETAL_2017 ) \
X(BEASOR_ETAL_2020 )
+#define POSTAGB_WIND_PRESCRIPTIONS_LIST \
+ X(GIANT) \
+ X(NONE) \
+ X(KRTICKA2020)
+
#define ROTATIONALLY_ENHANCED_MASSLOSS_LIST \
X(NONE) \
X(LANGER_FORMULA) \
diff --git a/src/wind/wind_prescriptions.h b/src/wind/wind_prescriptions.h
index 32dcc3632..5405dfe7e 100644
--- a/src/wind/wind_prescriptions.h
+++ b/src/wind/wind_prescriptions.h
@@ -67,6 +67,12 @@ enum { EAGB_WIND_PRESCRIPTIONS_LIST };
#define X(CODE) TPAGB_WIND_##CODE ,
enum { TPAGB_WIND_PRESCRIPTIONS_LIST };
+/* different post-(A)GB winds */
+#undef X
+#define X(CODE) POSTAGB_WIND_##CODE ,
+enum { POSTAGB_WIND_PRESCRIPTIONS_LIST };
+
+
/* rotationally enhanced mass-loss prescriptions */
#undef X
diff --git a/src/wind/wind_prototypes.h b/src/wind/wind_prototypes.h
index 9907dece7..101ac977b 100644
--- a/src/wind/wind_prototypes.h
+++ b/src/wind/wind_prototypes.h
@@ -7,6 +7,9 @@ double wind_mass_loss_rate(WIND_PROTOTYPE_ARGS);
double wind_other(WIND_PROTOTYPE_ARGS) Constant_function;
double wind_GB(WIND_PROTOTYPE_ARGS) Pure_function;
double wind_AGB(WIND_PROTOTYPE_ARGS);
+double wind_WD(WIND_PROTOTYPE_ARGS);
+double wind_postAGB(WIND_PROTOTYPE_ARGS);
+
double Constant_function_if_no_debug wind_LBV(WIND_PROTOTYPE_ARGS);
double Pure_function wind_WR(WIND_PROTOTYPE_ARGS);
void wind_multiplicative_factors(double * const mdot,
diff --git a/tbse b/tbse
index ddadd207b..7472c373f 100755
--- a/tbse
+++ b/tbse
@@ -1048,6 +1048,15 @@ EAGBWIND=EAGB_WIND_BSE
# TPAGB_WIND_BEASOR_ETAL_2020 = 10 = Beasor et al. (2020)
TPAGBWIND=TPAGB_WIND_VW93_KARAKAS
+# post-(A)GB wind:
+# POSTAGB_WIND_GIANT = 0 : use the prescription used on the giant branch
+# POSTAGB_WIND_NONE = 1 : no wind on the post-AGB
+# POSTAGB_WIND_KRTICKA2020 = 2 : Krticka, Kubát and Krticková (2020, A&A 635, A173)
+#
+# The post-AGB is identified as having Teff > Teff_postAGB_min
+POSTAGBWIND=POSTAGB_WIND_GIANT
+TEFF_POSTAGB_MIN=6000.0
+
# implement a fixed, fast wind for the PN phase
#
# ...DM... are the envelope masses below which these are applied
@@ -2224,6 +2233,8 @@ DEFAULT_ARGS=" \
--gbwind $GBWIND \
--eagbwind $EAGBWIND \
--tpagbwind $TPAGBWIND \
+--postagbwind $POSTAGBWIND \
+--Teff_postAGB_min $TEFF_POSTAGB_MIN \
--wind_gas_to_dust_ratio $WIND_GAS_TO_DUST_RATIO \
--gb_reimers_eta $GB_REIMERS_ETA \
--probability $PROBABILITY \
--
GitLab