From e406e75b0745f43a7515f28ae4922aafc1a3b87c Mon Sep 17 00:00:00 2001
From: Robert Izzard <r.izzard@surrey.ac.uk>
Date: Mon, 7 Aug 2017 16:30:48 +0000
Subject: [PATCH] stillbroken
---
src/binary_maths.h | 9 ++
src/binary_parameters.h | 2 +-
...ar_mass_and_angular_momentum_derivatives.c | 41 ++++++--
src/disc/cbdisc_eccentricity_pumping_rate.c | 31 ++++++-
src/disc/disc.h | 2 +-
src/single_star_functions/giant_age.c | 2 +-
src/stellar_structure/stellar_structure_BSE.c | 83 ++++++++++++-----
.../stellar_structure_CHeB.c | 93 ++++++++++++++++---
.../stellar_structure_EAGB.c | 9 +-
src/stellar_structure/stellar_structure_HG.c | 3 -
src/stellar_structure/stellar_structure_RG.c | 21 ++---
src/timestep/timestep_CHeB.c | 29 +++++-
src/zfuncs/lgbtf.c | 22 -----
tbse | 11 +--
test_new.pl | 14 ++-
15 files changed, 274 insertions(+), 98 deletions(-)
diff --git a/src/binary_maths.h b/src/binary_maths.h
index 1bbdcd532..517c351a0 100644
--- a/src/binary_maths.h
+++ b/src/binary_maths.h
@@ -333,6 +333,8 @@ static inline int func_is_really_zero(const double a){return(IS_REALLY_ZERO_MACR
#define POW5d4_MACRO(X) ((X)*sqrt(sqrt(X)))
#define POW5p5_MACRO(X) (POW5_MACRO(X)*sqrt(X))
#define POW3d4_MACRO(X) ((X)/(POW1d4_MACRO(X)))
+#define POW2d3_MACRO(X) (POW2_MACRO(cbrt(X)))
+#define POW4d3_MACRO(X) (POW4_MACRO(cbrt(X)))
#ifdef POWER_OPERATIONS_ARE_FUNCTIONS
/* use functions for simple power operations */
@@ -351,6 +353,8 @@ static inline int func_is_really_zero(const double a){return(IS_REALLY_ZERO_MACR
#define POW1p5(A) func_pow1p5(A)
#define POW5p5(A) func_pow5p5(A)
#define POW3d4(A) func_pow3d4(A)
+#define POW2d3(A) func_pow2d3(A)
+#define POW4d3(A) func_pow4d3(A)
static inline double func_pow2(const double x) pure_function;
static inline double func_pow3(const double x) pure_function;
@@ -377,6 +381,8 @@ static inline double func_pow1d4(const double x) {return POW1d4_MACRO(x);}
static inline double func_pow5d4(const double x) {return POW5d4_MACRO(x);}
static inline double func_pow5p5(const double x) {return POW5p5_MACRO(x);}
static inline double func_pow3d4(const double x) {return POW3d4_MACRO(x);}
+static inline double func_pow2d3(const double x) {return POW2d3_MACRO(x);}
+static inline double func_pow4d3(const double x) {return POW4d3_MACRO(x);}
#else
/* use macros for simple power operations */
#define POW2(A) POW2_MACRO(A)
@@ -394,6 +400,9 @@ static inline double func_pow3d4(const double x) {return POW3d4_MACRO(x);}
#define POW1p5(A) POW1p5_MACRO(A)
#define POW5p5(A) POW5p5_MACRO(A)
#define POW3d4(A) POW3d4_MACRO(A)
+#define POW2d3(A) POW2d3_MACRO(A)
+#define POW4d3(A) POW4d3_MACRO(A)
+
#endif //POWER_OPERATIONS_ARE_FUNCTIONS
diff --git a/src/binary_parameters.h b/src/binary_parameters.h
index 572fecab6..9fddb4d87 100644
--- a/src/binary_parameters.h
+++ b/src/binary_parameters.h
@@ -1326,7 +1326,7 @@
/*
* Logging for individual systems with discs
*/
-//#define DISC_LOG
+#define DISC_LOG
//#define DISC_LOG_2D
diff --git a/src/binary_star_functions/apply_stellar_mass_and_angular_momentum_derivatives.c b/src/binary_star_functions/apply_stellar_mass_and_angular_momentum_derivatives.c
index 541fd3f12..b6b12480b 100644
--- a/src/binary_star_functions/apply_stellar_mass_and_angular_momentum_derivatives.c
+++ b/src/binary_star_functions/apply_stellar_mass_and_angular_momentum_derivatives.c
@@ -47,10 +47,17 @@ void apply_stellar_mass_and_angular_momentum_derivatives(struct stardata_t * sta
star->mass += star->dM_in_timestep;
- Dprint("NEWMASS %g, Mc = %g, stellar type %d\n",
+ if(star->starnum == 0)
+ printf("NEWMASS %g, Mc = %g, stellar type %d, Menv = %g, R = %g (mdot wind loss %g rlof gain %g) dt = %g Myr\n",
star->mass,
star->core_mass,
- star->stellar_type);
+ star->stellar_type,
+ star->mass - star->core_mass,
+ star->radius,
+ star->derivative[DERIVATIVE_STELLAR_MASS_WIND_LOSS],
+ star->derivative[DERIVATIVE_STELLAR_MASS_RLOF_GAIN],
+ stardata->model.dtm
+ );
#ifdef MINIMUM_STELLAR_MASS
star->mass = MAX(MINIMUM_STELLAR_MASS,
@@ -193,25 +200,32 @@ void apply_stellar_mass_and_angular_momentum_derivatives(struct stardata_t * sta
accretor->stellar_type == HeMS ? CHeB :
accretor->stellar_type == HeHG ? TPAGB :
TPAGB; //HeGB
- Dprint("Convert accreting star %d from stellar type %d to %d\n",
+ printf("Convert accreting star %d from stellar type %d to %d, accretion rate %g\n",
accretor->starnum,
was,
- accretor->stellar_type);
+ accretor->stellar_type,
+ accretion_rate);
double mcx;
+ printf("PRE AGE %g\n",accretor->age);
if(accretor->stellar_type==CHeB)
{
accretor->age /= accretor->tms;
}
/*
- * Core mass is the whole of the ex-helium star.
+ * Core mass is the whole of the helium star.
* The envelope is the accreted hydrogen.
*/
mcx = accretor->mass - star->dM_in_timestep;
- Dprint("New core mass is %g, M = %g\n",mcx,accretor->mass);
-
-
+ printf("New core mass is %g, M = %g, Menv = %g\n",
+ mcx,
+ accretor->mass,
+ accretor->mass - mcx
+ );
+
+ accretor->phase_start_mass = MAX(accretor->mass,
+ accretor->phase_start_mass);
giant_age(&mcx,
accretor->mass,
& accretor->stellar_type,
@@ -219,14 +233,21 @@ void apply_stellar_mass_and_angular_momentum_derivatives(struct stardata_t * sta
& accretor->age,
stardata,
accretor);
+ printf("post gnt age %d\n",accretor->stellar_type);
+ printf("POST AGE %g M = %g Mc = %g mcx = %g\n",
+ accretor->age,
+ accretor->mass,
+ accretor->core_mass,
+ mcx);
accretor->core_mass = mcx;
accretor->epoch = stardata->model.time +
stardata->model.dtm - accretor->age;
- Dprint("Set secondary epoch = %g + %g - %g = %g\n",
+ printf("Set secondary epoch = t=%g + dtm=%g - age=%g = epoch=%g\n",
stardata->model.time ,
- stardata->model.dtm , accretor->age,
+ stardata->model.dtm ,
+ accretor->age,
accretor->epoch);
}
#ifdef DISCS_ACCRETION_ONTO_RGB_FUDGE_BROKEN
diff --git a/src/disc/cbdisc_eccentricity_pumping_rate.c b/src/disc/cbdisc_eccentricity_pumping_rate.c
index 9fe8c5196..d48831a0d 100644
--- a/src/disc/cbdisc_eccentricity_pumping_rate.c
+++ b/src/disc/cbdisc_eccentricity_pumping_rate.c
@@ -4,7 +4,12 @@
#include "disc.h"
-
+static void resonance_radii(struct binary_system_t * binary,
+ double l,
+ double m,
+ double *r_corotation,
+ double *r_inner,
+ double *r_outer);
void cbdisc_eccentricity_pumping_rate(struct stardata_t * stardata,
struct disc_t * disc,
@@ -91,7 +96,13 @@ void cbdisc_eccentricity_pumping_rate(struct stardata_t * stardata,
stardata->preferences->cbdisc_resonance_multiplier;
+ if(DISC_DEBUG)
{
+
+ double rc,rin,rout;
+ resonance_radii(binary,l,m,&rc,&rin,&rout);
+ Discdebug(1,"ECC PUMP radii %g %g %g\n",
+ rc,rin,rout);
Discdebug(2,
"a=%g RJhalf=%g H/R(RJhalf)=%g\n",a,R,disc->HRJhalf);
Discdebug(2,
@@ -199,4 +210,22 @@ void cbdisc_eccentricity_pumping_rate(struct stardata_t * stardata,
#endif //NANCHECKS
}
+
+static void resonance_radii(struct binary_system_t * binary,
+ double l,
+ double m,
+ double *r_corotation,
+ double *r_inner,
+ double *r_outer)
+{
+ /*
+ * Calculate the radius at which a resonance occurs
+ */
+ *r_corotation = POW2d3( m / l ) * binary->separation;
+ *r_inner = POW2d3( (m - 1.0)/ l ) * binary->separation;
+ *r_outer = POW2d3( (m + 1.0)/ l ) * binary->separation;
+}
+
+
+
#endif // DISCS
diff --git a/src/disc/disc.h b/src/disc/disc.h
index 2935d35ff..997f9d2a1 100644
--- a/src/disc/disc.h
+++ b/src/disc/disc.h
@@ -5,7 +5,7 @@
* Header file for the binary_c disc library
*/
-#define DISC_DEBUG 0
+#define DISC_DEBUG 1
#include "disc_macros.h"
#include "disc_power_laws.h"
diff --git a/src/single_star_functions/giant_age.c b/src/single_star_functions/giant_age.c
index 24252716b..be8a7e4ee 100644
--- a/src/single_star_functions/giant_age.c
+++ b/src/single_star_functions/giant_age.c
@@ -176,7 +176,7 @@ double giant_age(double *mc, /** core mass in/out **/
{
/** we cannot have a fractional age that's not between 0 and 1 ! **/
CLAMP(*age,0.0,1.0);
-#ifdef GIANT_AGE_FATAL_ERRORS
+#ifdef GIANT_AGE_FATAL_ERRORS
Exit_binary_c(GIANT_AGE_FRACTIONAL_AGE_OUT_OF_BOUNDS,
" FATAL ERROR! giant_age: fractional age out of bounds : FAGE %f\n",
*age);
diff --git a/src/stellar_structure/stellar_structure_BSE.c b/src/stellar_structure/stellar_structure_BSE.c
index d8e405052..da6a55b1b 100644
--- a/src/stellar_structure/stellar_structure_BSE.c
+++ b/src/stellar_structure/stellar_structure_BSE.c
@@ -27,7 +27,10 @@
*/
#include "../binary_c.h"
#include "stellar_structure_debug.h"
-
+#ifdef NEW_CORE_GROWTH
+static Boolean helium_ignition_test(struct new_stellar_structure_t * s,
+ struct stardata_t *stardata);
+#endif
void stellar_structure_BSE(struct new_stellar_structure_t * s,
struct star_t *star,
struct stardata_t *stardata,
@@ -178,7 +181,13 @@ void stellar_structure_BSE(struct new_stellar_structure_t * s,
s->timescales[T_GIANT_TINF_1],
s->timescales[T_GIANT_TINF_2]
);
-
+
+
+#ifdef NEW_CORE_GROWTH
+ Boolean helium_ignition =
+ helium_ignition_test(s,stardata);
+#endif
+
/* do not allow CHeB stars to reverse their evolution */
if (s->age < s->timescales[T_BGB] &&
(s->stellar_type<CHeB ||
@@ -241,38 +250,25 @@ void stellar_structure_BSE(struct new_stellar_structure_t * s,
Dprint("Call stellar_structure_HG\n");
stellar_structure_HG(s,rg,star,stardata,caller_id);
}
- helium_burning = FALSE;
}
/*
* Now the GB, CHeB and AGB evolution.
*/
- else if(s->GB_core_mass
- <
- mcgbf(s->luminosities[L_HE_IGNITION],
- s->GB,
- s->luminosities[L_LMX]))
-
- //s->age < s->timescales[T_HE_IGNITION])
+
+
+#ifdef NEW_CORE_GROWTH
+ else if(s->stellar_type <= GIANT_BRANCH && IS_ZERO(s->CO_core_mass) && helium_ignition == FALSE &&
+ ! skip to CHeB
+ //else if(s->age < s->timescales[T_HE_IGNITION])
+#else
+ //else if(s->age < s->timescales[T_HE_IGNITION])
+#endif
{
/*
* Red Giant.
*/
Dprint("Call stellar_structure_RG\n");
-
- double Lig = s->luminosities[L_HE_IGNITION];
stellar_structure_RG(s,&rg,star,stardata,caller_id);
- double Mcig = mcgbf(s->luminosities[L_HE_IGNITION],
- s->GB,
- s->luminosities[L_LMX]);
- printf("HEIG %d L = %g cf %g Mc = %g CO %g GB %g MCig %g\n",
- s->stellar_type,
- s->luminosity,
- Lig,
- s->core_mass,
- s->CO_core_mass,
- s->GB_core_mass,
- Mcig
- );
}
else if(s->age < tbagb)
{
@@ -544,3 +540,42 @@ void stellar_structure_BSE(struct new_stellar_structure_t * s,
}
}
+#ifdef NEW_CORE_GROWTH
+static Boolean helium_ignition_test(struct new_stellar_structure_t * s,
+ struct stardata_t *stardata)
+{
+ /*
+ * Determine whether a HG/GB star has ignited helium
+ */
+ Boolean degenerate_core = LESS_OR_EQUAL(s->phase_start_mass,
+ stardata->common.metallicity_parameters[ZPAR_MASS_HE_FLASH]);
+
+ double Mc_helium_ignition =
+ mcgbf(s->luminosities[L_HE_IGNITION],
+ s->GB,
+ s->luminosities[L_LMX]);
+
+ double test_Mc =
+ degenerate_core ? s->core_mass :
+ (IS_ZERO(s->GB_core_mass) ? s->core_mass : s->GB_core_mass);
+
+ Boolean helium_ignition = Boolean_(test_Mc > Mc_helium_ignition);
+
+ if(0 && (s->stellar_type==HERTZSPRUNG_GAP ||
+ s->stellar_type==GIANT_BRANCH))
+ {
+ printf("MCHE %d Mc=%g McGB=%g : testMc=%g vs McHeIG=%g (Mx=%g) : AGE %g vs %g : ignite? %d\n",
+ s->stellar_type,
+ s->core_mass,
+ s->GB_core_mass,
+ test_Mc,
+ Mc_helium_ignition,
+ s->GB[GB_Mx],
+ s->age,
+ s->timescales[T_HE_IGNITION],
+ helium_ignition);
+ }
+
+ return helium_ignition;
+}
+#endif
diff --git a/src/stellar_structure/stellar_structure_CHeB.c b/src/stellar_structure/stellar_structure_CHeB.c
index b1a9bf97a..6ae2e712b 100644
--- a/src/stellar_structure/stellar_structure_CHeB.c
+++ b/src/stellar_structure/stellar_structure_CHeB.c
@@ -69,25 +69,95 @@ Stellar_type stellar_structure_CHeB(struct new_stellar_structure_t *s,
#endif
#endif /* NUCSYN_FIRST_DREDGE_UP and NUCSYN */
+ double mcbagb = mcagbf(s->phase_start_mass,giant_branch_parameters);
+
if(LESS_OR_EQUAL(s->phase_start_mass,
metallicity_parameters[ZPAR_MASS_HE_FLASH]))
{
- mcx=mcgbf(s->luminosities[L_HE_IGNITION],
- s->GB,
- s->luminosities[L_LMX]);
+ mcx = mcgbf(s->luminosities[L_HE_IGNITION],
+ s->GB,
+ s->luminosities[L_LMX]);
}
else
{
- mcx=mcheif(s->phase_start_mass,
- metallicity_parameters[ZPAR_MASS_HE_FLASH],
- metallicity_parameters[10],
- giant_branch_parameters);
+ mcx = mcheif(s->phase_start_mass,
+ metallicity_parameters[ZPAR_MASS_HE_FLASH],
+ metallicity_parameters[10],
+ giant_branch_parameters);
}
- double mcbagb = mcagbf(s->phase_start_mass,giant_branch_parameters);
-
- tau=(s->age - s->timescales[T_HE_IGNITION])/s->timescales[T_HE_BURNING];
+
+ tau = (s->age - s->timescales[T_HE_IGNITION])/s->timescales[T_HE_BURNING];
+ CLAMP(tau,0.0,1.0);
+
+
+#ifdef NEW_CORE_GROWTH
+ printf("OLD CORE %g\n",s->core_mass);
+ printf("New CORE CHEB Mc = %g vs M = %g, mcx = %g, mcbagb = %g\n",
+ s->core_mass,
+ s->mass,
+ mcx,
+ mcbagb);
+ fflush(stdout);
+ _exit(0);
+
+ double accretion_rate = Mdot_net(star);
+
+ /*
+ * Core growth rate according to normal
+ * (single-star) evolution
+ */
+ double dMc_dt_stellar_ev = 1e-6 * (mcbagb - mcx) / s->timescales[T_HE_BURNING];
+ double dMc_dt;
+
+ /*
+ * Stars with low-mass envelopes that are accreting likely
+ * formed their envelope by accretion. Limit the core growth
+ * rate to the accretion rate in this case, on the assumption
+ * the core grows as the amount of hydrogen available to it
+ * also grows.
+ */
+ if(accretion_rate > VERY_TINY &&
+ star->mass - star->core_mass < 0.05)
+ {
+ dMc_dt = MIN(accretion_rate, dMc_dt_stellar_ev);
+
+ /*
+ * The star's helium burning lifetime is actually
+ * the accretion lifetime. We don't know this, so set
+ * the timescale for nuclear burning to be long so it
+ * cannot limit the timestep.
+ */
+ s->timescales[T_HE_BURNING] = 1e10;
+ s->tn = 1e10;
+ }
+ else
+ {
+ /*
+ * No accretion : core grows at the normal stellar-evolution
+ * rate.
+ */
+ dMc_dt = dMc_dt_stellar_ev;
+ }
+
+ /*
+ * Hence core growth
+ */
+ s->core_mass += stardata->model.dtm * dMc_dt;
+
+ if(star->starnum==0)
+ printf("dMc/dt stellar ev wants %g, accretion gives %g : dMc/dt = %g, Mc = %20.10g\n",
+ dMc_dt_stellar_ev,
+ accretion_rate,
+ dMc_dt,
+ s->core_mass
+ );
+#else
s->core_mass = mcx + (mcbagb - mcx)*tau;
+#endif // NEW_CORE_GROWTH
+
+ s->GB_core_mass = 0.0;
+ s->CO_core_mass = 0.0;
Dprint("Core mass set %12.12e from mcx=%12.12e + bit=%12.12e (tau=%12.12e) : mass = %20.12e phase_start_mass = %20.12e\n",
s->core_mass,
@@ -123,8 +193,7 @@ Stellar_type stellar_structure_CHeB(struct new_stellar_structure_t *s,
,stardata->common.metallicity
#endif
);
-
-
+
if(rmin<rx)
{
taul = cbrt(log(rx/rmin));
diff --git a/src/stellar_structure/stellar_structure_EAGB.c b/src/stellar_structure/stellar_structure_EAGB.c
index bcbe7a9d0..79a4eaecf 100644
--- a/src/stellar_structure/stellar_structure_EAGB.c
+++ b/src/stellar_structure/stellar_structure_EAGB.c
@@ -69,7 +69,7 @@ Boolean stellar_structure_EAGB(struct new_stellar_structure_t * s,
A * (use == LOW ? D : B) *
pow(s->CO_core_mass, use == LOW ? p : q);
}
- printf("CO dmc/dt = %g\n",CO_dMc_dt);
+ //printf("CO dmc/dt = %g\n",CO_dMc_dt);
/*
* The helium core mass is always the helium core mass at the base of the AGB:
@@ -78,6 +78,11 @@ Boolean stellar_structure_EAGB(struct new_stellar_structure_t * s,
*/
s->core_mass = mcbagb;
+ /*
+ printf("CHEB > EAGB Mc %g McCO %g\n",
+ s->core_mass,
+ s->CO_core_mass);
+ */
Dprint("Set core from MAX(mcbagb=%12.12e,mcx=%12.12e) -> %g\n",
mcbagb,star->mcx_EAGB,s->core_mass);
@@ -139,7 +144,7 @@ Boolean stellar_structure_EAGB(struct new_stellar_structure_t * s,
if(count>10)
{
fflush(stdout);
- _exit(0);
+ //_exit(0);
}
#endif // NEW_CORE_GROWTH
diff --git a/src/stellar_structure/stellar_structure_HG.c b/src/stellar_structure/stellar_structure_HG.c
index 255dcae96..c5bd011a5 100644
--- a/src/stellar_structure/stellar_structure_HG.c
+++ b/src/stellar_structure/stellar_structure_HG.c
@@ -38,9 +38,6 @@ Stellar_type stellar_structure_HG(struct new_stellar_structure_t * s,
double mcbgbf = mcbgb(s,stardata);
s->core_mass = mcbgbf * ((1.0-tau)*rho+tau);
s->core_mass = MAX(s->core_mass,mcx);
- printf("MC HG %g (mcbgb = %g)\n",
- s->core_mass,
- mcbgbf);
#ifdef NELEMANS_MASS_PLOT
if(gijs_out==TRUE)
diff --git a/src/stellar_structure/stellar_structure_RG.c b/src/stellar_structure/stellar_structure_RG.c
index 5fb99fd56..3ec285246 100644
--- a/src/stellar_structure/stellar_structure_RG.c
+++ b/src/stellar_structure/stellar_structure_RG.c
@@ -19,7 +19,6 @@ Stellar_type stellar_structure_RG(struct new_stellar_structure_t *s,
*/
const double * metallicity_parameters=stardata->common.metallicity_parameters;
const double * giant_branch_parameters=stardata->common.giant_branch_parameters;
-
Stellar_type stellar_type=FIRST_GIANT_BRANCH;
/*
@@ -127,16 +126,14 @@ Stellar_type stellar_structure_RG(struct new_stellar_structure_t *s,
}
}
+ /*
double McBGB = mcbgb(s,stardata);
double LBGB = lbgbf(s->phase_start_mass,
stardata->common.giant_branch_parameters);
-
-
printf("Age %g, tBGB = %g, time along RGB %g\n",
s->age,
s->timescales[T_BGB],
s->age - s->timescales[T_BGB]);
-
printf("s->phase_start_core_mass = %g\n",s->core_mass);
printf("McBGB = %g\n",McBGB);
printf("s->core_mass = %g\n",s->core_mass);
@@ -146,6 +143,7 @@ Stellar_type stellar_structure_RG(struct new_stellar_structure_t *s,
s->phase_start_mass,
LBGB
);
+ */
/*
* Luminosity
@@ -156,14 +154,16 @@ Stellar_type stellar_structure_RG(struct new_stellar_structure_t *s,
/*
* Update the core masses
*/
+ /*
printf("Core masses were %g %g, growth rates %g %g\n",
s->GB_core_mass,
s->alt_core_mass,
GB_dMcdt,
dMcdt);
+ */
s->GB_core_mass += GB_dMcdt * stardata->model.dtm;
s->alt_core_mass += dMcdt * stardata->model.dtm;
- printf("Now %g %g\n",s->GB_core_mass,s->alt_core_mass);
+ //printf("Now %g %g\n",s->GB_core_mass,s->alt_core_mass);
/*
* Hence the luminosity as a function of core mass
@@ -172,10 +172,14 @@ Stellar_type stellar_structure_RG(struct new_stellar_structure_t *s,
B * pow(s->GB_core_mass,q),
D * pow(s->GB_core_mass,p)
);
-
+
+ /*
printf("NEW_CORE LGB BSE %g NEW %g\n",
s->luminosity,
s->alt_luminosity);
+
+ */
+
#endif // NEW_CORE_GROWTH
@@ -216,21 +220,16 @@ Stellar_type stellar_structure_RG(struct new_stellar_structure_t *s,
}
#ifdef NEW_CORE_GROWTH
-
printf("NEW_CORE Mc BSE %g alt %g GB_core_mass %g\n",
s->core_mass,
s->alt_core_mass,
s->GB_core_mass);
-
fflush(NULL);
-
static int count = 0;
count++;
//if(count>1)_exit(0);
#endif
-
-
Dprint("setting r from rgbf to ");
s->radius = rgbf(s->mass,
s->luminosity,
diff --git a/src/timestep/timestep_CHeB.c b/src/timestep/timestep_CHeB.c
index 879949194..603578b0b 100644
--- a/src/timestep/timestep_CHeB.c
+++ b/src/timestep/timestep_CHeB.c
@@ -11,7 +11,30 @@ void timestep_CHeB(Timestep_prototype_args)
#endif // HRDIAG_DTLOG
HRDIAG_DT_NEGATIVE_CHECK;
#endif // HRDIAG
- *time_remaining = MIN(tn,
- (tscls[T_HE_BURNING] + tscls[T_HE_IGNITION]))
- - age;
+
+ printf("CF tn = %g vs (%g + %g) vs age = %g\n",
+ tn,
+ tscls[T_HE_BURNING],
+ tscls[T_HE_IGNITION],
+ age);
+
+ double accretion_rate = Mdot_net(star);
+ if(IS_ZERO(accretion_rate))
+ {
+ *time_remaining =
+ MIN(tn,
+ (tscls[T_HE_BURNING] + tscls[T_HE_IGNITION]))
+ - age;
+ }
+ else
+ {
+ /*
+ * Accretion CHeB star has an unknown lifetime,
+ * because it depends on the accretion timescale,
+ * which is generally unknown. Assume the mass
+ * transfer rate timescale limiter will resolve
+ * the timescale.
+ */
+ *time_remaining = 1e10;
+ }
}
diff --git a/src/zfuncs/lgbtf.c b/src/zfuncs/lgbtf.c
index e5329293f..79dc9b007 100644
--- a/src/zfuncs/lgbtf.c
+++ b/src/zfuncs/lgbtf.c
@@ -16,34 +16,13 @@ double lgbtf(const double t,
*/
double result,x;
- printf("LGBTF : t=%g A=%g D=%g p=%g tinf1=%g tinf2=%g tx=%g \n",
- t,A,D,p,tinf1,tinf2,tx);
-
- double dummyMc = mcgbtf(t,
- A,
- GB,
- tinf1,
- tinf2,
- tx);
- printf("Dummy Mc %g\n",dummyMc);
-
if(LESS_OR_EQUAL(t,tx))
{
result = D * pow((p-1) * A * D * (tinf1 - t),
-p/(p-1));
- printf("use low : L=%g : tinf1-t=%g exponent = %g mult = %g\n",
- result,
- tinf1 - t,
- p,
- D
- );
- printf("equiv Mc = %g\n",
- pow(A*D*(p-1)*(tinf1-t),1.0/(1.0-p)));
-
}
else
{
- printf("use high\n");
x = q - 1.0;
double y = tinf2-t;
//y=MAX(1e-2,y); // required any more?
@@ -51,7 +30,6 @@ double lgbtf(const double t,
-q/x);
}
- printf("LGBTF gives L = %g\n",result);
return (result);
}
diff --git a/tbse b/tbse
index fe361c303..fe039c208 100755
--- a/tbse
+++ b/tbse
@@ -57,7 +57,7 @@ STELLAR_TYPE1=1
STELLAR_TYPE2=1
# stellar masses (solar units)
-M1=1
+M1=20
M2=0.1
# orbit:
@@ -85,7 +85,6 @@ METALLICITY=0.02
#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
@@ -94,7 +93,7 @@ METALLICITY=0.02
#M1=3; M2=1; ORBITAL_PERIOD=30
# AGB test
-#M1=4; M2=1; ORBITAL_PERIOD=1000000
+M1=4; M2=1; ORBITAL_PERIOD=4000
# CBDISC
#M1=2; M2=1; ORBITAL_PERIOD=1000; MAX_EVOLUTION_TIME=2000; METALLICITY=0.014
@@ -355,7 +354,7 @@ MINIMUM_ORBITAL_PERIOD_FOR_INSTANT_RLOF=0
COMENV_PRESCRIPTION=0
# Alpha (energy from orbit put into envelope, 1.0)
-ALPHA_CE=10.0
+ALPHA_CE=0.2
# Lambda: structure parameter
# if >0 (e.g. 0.5) then use the given value
@@ -407,7 +406,7 @@ COMENV_DISC_MASS_FRACTION=1e-2
# -2 : use comenv specific angular momentum (minimum?)
# -3 : use L2 specific angular momentum (maximum?)
#
-COMENV_DISC_ANGMOM_FRACTION=-3
+COMENV_DISC_ANGMOM_FRACTION=0.01
# Circumbinary disc physical properties
CBDISC_GAMMA=1.6666666 # 1.6666666
@@ -444,7 +443,7 @@ CBDISC_MASS_LOSS_ISM_RAM_PRESSURE_MULTIPLIER=1.0
# multipliers for X-ray and FUV winds
CBDISC_MASS_LOSS_FUV_MULTIPLIER=0.0
-CBDISC_MASS_LOSS_XRAY_MULTIPLIER=1.0
+CBDISC_MASS_LOSS_XRAY_MULTIPLIER=1.0e-2
# couple photoevaporation and viscous losses at the inner edge
CBDISC_VISCOUS_PHOTOEVAPORATION_COUPLING=0
diff --git a/test_new.pl b/test_new.pl
index 2efb56454..1e37df5cd 100755
--- a/test_new.pl
+++ b/test_new.pl
@@ -4,11 +4,23 @@ use strict;
my $altcode = '/home/rgi/svn/binary_c/releases/2.0pre27';
+#
+# Test two versions of binary_c
+# to see if their output is the same.
+#
+# To make it work, compile binary_c to
+# two executables.
+#
+# The "reference" version should be called binary_c-ref
+# The "new" version should be called binary_c
+#
+
+
while(1)
{
# run a random system with this code
print "Running reference system\n";
- my $r = `tbse --random_systems 1 --log_filename /tmp/c_log.ref 2>\&1`;
+ my $r = `env BIN_C=binary_c-ref tbse --random_systems 1 --log_filename /tmp/c_log.ref 2>\&1`;
my $args = ($r=~/^(--M_1.*)/)[0];
print "ARGS:\n\n$args\n\n";
--
GitLab