diff --git a/src/binary_macros.h b/src/binary_macros.h
index 7e41017cd46839d4bbb2e357b7b11af71fbdbc18..8b6d6edfb816108eee6bacf370e30ed5b1f93689 100644
--- a/src/binary_macros.h
+++ b/src/binary_macros.h
@@ -19,7 +19,7 @@
#include "wind/wind_macros.h"
#include "derivatives.h"
#include "stellar_colours/stellar_colour_macros.h"
-
+#include "sn.h"
/* buffering states */
#define INTERNAL_BUFFERING_OFF 0
@@ -70,9 +70,9 @@
#define MEMMAP_STRING_LENGTH 82
#endif /* MEMMAP */
-#ifdef LOG_SUPERNOVAE
-#include "sn.h"
-#endif //LOG_SUPERNOVAE
+
+
+
#define HRDIAG_FILENAME "/tmp/c_binary_hrdiag.dat"
#include "stellar_structure/stellar_structure_macros.h"
diff --git a/src/binary_maths.h b/src/binary_maths.h
index 66041e2b813ddcb7ce7ab984bebdbe93d11a5d61..485820ad05beffccf697ad1053b2ee8297b8ee0f 100644
--- a/src/binary_maths.h
+++ b/src/binary_maths.h
@@ -95,6 +95,7 @@ static inline double func_min(const double a,const double b){return(MIN_MACRO(a,
/* max, min of several floating point numbers */
#define MAX3(A,B,C) MAX((A),MAX((B),(C)))
#define MAX4(A,B,C,D) MAX((A),MAX3((A),(B),(C)))
+#define MAX5(A,B,C,D,E) MAX((A),MAX4((A),(B),(C),(D)))
#define MIN3(A,B,C) MIN((A),MIN((B),(C)))
#define MIN4(A,B,C,D) MIN((D),MIN((A),MIN((B),(C))))
#define MIN5(A,B,C,D,E) MIN((E),MIN((D),MIN((A),MIN((B),(C)))))
diff --git a/src/binary_parameters.h b/src/binary_parameters.h
index 5b8f38fe627cd8c029958c3414197342cd92508a..1393dae3b09b9ed0ea381796557251cada8687a5 100644
--- a/src/binary_parameters.h
+++ b/src/binary_parameters.h
@@ -53,9 +53,6 @@
*/
//#define LIKE_AMANDA
-// Enable Axel's changes : WILL FAIL
-//#define AXEL
-
/*
* Enable SMOOTH_AGB to smooth the HR-diagram through the transition
* from overshooting EAGB to non-overshooting TPAGB model fits.
@@ -639,10 +636,6 @@
/* selma's stuff for MS binary RLOF*/
//#define SELMA
-/* If AXEL is defined, AGB smoothing should not be defined */
-#ifdef AXEL
-#undef SMOOTH_AGB
-#endif
/*
* On the cluster you do not want a FILE_LOG, it just slows
@@ -664,7 +657,6 @@
/*
* Enable this to use Axel's updated (and much better!) RGB radius function.
- * (requries rebuild if you change it)
*/
#define AXEL_RGBF_FIX
diff --git a/src/binary_structures.h b/src/binary_structures.h
index 03b83a8383317bc8694cdc781041aa068160ed57..2f86386e89ea67af0912d3cf4b59711a999f5c89 100644
--- a/src/binary_structures.h
+++ b/src/binary_structures.h
@@ -523,10 +523,6 @@ struct preferences_t {
int batch_submode;
Boolean no_arg_checking;
#endif//BATCHMODE
-#ifdef AXEL
- // overshooting TPAGB
- Boolean overshooting_tpagb;
-#endif//AXEL
Boolean rlperi;
#ifdef ROTATIONALLY_ENHANCED_MASS_LOSS
int rotationally_enhanced_mass_loss;
diff --git a/src/common_envelope/common_envelope_evolution.c b/src/common_envelope/common_envelope_evolution.c
index 00a8a28ed541f29b579de3950e98cbb5922eb1d4..01d9fab14efdcadb2d0732f11d3c85855941fae4 100644
--- a/src/common_envelope/common_envelope_evolution.c
+++ b/src/common_envelope/common_envelope_evolution.c
@@ -238,7 +238,8 @@ int common_envelope_evolution (double *m01,
stardata->star[nstar1].Xenv,
#endif//NUCSYN
stardata);
- //_exit(0);
+ fflush(NULL);
+ _exit(0);
#endif
Append_logstring(LOG_COMENV_PRE,
diff --git a/src/common_envelope/common_envelope_polytrope.c b/src/common_envelope/common_envelope_polytrope.c
index 0d131e691435d2721f71bc2861f905e46975b256..29940255a3306ffe1b36dd16783b707a10bce561 100644
--- a/src/common_envelope/common_envelope_polytrope.c
+++ b/src/common_envelope/common_envelope_polytrope.c
@@ -14,7 +14,7 @@
#ifdef COMENV_POLYTROPES
-//#define COMENV_POLYTROPES_DEBUG
+#define COMENV_POLYTROPES_DEBUG
#include "common_envelope_polytrope.h"
@@ -50,7 +50,7 @@ double common_envelope_polytrope(
#ifdef NUCSYN
nucsyn_recombination_energy_per_gram(stardata,X);
#else
- RYDBERG_CONSTANT / PROTON_MASS;
+ 1.3 * RYDBERG_CONSTANT / PROTON_MASS;
#endif
double GMMR = - G * m * me / r;
@@ -79,7 +79,7 @@ double common_envelope_polytrope(
#endif
double MAYBE_UNUSED err = x[3];
- double q = mc / (4.0 * PI * beta * POW3(alpha) );
+ double q = mc / (4.0 * PI * beta * POW3(alpha) ); // core density
double K = 4.0 * PI * G * POW2(alpha) / ((1.0 + n) * pow(q, 1.0/n - 1.0));
diff --git a/src/logging/log_every_timestep.c b/src/logging/log_every_timestep.c
index 04f696897f051dd3f5b500e5e114a02c1d211fa4..d3b15c88e9067da4ec25886ad7452b45426b6606 100644
--- a/src/logging/log_every_timestep.c
+++ b/src/logging/log_every_timestep.c
@@ -2205,13 +2205,14 @@ void log_every_timestep(struct stardata_t * RESTRICT stardata)
#endif
stardata);
- PRINTF("COMENVLOG %g %d %g %g %g %g\n",
+ PRINTF("COMENVLOG %g %d %g %g %g %g %g\n",
stardata->model.time,
star->stellar_type,
star->mass,
star->radius,
lambda,
- star->k2
+ star->k2,
+ star->core_mass
);
}
}
diff --git a/src/nucsyn/nucsyn_init_first_pulse.c b/src/nucsyn/nucsyn_init_first_pulse.c
index a84b88e87a1ecabf7cd9a0dc7717bb1c345d21ae..08fea18af3a535b4397e02d4d535bfc3a2eed805 100644
--- a/src/nucsyn/nucsyn_init_first_pulse.c
+++ b/src/nucsyn/nucsyn_init_first_pulse.c
@@ -76,11 +76,7 @@ void nucsyn_init_first_pulse(double * RESTRICT tagb,
/* save initial envelope mass and initial mass */
star->menv0 = (*mt)-(*mc);
-#ifdef AXEL
- *m0 = *mt;
-#else
*m0 = *mass;
-#endif
/*
* initially set interpulse period (assuming lambda=0 - which is true
diff --git a/src/nucsyn/nucsyn_mcmin_for_dredgeup.c b/src/nucsyn/nucsyn_mcmin_for_dredgeup.c
index e7734ec5664f8e22a4a00400949f1bd1952ac605..4eb82b3fa47d34936fcf572f2cfdb018c0618b7c 100644
--- a/src/nucsyn/nucsyn_mcmin_for_dredgeup.c
+++ b/src/nucsyn/nucsyn_mcmin_for_dredgeup.c
@@ -99,11 +99,8 @@ double nucsyn_mcmin_for_dredgeup(const double m,
/* Now to fix for M>Msdu-0.5 */
-#ifdef AXEL
- y=MAX(star->mc0,MIN(0.7,y));
-#else
y=MAX(Karakas2002_mc1tp(m,z,star,stardata),MIN(0.7,y));
-#endif
+
/* DUP calibration */
y += stardata->preferences->delta_mcmin;
diff --git a/src/nucsyn/nucsyn_recombination_energy_per_gram.c b/src/nucsyn/nucsyn_recombination_energy_per_gram.c
index 58bb242bb43d1b2124a594b8a11e37c1f4d43974..85744770a62155a8899b761daa5a9c0c5b6fa525 100644
--- a/src/nucsyn/nucsyn_recombination_energy_per_gram.c
+++ b/src/nucsyn/nucsyn_recombination_energy_per_gram.c
@@ -1,4 +1,7 @@
#include "../binary_c.h"
+
+#ifdef NUCSYN
+
#include "nucsyn_elements.h"
/*
* Calculate the recombination energy per gram
@@ -195,3 +198,4 @@ year = {2015},
return (E/M);
}
+#endif // NUCSYN
diff --git a/src/perl/modules_targz/Binning-0.03.tar.gz b/src/perl/modules_targz/Binning-0.03.tar.gz
index d54f9dacb75f1c93c9e4db1a4f3a00ce39f0fde2..426a21f0cf827a5f65c06cc87a10e244ef12ec19 100644
Binary files a/src/perl/modules_targz/Binning-0.03.tar.gz and b/src/perl/modules_targz/Binning-0.03.tar.gz differ
diff --git a/src/perl/modules_targz/Data-Serializer-RobJSON-0.05.tar.gz b/src/perl/modules_targz/Data-Serializer-RobJSON-0.05.tar.gz
index e5f56f54cbe51630d41df514f5c0c44626c3d337..e58751bdd59cc6d50f7e28acc8be216f7c03f377 100644
Binary files a/src/perl/modules_targz/Data-Serializer-RobJSON-0.05.tar.gz and b/src/perl/modules_targz/Data-Serializer-RobJSON-0.05.tar.gz differ
diff --git a/src/perl/modules_targz/Hash-RobMerge-0.14.tar.gz b/src/perl/modules_targz/Hash-RobMerge-0.14.tar.gz
index 9641f42045a46a45bb6488f164264a4d75541219..6e9b42e220b4c0e93f5d985603ee87dda2eea259 100644
Binary files a/src/perl/modules_targz/Hash-RobMerge-0.14.tar.gz and b/src/perl/modules_targz/Hash-RobMerge-0.14.tar.gz differ
diff --git a/src/perl/modules_targz/Histogram-0.01.tar.gz b/src/perl/modules_targz/Histogram-0.01.tar.gz
index f19e0ab811fc32cd1fa1e1401ee59c3f59136138..88e4947c8508e1395a0c84ce7aee4141fe2d5d68 100644
Binary files a/src/perl/modules_targz/Histogram-0.01.tar.gz and b/src/perl/modules_targz/Histogram-0.01.tar.gz differ
diff --git a/src/perl/modules_targz/IMF-0.05.tar.gz b/src/perl/modules_targz/IMF-0.05.tar.gz
index 6aeaa2a1508879e49c7fe93d0b365bface520789..d0b4af4d76b05a7d8b526fbae08102bfa25c416f 100644
Binary files a/src/perl/modules_targz/IMF-0.05.tar.gz and b/src/perl/modules_targz/IMF-0.05.tar.gz differ
diff --git a/src/perl/modules_targz/Maths_Double-0.01.tar.gz b/src/perl/modules_targz/Maths_Double-0.01.tar.gz
index 3f2d8a98f14a7d159427f265f5fb1a3696d44ab7..dca9333eb2122fea948dcd464f81e809dddadf08 100644
Binary files a/src/perl/modules_targz/Maths_Double-0.01.tar.gz and b/src/perl/modules_targz/Maths_Double-0.01.tar.gz differ
diff --git a/src/perl/modules_targz/RobInterpolation-0.04.tar.gz b/src/perl/modules_targz/RobInterpolation-0.04.tar.gz
index b30ac2576c9761249b2f9eb8a7b89960c35e4ce0..eb2658bb6cb45f9be9c5a149d63901b45530a524 100644
Binary files a/src/perl/modules_targz/RobInterpolation-0.04.tar.gz and b/src/perl/modules_targz/RobInterpolation-0.04.tar.gz differ
diff --git a/src/perl/modules_targz/binary_grid-2.00.tar.gz b/src/perl/modules_targz/binary_grid-2.00.tar.gz
index 38597d94bd52e3292ad37b76ce85fa345ede2a53..19739762d74f3849d1bdad7f5052264047c52aff 100644
Binary files a/src/perl/modules_targz/binary_grid-2.00.tar.gz and b/src/perl/modules_targz/binary_grid-2.00.tar.gz differ
diff --git a/src/perl/modules_targz/binary_stars-0.04.tar.gz b/src/perl/modules_targz/binary_stars-0.04.tar.gz
index 7e42141ed96299ef8c2561d18f7895edbf404f6f..c5ca6704fded0b58dcb3d111234b982bde1e467b 100644
Binary files a/src/perl/modules_targz/binary_stars-0.04.tar.gz and b/src/perl/modules_targz/binary_stars-0.04.tar.gz differ
diff --git a/src/perl/modules_targz/cosmology-0.01.tar.gz b/src/perl/modules_targz/cosmology-0.01.tar.gz
index 08aa5136b936d5cfd1059318d09bdbd5b5a2c2fd..b073399c0f8113f702823c4db8f533dcc58686ee 100644
Binary files a/src/perl/modules_targz/cosmology-0.01.tar.gz and b/src/perl/modules_targz/cosmology-0.01.tar.gz differ
diff --git a/src/perl/modules_targz/distribution_functions-0.05.tar.gz b/src/perl/modules_targz/distribution_functions-0.05.tar.gz
index 397187ec9dd100915283ef73a1d064956418a587..589d5c5ae2d11ec7f880fc9256a85b03efea545a 100644
Binary files a/src/perl/modules_targz/distribution_functions-0.05.tar.gz and b/src/perl/modules_targz/distribution_functions-0.05.tar.gz differ
diff --git a/src/perl/modules_targz/rob_misc-0.13.tar.gz b/src/perl/modules_targz/rob_misc-0.13.tar.gz
index 63b25e047424929ff1ea2cc5a3ffc0a5dd8acd79..59772d5a9801dc54025e83ceb9043b0b62d289f4 100644
Binary files a/src/perl/modules_targz/rob_misc-0.13.tar.gz and b/src/perl/modules_targz/rob_misc-0.13.tar.gz differ
diff --git a/src/perl/modules_targz/robqueue-0.05.tar.gz b/src/perl/modules_targz/robqueue-0.05.tar.gz
index 98ba955c6d1717226e8e1cb83cd2d4d9a79a04ee..5828be7869c9a40647c2d32de4f2d2c679b24d3f 100644
Binary files a/src/perl/modules_targz/robqueue-0.05.tar.gz and b/src/perl/modules_targz/robqueue-0.05.tar.gz differ
diff --git a/src/perl/modules_targz/spacing_functions-0.02.tar.gz b/src/perl/modules_targz/spacing_functions-0.02.tar.gz
index 348e0db2350150aeb390238ceb48d0d9731752b2..056015541e42f87b9fd5faba03d18457b89dd264 100644
Binary files a/src/perl/modules_targz/spacing_functions-0.02.tar.gz and b/src/perl/modules_targz/spacing_functions-0.02.tar.gz differ
diff --git a/src/setup/reset_preferences.c b/src/setup/reset_preferences.c
index 5c060bc0e17c967d2b4d355359287454e65f6e24..fd40f876aca1661677367ebd1293aec8cb06044f 100644
--- a/src/setup/reset_preferences.c
+++ b/src/setup/reset_preferences.c
@@ -185,9 +185,6 @@ void reset_preferences(struct preferences_t * RESTRICT preferences)
#endif
-#ifdef AXEL
- preferences->overshooting_tpagb=TRUE;
-#endif
#ifdef VW93_MIRA_SHIFT
preferences->vw93_mira_shift=0.0;
#endif
diff --git a/src/setup/version.c b/src/setup/version.c
index 82f8dd589601b446ef816afd61411871cd688b4a..d252ee6cb500e8773b9ecaa22c589745fd475b87 100644
--- a/src/setup/version.c
+++ b/src/setup/version.c
@@ -507,11 +507,8 @@ void version(struct stardata_t * RESTRICT stardata)
PRINTF("ACQUIRE_TEST_FOR_NEGATIVE_AGES is not defined\n");
#endif
-#ifdef AXEL
- PRINTF("AXEL is defined\n");
-#else
- PRINTF("AXEL is not defined\n");
-#endif
+ PRINTF("AXEL is not defined (it is depreacted)\n");
+
#ifdef ROBSMOOTH
PRINTF("ROBSMOOTH is defined\n");
#ifdef ROBSMOOTH_STOP_AT_TPAGB
diff --git a/src/stellar_structure/stellar_structure_AGB.c b/src/stellar_structure/stellar_structure_AGB.c
index 9bc7f8d06f251970170efd3787f58e85986747da..36f0157995953bc727fa88a5e0270ad7ea01ef6f 100644
--- a/src/stellar_structure/stellar_structure_AGB.c
+++ b/src/stellar_structure/stellar_structure_AGB.c
@@ -1,5 +1,7 @@
#include "../binary_c.h"
#include "stellar_structure_debug.h"
+//#undef DEBUG
+//#define DEBUG 0
Stellar_type stellar_structure_AGB(struct new_stellar_structure_t * s,
double *mcx,
@@ -36,10 +38,7 @@ Stellar_type stellar_structure_AGB(struct new_stellar_structure_t * s,
#ifdef SMOOTH_AGB
*mcx = mass_CO_core_at_BAGB(s->phase_start_mass,stardata->common.metallicity);
#else //SMOOTH_AGB
-#ifdef AXEL
- *mcx = mass_CO_core_at_BAGB(s->phase_start_mass,stardata->common.metallicity);
- Dprint("Axel's CO core mass at BAGB = %12.12e\n",*mcx);
-#else //AXEL
+
*mcx = mcgbtf(tbagb,
s->GB[8],
s->GB,
@@ -48,8 +47,6 @@ Stellar_type stellar_structure_AGB(struct new_stellar_structure_t * s,
s->timescales[T_EAGB_T],
s->phase_start_mass
);
-
-#endif //AXEL
#endif //SMOOTH_AGB
@@ -59,12 +56,15 @@ Stellar_type stellar_structure_AGB(struct new_stellar_structure_t * s,
* of stardata->preferences->chandrasekhar_mass. Of course these might be exactly the
* same thing...
*/
- mcmax = MAX(MAX(stardata->preferences->max_tpagb_core_mass,(0.7730*mcbagb-0.350)),1.050*(*mcx));
+ mcmax = MAX4(
+ stardata->preferences->max_tpagb_core_mass,
+ 0.7730*mcbagb-0.350,
+ 1.050* *mcx,
+ (s->stellar_type == TPAGB ? stardata->preferences->max_tpagb_core_mass : 0.0)
+ );
Dprint("Compare mcx=%g to mcxmax=%g\n",*mcx,mcmax);
- if(s->stellar_type==TPAGB) mcmax=stardata->preferences->max_tpagb_core_mass;
-
/*
* NB real TPAGB evolution is now done separately by the nucsyn
* routine nucsyn_tpagb (if NUCSYN is defined),
@@ -75,15 +75,13 @@ Stellar_type stellar_structure_AGB(struct new_stellar_structure_t * s,
/*
* Set core mass, luminosity and age
*/
-#if (defined(NUCSYN) && defined(ROBSMOOTHDEBUG))
+#if defined NUCSYN && defined ROBSMOOTHDEBUG
printf("ROBSMOOTH Check TPAGB age=%g tscl=%g mcx=%g core_mass=%g\n",
s->age,s->timescales[T_TPAGB_FIRST_PULSE],*mcx,s->core_mass);
#endif
- if(((s->age<s->timescales[T_TPAGB_FIRST_PULSE])
-// ||(s->stellar_type==CHeB)
- )&&
- (s->stellar_type!=TPAGB))
+ if(s->age < s->timescales[T_TPAGB_FIRST_PULSE] &&
+ s->stellar_type!=TPAGB)
{
/* EAGB star */
convert_to_helium_star=
@@ -126,50 +124,30 @@ Stellar_type stellar_structure_AGB(struct new_stellar_structure_t * s,
{
#ifndef NUCSYN
/* NUCSYN is not defined */
-
s->radius = ragbf(s->mass,
- s->luminosity,
- metallicity_parameters[ZPAR_MASS_HE_FLASH],
- giant_branch_parameters
+ s->luminosity,
+ metallicity_parameters[ZPAR_MASS_HE_FLASH],
+ giant_branch_parameters
#ifdef AXEL_RGBF_FIX
- ,stardata->common.metallicity
+ ,stardata->common.metallicity
#endif // AXEL_RGBF_FIX
);
+
#else // NUCSYN
/* NUCSYN is defined */
-#ifdef AXEL
- s->radius = ragbf(s->phase_start_mass,
- s->luminosity,
- metallicity_parameters[ZPAR_MASS_HE_FLASH],
- giant_branch_parameters
-#ifdef AXEL_RGBF_FIX
- ,stardata->common.metallicity
-#endif // AXEL_RGBF_FIX
- );
-
- s->radius *= pow(s->mass/ s->phase_start_mass,-0.3333333333);
-
- Dprint("Axel radius %12.12e from mass=%12.12e lum=%12.12e zpar[2]=%12.12e\n",s->radius,
- s->phase_start_mass,s->luminosity,metallicity_parameters[ZPAR_MASS_HE_FLASH]);
-
-#else // AXEL
-
- /* AXEL not defined */
-
-
if(s->stellar_type==EAGB)
{
/* EAGB */
s->radius = ragbf(s->mass,
- s->luminosity,
- metallicity_parameters[ZPAR_MASS_HE_FLASH],
- giant_branch_parameters
+ s->luminosity,
+ metallicity_parameters[ZPAR_MASS_HE_FLASH],
+ giant_branch_parameters
#ifdef AXEL_RGBF_FIX
- ,stardata->common.metallicity
+ ,stardata->common.metallicity
#endif // AXEL_RGBF_FIX
);
@@ -216,14 +194,14 @@ Stellar_type stellar_structure_AGB(struct new_stellar_structure_t * s,
/* estimate the radius at the first thermal pulse (from Amanda) */
double r1tp_Karakas=Karakas2002_ragb_teff(s->mass,
- stardata->common.metallicity,
- star->l1tp_guess, // set in stellar_structure_EAGB
- s->phase_start_mass,
- star->mc1tp_guess,
- s->mass-star->mc1tp_guess,
- star->mc1tp_guess,
- star,
- stardata);
+ stardata->common.metallicity,
+ star->l1tp_guess, // set in stellar_structure_EAGB
+ s->phase_start_mass,
+ star->mc1tp_guess,
+ s->mass-star->mc1tp_guess,
+ star->mc1tp_guess,
+ star,
+ stardata);
{
/* scale the radius to match at the first thermal pulse */
@@ -266,10 +244,10 @@ Stellar_type stellar_structure_AGB(struct new_stellar_structure_t * s,
s->radius *=
pow((1-((s->age-tbagb)/(s->timescales[T_TPAGB_FIRST_PULSE]-tbagb))),1.0)+
Karakas2002_ragb_teff(s->mass,
- stardata->common.metallicity,
- s->luminosity,s->phase_start_mass,s->core_mass,
- s->mass-s->core_mass,
- s->core_mass,star,stardata) *
+ stardata->common.metallicity,
+ s->luminosity,s->phase_start_mass,s->core_mass,
+ s->mass-s->core_mass,
+ s->core_mass,star,stardata) *
pow(((s->age-tbagb)/(s->timescales[T_TPAGB_FIRST_PULSE]-tbagb)),1.0);
#endif //SMOOTH_AGB
@@ -279,13 +257,12 @@ Stellar_type stellar_structure_AGB(struct new_stellar_structure_t * s,
{
/* TPAGB : always use fit to Amanda's models */
s->radius = Karakas2002_ragb_teff(s->mass,
- stardata->common.metallicity,
- s->luminosity,s->phase_start_mass,star->mc0,
- s->mass-s->core_mass,
- star->mc_nodup,star,stardata);
-
- }
-#endif//AXEL
+ stardata->common.metallicity,
+ s->luminosity,s->phase_start_mass,star->mc0,
+ s->mass-s->core_mass,
+ star->mc_nodup,star,stardata);
+
+ }
Dprint("stellar_structure_AGB: stellar type = %d, r=%12.12e (L=%g)\n",s->stellar_type,s->radius,s->luminosity);
diff --git a/src/stellar_structure/stellar_structure_EAGB.c b/src/stellar_structure/stellar_structure_EAGB.c
index 0b07845866c3e0ee606c451ba72a64e6e1699136..3b0a6c69a9943159b9c666346039247b92c4e5d0 100644
--- a/src/stellar_structure/stellar_structure_EAGB.c
+++ b/src/stellar_structure/stellar_structure_EAGB.c
@@ -33,19 +33,6 @@ Boolean stellar_structure_EAGB(struct new_stellar_structure_t * s,
s->timescales[8],
s->timescales[9]);
#else //SMOOTH_AGB
-#ifdef AXEL
- star->mcx_EAGB = mcx_eagbf(s->age,
- s->GB[8],
- GB,
- s->timescales[7],
- s->timescales[8],
- s->timescales[9]);
- *mcx = star->mcx_EAGB;
-
- s->core_mass = MAX(mcbagb,star->mcx_EAGB);
- Dprint("Axel's CO core mass as f(t) = %12.12e so core mass = %12.12e\n",star->mcx_EAGB,s->core_mass);
-
-#else//AXEL
*mcx = mcgbtf(s->age,
s->GB[8],
s->GB,
@@ -97,7 +84,6 @@ Boolean stellar_structure_EAGB(struct new_stellar_structure_t * s,
s->core_mass=mcbagb;
#endif // NUCSYN
-#endif // AXEL
#endif //SMOOTH_AGB
Dprint("STELLAR_STRUCTURE type=%d mc=mcbagb=%12.12e mcx=%12.12e\n",
@@ -106,11 +92,11 @@ Boolean stellar_structure_EAGB(struct new_stellar_structure_t * s,
*mcx);
-#if (defined(SMOOTH_AGB)||defined(AXEL))
+#if defined SMOOTH_AGB
s->luminosity = lum_mc_EAGB(star->mcx_EAGB,GB);
#else
s->luminosity = lmcgbf(*mcx,s->GB);
-#endif // SMOOTH_AGB or AXEL
+#endif // SMOOTH_AGB
//Dprint("STAR Set luminosity=%12.12e (stellar type %d mcx=%g)\n",s->luminosity,s->stellar_type,*mcx);
star->mcx_EAGB = *mcx;
@@ -199,7 +185,7 @@ Boolean stellar_structure_EAGB(struct new_stellar_structure_t * s,
Dprint("mt=%12.12e core=%12.12e (if mt<core > TPAGB)\n",s->mass,s->core_mass);
- if(LESS_OR_EQUAL(s->mass,s->core_mass))
+ if(LESS_OR_EQUAL(s->mass, s->core_mass))
{
/*
* Evolved naked helium star as the envelope is lost but the
@@ -213,36 +199,18 @@ Boolean stellar_structure_EAGB(struct new_stellar_structure_t * s,
#ifndef NO_UNEXPECTED_MASS_CHANGES
s->mass = s->core_mass;
#endif
-
s->phase_start_mass = s->mass;
#ifdef SMOOTH_AGB
s->core_mass = star->mcx_EAGB;
#else //SMOOTH_AGB
-#ifdef AXEL
- s->core_mass = star->mcx_EAGB;
-
- Dprint("Axel's CO core mass on the EAGB %12.12e\n",
- s->core_mass);
-#else//AXEL
-
- // NB use star->stellar NOT s->stellar_type, as that is
- // definitely post-AGB!
- if(star->stellar_type>TPAGB)
- {
- /*
- * The core cannot shrink if we're post-AGB as
- * the CO core has already been built up!
- */
- s->core_mass = MAX(s->core_mass,*mcx);
- }
- else
- {
- s->core_mass = *mcx;
- }
-
-#endif//AXEL
+ /*
+ * The core cannot shrink if we're post-AGB as
+ * the CO core has already been built up!
+ */
+ s->core_mass = star->stellar_type>TPAGB ? MAX(s->core_mass,*mcx) : *mcx;
+
#endif //SMOOTH_AGB
call_calc_lt2(s->stellar_type,s->phase_start_mass,s->mass);
diff --git a/src/stellar_structure/stellar_structure_HeStar.c b/src/stellar_structure/stellar_structure_HeStar.c
index b97fe69d93d3e46f5b6f36a08f4ab0e7b2a67e17..48800b3993ee8c8f23070de6b6275424b24a37c8 100644
--- a/src/stellar_structure/stellar_structure_HeStar.c
+++ b/src/stellar_structure/stellar_structure_HeStar.c
@@ -351,7 +351,9 @@ void stellar_structure_HeStar(struct new_stellar_structure_t * s,
*/
if(j>jlim)
{
+#ifdef NUCSYN
star->sn_last_time=SN_GRB_COLLAPSAR;
+#endif
star->SN_type = SN_GRB_COLLAPSAR;
}
}
diff --git a/src/stellar_structure/stellar_structure_TPAGB.c b/src/stellar_structure/stellar_structure_TPAGB.c
index 4dbd9ffced4fd6df0dc5fc980a85d79bd1797292..3c6ac493086c9ecda13f813d810dd5e38afadcef 100644
--- a/src/stellar_structure/stellar_structure_TPAGB.c
+++ b/src/stellar_structure/stellar_structure_TPAGB.c
@@ -22,14 +22,6 @@ void stellar_structure_TPAGB(struct new_stellar_structure_t * s,
double mc1tp = Karakas2002_mc1tp(s->phase_start_mass,stardata->common.metallicity,star,stardata);
#ifdef NUCSYN_SECOND_DREDGE_UP
-#ifdef AXEL
- // check this!
- if(!RLOF_CALL &&
- star->second_dredge_up == FALSE &&
- mc1tp < mcbagb &&
- mc1tp < stardata->preferences->chandrasekhar_mass &&
- s->stellar_type == EAGB)
-#else
Dprint("2DUP check : if (RLOF_CALL=%d==0) && (star->second_dredge_up=%d=0) && %12.12e < mcbagb=%12.12e < %12.12e (st=%d)\n",
RLOF_CALL,
star->second_dredge_up,
@@ -39,7 +31,6 @@ void stellar_structure_TPAGB(struct new_stellar_structure_t * s,
mcbagb > mc1tp &&
mcbagb < stardata->common.max_mass_for_second_dredgeup &&
(s->stellar_type == EAGB || s->stellar_type == TPAGB))
-#endif
{
if(star->second_dredge_up==FALSE)
@@ -59,9 +50,7 @@ void stellar_structure_TPAGB(struct new_stellar_structure_t * s,
}
#endif /* NUCSYN_SECOND_DREDGE_UP */
-#endif /* NUCSYN */
-
-#ifdef NUCSYN
+
/*
* Call the tpagb function to deal with the evolution, since
* it is closely linked to the nucleosynthesis! *
@@ -81,40 +70,11 @@ void stellar_structure_TPAGB(struct new_stellar_structure_t * s,
s->timescales[7],
s->timescales[7],
s->timescales[9]);
-#else //SMOOTH_AGB
-#ifdef AXEL
- /*
- * core mass evolution is dealt with on the TPAGB by the
- * nucleosynthesis TPAGB code, however it still has to be
- * set at the beginning of the TPAGB. For continuity we use
- * the same function as for the overshooting model.
- */
- s->core_mass = mcx_eagbf(s->age,
- s->GB[8],
- s->GB,
- s->timescales[7],
- s->timescales[7],
- s->timescales[9]);
-
-
- Dprint("Axel's core mass on EAGB (post2dup!) = %12.12e from age=%12.12e s->GB[8]=%12.12e s->timescales7,8,9=%12.12e %12.12e %12.12e\n",
- s->core_mass,
- s->age,s->GB[8],
- s->timescales[7],
- s->timescales[7],
- s->timescales[9]
- );
-
-#endif // AXEL
#endif // SMOOTH_AGB
Dprint("COREDIFF mc1tp=%12.12e core_mass=%12.12e\n",
mc1tp,s->core_mass);
-
-
-#ifndef AXEL
-
if(star->num_thermal_pulses<1)
{
/*
@@ -127,11 +87,13 @@ void stellar_structure_TPAGB(struct new_stellar_structure_t * s,
/* Use Hurley core mass at the 1st TP */
#else
/* Enforce Amanda's core masses at the 1st TP */
- s->core_mass = MIN(mc1tp,s->core_mass);
+ s->core_mass = MIN(mc1tp,
+ s->core_mass);
#endif
/* Second dredge up may reduce the core mass */
- s->core_mass = MIN(s->core_mass, guess_mc1tp(s->phase_start_mass,star,stardata));
+ s->core_mass = MIN(s->core_mass,
+ guess_mc1tp(s->phase_start_mass,star,stardata));
/*
printf("COMPACT%d %g %g %d compact=%g mc=%g\n",
@@ -146,11 +108,11 @@ void stellar_structure_TPAGB(struct new_stellar_structure_t * s,
/* prevent 2nd DUP into the CO core */
if(stardata->model.intpol==0) s->core_mass = MAX(star->compact_core_mass,s->core_mass);
}
-#endif //AXEL
+
Dprint("TPAGB BEGIN mc=%12.12e\n",s->core_mass);
- /* now set the stellar type and evolve on the TPAGB */
+ /* set the stellar type to TPAGB */
s->stellar_type = TPAGB;
if(!RLOF_CALL)
@@ -192,8 +154,6 @@ void stellar_structure_TPAGB(struct new_stellar_structure_t * s,
#else
/* Old way of doing things - somewhat dubious */
s->stellar_type = TPAGB;
-
- double mcwas=s->core_mass;
s->core_mass = mcgbtf(s->age,
s->GB[2],
s->GB,
@@ -202,14 +162,12 @@ void stellar_structure_TPAGB(struct new_stellar_structure_t * s,
s->timescales[T_TPAGB_T],
s->phase_start_mass
);
-
s->luminosity = lmcgbf(s->core_mass,s->GB);
/*
* Approximate 3rd Dredge-up on AGB by limiting Core_Mass.
*/
double lambda = MIN(0.90,(0.30+0.0010*POW5(s->phase_start_mass)));
-
*mcx = mcgbtf(s->timescales[T_TPAGB_FIRST_PULSE],
s->GB[2],
s->GB,
@@ -225,5 +183,4 @@ void stellar_structure_TPAGB(struct new_stellar_structure_t * s,
#endif /* NUCSYN */
*mcx = s->core_mass;
*mcmax = MIN(s->mass,*mcmax);
-
}
diff --git a/src/stellar_structure/stellar_structure_remnant_and_perturbations.c b/src/stellar_structure/stellar_structure_remnant_and_perturbations.c
index 9a90d319923eebf90d9ddc3185dbac779de94de8..5b1ba334bd9cd3c1843300d0db522bd8bcca6b8c 100644
--- a/src/stellar_structure/stellar_structure_remnant_and_perturbations.c
+++ b/src/stellar_structure/stellar_structure_remnant_and_perturbations.c
@@ -92,13 +92,7 @@ void stellar_structure_remnant_and_perturbations(struct new_stellar_structure_t
#ifdef SMOOTH_AGB
remnant_luminosity = lum_mc_EAGB(star->mcx_EAGB,GB);
#else //SMOOTH_AGB
-#ifdef AXEL
- remnant_luminosity = lum_mc_EAGB(star->mcx_EAGB,GB);
- Dprint("Axel's EAGB luminosity %12.12e\n",remnant_luminosity);
-
-#else //AXEL
remnant_luminosity = lmcgbf(mcx,s->GB);
-#endif//AXEL
#endif//SMOOTH_AGB
if(tau<1.0)
diff --git a/src/stellar_timescales/stellar_timescales_massive_AGB.c b/src/stellar_timescales/stellar_timescales_massive_AGB.c
index e2d37ea2f782c720d680ae50315832d209450e0a..83a4a20f644f2b44834098136eece3e2fb2c74a2 100644
--- a/src/stellar_timescales/stellar_timescales_massive_AGB.c
+++ b/src/stellar_timescales/stellar_timescales_massive_AGB.c
@@ -43,9 +43,6 @@ void stellar_timescales_massive_AGB_tn(Stellar_type stellar_type,
*tn = timescales[T_TPAGB_TINF_1] - (1.0/((GB[GB_SMOOTH_2]-1.0)*GB[GB_A_H_HE]*GB[GB_SMOOTH_1]))*pow(MAX(mc1,mc1_jarr),(1.0-GB[GB_SMOOTH_2]));
#else //SMOOTH_AGB
-#ifdef AXEL
- *tn = timescales[T_TPAGB_TINF_1] - (1.0/((GB[GB_AXEL_2]-1.0)*GB[GB_A_H_HE]*GB[GB_AXEL_1]))*pow(mc1,(1.0-GB[GB_AXEL_2]));
-#else
if(LESS_OR_EQUAL(mc1,GB[GB_Mx]))
{
*tn = timescales[T_TPAGB_TINF_1] - (ip/(GB[GB_A_H_HE]*GB[GB_D]))*pow(mc1,mp);
@@ -54,7 +51,5 @@ void stellar_timescales_massive_AGB_tn(Stellar_type stellar_type,
{
*tn = timescales[T_TPAGB_TINF_2] - (iq/(GB[GB_A_H_HE]*GB[GB_B]))*pow(mc1,mq);
}
-
-#endif // AXEL
#endif // SMOOTH_AGB
}
diff --git a/src/stellar_timescales/stellar_timescales_post_HG.c b/src/stellar_timescales/stellar_timescales_post_HG.c
index 636f670ea5ff6b3c09920af4d1af7b0a1b3e3d64..ee9f10b6254edf0f7e171a1e395f7b142c772f9c 100644
--- a/src/stellar_timescales/stellar_timescales_post_HG.c
+++ b/src/stellar_timescales/stellar_timescales_post_HG.c
@@ -199,54 +199,6 @@ int stellar_timescales_post_HG(double tbagb,
#else //SMOOTH_AGB
-#ifdef AXEL
- /*
- * Use phase_start_mass NOT mass because otherwise there is
- * a problem with the start of the TPAGB. This
- * is a reasonable assumption: there is not usually
- * much mass loss on the EAGB due to winds and if
- * there is RLOF then we'll probably not get to the
- * TPAGB anyway...
- */
- Dprint("L 1TP (AXEL) from phase_start_mass=%g mc1=%g menv=%g\n",phase_start_mass,mc1,phase_start_mass-mc1);
- luminosities[L_TPAGB_FIRST_PULSE] =
- Karakas2002_lumfunc(
- phase_start_mass,
- mc1,
- phase_start_mass - mc1,
- mc1,
- 0.0,
- metallicity,
- stardata,
- star
- );
-
- /*
- * GB numbers for EAGB : used to follow L as f(t) on the EAGB (assuming little
- * mass loss during that stage)
- */
- GB[GB_AXEL_2]=log(luminosities[L_TPAGB_FIRST_PULSE]/luminosities[L_BAGB])/
- log(mc1/mass_CO_core_at_BAGB(phase_start_mass,metallicity));
- GB[GB_AXEL_1]=luminosities[L_TPAGB_FIRST_PULSE]*pow(mc1,-GB[GB_AXEL_2]);
-
- /*
- * EAGB time parameters
- */
- tbagb = timescales[T_HE_IGNITION] + timescales[T_HE_BURNING];
-
- y2=GB[GB_AXEL_2]-1.0;
- y=y2/GB[GB_AXEL_2];
-
-#ifndef USB_BSE_TIMESCALES_H
- // NB these are overridden from the values given above
-
- timescales[T_EAGB_TINF_2]=timescales[T_EAGB_TINF_1] = tbagb + (1.0/(y2*GB[GB_A_HE]*GB[GB_AXEL_1]))*
- pow(GB[GB_AXEL_1]/luminosities[L_BAGB],y);
-
- timescales[T_EAGB_T] = timescales[T_EAGB_TINF_1] - (timescales[T_EAGB_TINF_1] - tbagb)*
- pow(luminosities[L_BAGB]/luminosities[L_LMX],y);
-#endif
-#else//AXEL
/* AXEL not defined */
Dprint("No AXEL defined\n");
@@ -285,7 +237,7 @@ int stellar_timescales_post_HG(double tbagb,
Dprint("L1tp = %g from mass=%g mc1=%g c.f. jarrod=lmcgbf(%g)=%g\n",
luminosities[L_TPAGB_FIRST_PULSE],
mass,mc1,mc1,lmcgbf(mc1,GB));
-#endif //AXEL
+
#endif //SMOOTH_AGB
#endif // NUCSYN
@@ -296,12 +248,7 @@ int stellar_timescales_post_HG(double tbagb,
#else //SMOOTH_AGB
-#ifdef AXEL
- timescales[T_TPAGB_FIRST_PULSE] = timescales[T_EAGB_TINF_1] -
- (1.0/(y2*GB[GB_A_HE]*GB[GB_AXEL_1]))*pow(mc1, 1.0-GB[GB_AXEL_2]);
-#else//AXEL
- /* jarrod method */
if(LESS_OR_EQUAL(mc1,GB[GB_Mx]))
{
timescales[T_TPAGB_FIRST_PULSE] = timescales[T_EAGB_TINF_1] - (ip/(GB[GB_A_HE]*GB[GB_D]))*pow(mc1,mp);
@@ -310,7 +257,7 @@ int stellar_timescales_post_HG(double tbagb,
{
timescales[T_TPAGB_FIRST_PULSE] = timescales[T_EAGB_TINF_2] - (iq/(GB[GB_A_HE]*GB[GB_B]))*pow(mc1,mq);
}
-#endif//AXEL
+
#endif //SMOOTH_AGB
Dprint("First pulse at %g\n",timescales[T_TPAGB_FIRST_PULSE]);
@@ -359,15 +306,6 @@ int stellar_timescales_post_HG(double tbagb,
timescales[T_EAGB_T]);
#else //SMOOTH_AGB
-#ifdef AXEL
- mc2 = mcx_eagbf(timescales[T_HE_IGNITION] + timescales[T_HE_BURNING],
- GB[GB_A_HE],
- GB,
- timescales[T_EAGB_TINF_1],
- timescales[T_EAGB_TINF_1],
- timescales[T_EAGB_T]);
-
-#else
mc2 = mcgbtf(timescales[T_HE_IGNITION] + timescales[T_HE_BURNING],
GB[GB_A_HE],
@@ -377,7 +315,7 @@ int stellar_timescales_post_HG(double tbagb,
timescales[T_EAGB_T],
phase_start_mass
);
-#endif//AXEL
+
#endif//SMOOTH_AGB
@@ -392,9 +330,6 @@ int stellar_timescales_post_HG(double tbagb,
#else //SMOOTH_AGB
-#ifdef AXEL
- timescales[T_TMCMAX] = timescales[T_EAGB_TINF_1] - (1.0/((GB[GB_AXEL_2]-1.0)*GB[GB_A_HE]*GB[GB_AXEL_1]))*pow(mcmax,(1.0-GB[GB_AXEL_2]));
-#else //AXEL
if(LESS_OR_EQUAL(mcmax,GB[GB_Mx]))
{
timescales[T_TMCMAX] = timescales[T_EAGB_TINF_1] - (ip/(GB[GB_A_HE]*GB[GB_D]))*pow(mcmax,mp);
@@ -403,7 +338,6 @@ int stellar_timescales_post_HG(double tbagb,
{
timescales[T_TMCMAX] = timescales[T_EAGB_TINF_2] - (iq/(GB[GB_A_HE]*GB[GB_B]))*pow(mcmax,mq);
}
-#endif //AXEL
#endif //SMOOTH_AGB
}
else
@@ -420,7 +354,6 @@ int stellar_timescales_post_HG(double tbagb,
mcmax = (mcmax - lambda*mc1)/(1.0 - lambda);
#ifndef SMOOTH_AGB
-#ifndef AXEL
if(LESS_OR_EQUAL(mcmax,GB[GB_Mx]))
{
@@ -430,7 +363,7 @@ int stellar_timescales_post_HG(double tbagb,
{
timescales[T_TMCMAX] = timescales[T_TPAGB_TINF_2] - (iq/(GB[GB_A_H_HE]*GB[GB_B]))*pow(mcmax,mq);
}
-#endif // AXEL
+
#endif //SMOOTH_AGB
}
timescales[T_TMCMAX] = MAX(tbagb,timescales[T_TMCMAX]);
diff --git a/src/timestep/timestep_EAGB.c b/src/timestep/timestep_EAGB.c
index 58cd294e38c1a8f92478331019a6ac6b940dde60..a890ecc842924680334366d949727ed16d9f56bf 100644
--- a/src/timestep/timestep_EAGB.c
+++ b/src/timestep/timestep_EAGB.c
@@ -31,17 +31,7 @@ void timestep_EAGB(Timestep_prototype_args)
#else // SMOOTH_AGB
- /* probably using ROBSMOOTH or AXEL */
-
-#ifdef AXEL
- *time_remaining = tscls[13] - age;
- if(*time_remaining<0.2)
- {
- Set_timestep(*dt,*time_remaining,star,DT_LIMIT_EAGB_AXEL);
- }
-#else // AXEL
*time_remaining = MIN(tn,tscls[13]) - age;
-#endif //AXEL
#endif //SMOOTH_AGB
diff --git a/src/zfuncs/Karakas2002_mc1tp.c b/src/zfuncs/Karakas2002_mc1tp.c
index 08f6c2712d1d419bc607f4f76ac0b68aad100f41..a549a9f692c65701321afa3de3b9107098336855 100644
--- a/src/zfuncs/Karakas2002_mc1tp.c
+++ b/src/zfuncs/Karakas2002_mc1tp.c
@@ -13,9 +13,6 @@ double Karakas2002_mc1tp(double m,
double m_in=m,metallicity_in=metallicity; // mass,metallicity passed in
double y;
-#ifdef AXEL
- double mcbagb,mc1tp_ov,mc1tp_nov;
-#endif
Dprint("find mc1tp from m=%g z=%g star->phase_start_mass=%g\n",m,metallicity,star->phase_start_mass);
@@ -128,43 +125,6 @@ double Karakas2002_mc1tp(double m,
Dprint("Amanda's interpolate gave y=%g (m=%g z=%g)\n",y,m,metallicity);
#endif
-#ifdef AXEL
- if(stardata->preferences->overshooting_tpagb==TRUE)
- {
- // Implement Axel's changes to include overshooting
- // cores on the TPAGB
-
- // helium core mass at the start of the EAGB
- // = helium core mass at the end of the EAGB
- mcbagb=mcagbf(star->phase_start_mass,star,stardata);
- Dprint("Axel mcbagb=%g from phase_start_mass=%g\n",mcbagb,star->phase_start_mass);
-
- // Amanda's non-overshooting fit to mc1tp
- mc1tp_nov=y;
-
- // Overshooting fit to mc1tp
- mc1tp_ov=0.44*mcbagb + 0.448;
-
- // now what?
- // The following expression works thus:
-
- // Choose the smaller of mcbagb and mc1tp_ov since
- // the TPAGB must be when the CO core reaches mcbagb
- // or perhaps smaller if there is 2nd dredge-up but it
- // cannot be larger.
-
- // If mc1tp_nov is larger than this value then believe it
- // because we trust Amanda's models more than we trust
- // the overshooting models at low mass (when this is
- // likely to occur).
- y=MAX(MIN(mcbagb,mc1tp_ov),
- mc1tp_nov);
-
- Dprint("Axel core mass %g from MAX(MIN(mcbagb=%g,mc1tp_ov=%g),mc1tp_nov=%g)\n",
- y,mcbagb,mc1tp_ov,mc1tp_nov);
- }
-#endif//AXEL
-
/* we have no fits above M~6.5, use Jarrod's */
double fl=1.0/(1.0+pow(0.0001,(star->phase_start_mass-7.0))); // smooth 'step' function
@@ -177,7 +137,6 @@ double Karakas2002_mc1tp(double m,
Dprint("MC1TP interpolate A=%g, J=%g -> ",y,mc1tpj);
y=(1.0-fl)*y + fl * mc1tpj;
Dprint("%g\n",y);
-
#ifdef NANCHECKS
if(isnan(y))
@@ -185,7 +144,6 @@ double Karakas2002_mc1tp(double m,
Exit_binary_c(EXIT_NAN,"Mc1tp NAN\n");
}
#endif
-
return(y);
}
diff --git a/src/zfuncs/ragbf.c b/src/zfuncs/ragbf.c
index 298794dc093df17e704b4e23a91f296d55e8df19..2741cde660949c6d346a666e75f35ee177e08c45 100644
--- a/src/zfuncs/ragbf.c
+++ b/src/zfuncs/ragbf.c
@@ -2,6 +2,11 @@
#ifdef AXEL_RGBF_FIX
+/*
+ * Axel's modified ragbf function, based on BSE
+ * but updated.
+ */
+
double ragbf(const double m, /* mass */
const double lum,/* luminosity */
const double mhelflsh, /* helium flash mass */
@@ -205,6 +210,11 @@ double ragbf_1(const double m, /* mass */
#else // AXEL_RGBF_FIX
+
+/*
+ * The original RGBF function from BSE
+ */
+
double ragbf(const double m, /* mass */
const double lum,/* luminosity */
const double mhelf, /* helium flash mass */
diff --git a/tbse b/tbse
index 66cd069ca994450756af3f837f8cb1eeb1aa9110..5db5a595d7397473d16bea293f9fb90291d13e0d 100755
--- a/tbse
+++ b/tbse
@@ -83,7 +83,7 @@ ECC=0.0
M1=3; M2=1; PER=30
# AGB test
-M1=3; M2=1; PER=1000000
+#M1=4; M2=1; PER=1000000
# max evolution time in Myr (WMAP suggests 13.7Gyr)
MAX_EVOLUTION_TIME=10000
@@ -1297,6 +1297,7 @@ $HE_IG_EXTRA \
--sn_kick_companion_AIC_BH $SN_KICK_COMPANION_AIC_BH \
--sn_kick_companion_BH_BH $SN_KICK_COMPANION_BH_BH \
--sn_kick_companion_BH_NS $SN_KICK_COMPANION_BH_NS \
+--AGB_method $AGB_METHOD \
$@"
fi
diff --git a/test_dewi_tauris.pl b/test_dewi_tauris.pl
index 8efd714832e64a67d68e2033561e01f637fd9175..24ed9c61d7839d150770f36feaeab2d3ca78dc8a 100755
--- a/test_dewi_tauris.pl
+++ b/test_dewi_tauris.pl
@@ -32,39 +32,67 @@ sub parse_bse
last if($header eq 'fin');
if($header eq 'COMENVLOG')
{
- print {$datafp} "@$la\n";
+ #if($la->[1] == 5 || $la->[1] == 6)
+ #if($la->[1] != 4)
+ {
+ print {$datafp} "@$la\n";
+ }
}
}
}
# reconstruct their Fig.1 which is lamdba as a
# function of R/Rsun for various masses
-my $gpfile = "$outdir/Fig1.plt";
-open(my $gp,'>',$gpfile)||die;
-print $gp "set terminal pdfcairo
-set output \"$outdir/Fig1.pdf\"
-set xlabel \"Radius / Rsun\"
+my $gpfile = {
+ dewi=>"$outdir/Dewi_Tauris_Fig1.plt",
+ vds=>"$outdir/VDSluys_Fig4.plt"
+};
+my $gp = {};
+
+foreach my $k (keys %$gpfile)
+{
+ open($gp->{$k},'>',$gpfile->{$k})||die;
+
+ my $pdf = $gpfile->{$k};
+ $pdf =~ s/plt$/pdf/;
+
+ print {$gp->{$k}} "set terminal pdfcairo
+set output \"$pdf\"
set ylabel \"Lambda\"
-set y2label \"k^{2}\"
-set xrange[0:500]
-set yrange[0.1:100]
-set logscale y
set xtics
set ytics nomirror
+set pointsize 0.2
+";
+}
+
+print {$gp->{dewi}} "
+set xrange[0:500]
+set yrange[0.1:100]
+set xlabel \"Radius / Rsun\"
+set y2label \"k^{2}\"
set y2tics nomirror
+set logscale y
";
-#foreach my $m (3..10)
-foreach my $m (3)
+print {$gp->{vds}} "
+set xrange[0.15:1.25]
+set yrange[0:5]
+set xlabel \"Core mass / Msun\"
+";
+
+# Dewi plots
+foreach my $m (3..10)
+#foreach my $m (4)
{
- print $gp "
+ print {$gp->{dewi}} "
set title \"M=$m Msun\"
";
my $plotcomma='plot ';
- foreach my $lambda_both (0,1)
+ foreach my $lambda_both (0.0, 1.0)
{
my $lambda_ionisation = $lambda_both;
my $lambda_enthalpy = $lambda_both;
- my $f = "$outdir/$m.$lambda_both.dat";
+ my $f = "$outdir/dewi-$m.$lambda_both.dat";
+
open($datafp, '>', $f)||die;
print "LAM $lambda_both\n";
@@ -85,14 +113,61 @@ set title \"M=$m Msun\"
max_evolution_time => 15000,
})
);
+
close $datafp;
- print $gp "$plotcomma \"$f\" u 4:5 axes x1y1 w l notitle , \"$f\" u 4:6 axes x1y2 w l dashtype \".-\" notitle ";
+ print {$gp->{dewi}} "$plotcomma \"$f\" u 4:5 axes x1y1 w lp notitle , \"$f\" u 4:6 axes x1y2 w lp dashtype \".-\" notitle ";
$plotcomma=',';
}
+}
+# VDS plots : masses are from Fig.3
+foreach my $m (0.91, 1.01, 1.14, 1.30, 1.48, 1.63, 1.81, 2.00, 2.46, 2.79, 3.17, 3.70, 4.09, 4.65, 5.28, 6.00, 6.82)
+{
+ print {$gp->{vds}} "
+set title \"M=$m Msun\"
+";
+ my $plotcomma='plot ';
+ foreach my $lambda_both (0.0, 1.0)
+ {
+ my $lambda_ionisation = $lambda_both;
+ my $lambda_enthalpy = $lambda_both;
+ my $f = "$outdir/vds-$m.$lambda_both.dat";
+
+ open($datafp, '>', $f)||die;
+ print "LAM $lambda_both\n";
+
+ $population->set(
+ lambda_ionisation=>$lambda_ionisation,
+ lambda_enthalpy=>$lambda_enthalpy,
+ );
+
+ $population->tbse(
+ $population->make_evcode_arghash
+ ({
+ probability=>1,
+ M_1 => $m,
+ M_2 => 0.1,
+ orbital_period => 1e9,
+ metallicity => 0.02,
+ eccentricity => 0.0,
+ max_evolution_time => 15000,
+ })
+ );
+
+ close $datafp;
+ print {$gp->{vds}} "$plotcomma \"$f\" u 7:5 axes x1y1 w lp notitle ";
+ $plotcomma=',';
+ }
}
-close $gp;
-`gnuplot $gpfile`;
+
+
+foreach my $k (keys %$gpfile)
+{
+ close $gp->{$k};
+ `gnuplot $gpfile->{$k}`;
+}
+
+
exit;
@@ -133,22 +208,11 @@ sub defaults
# The name of the binary_c executable
'prog' => 'binary_c-CND7',
- # VARIABLES TO BE CHANGED : maybe set on the command line!
'binary' => 1, # single stars -> 0, binary stars -> 1
- 'CRAP_parameter' => 1e3, # 1e3 = Tout
- 'WRLOF_method' => 1, # 1 = q-dep, 2=quadratic
+ 'CRAP_parameter' => 0, # 1e3 = Tout
+ 'WRLOF_method' => 0, # 1 = q-dep, 2=quadratic
'RLOF_method' => 3, # 3 = Claeys et al 2014
- # age range for star formation in the thick disc
- # Snaith et al 2015 suggest start=13Gyr, end=8Gyr
- # and is ~ constant (and large!)
- 'thick_disc_start_age' => 10e3,
- 'thick_disc_end_age' => 4e3,
-
- # logg range : set to -+100.0 for all stars
- 'thick_disc_logg_min'=>-100.0,
- 'thick_disc_logg_max'=>+100.0,
-
# we want array refs, not lines, from tbse_line
return_array_refs => 1,
@@ -158,5 +222,7 @@ sub defaults
lambda_enthalpy => 1.0,
comenv_splitmass=>1.0,
+ AGB_method => 0
+
);
}