diff --git a/configure b/configure
index 52ea70cd918f4fcb850615cfa0ac81be07920a3a..f21bf2e630b92a34df521e7002cc5d818c111e5d 100755
--- a/configure
+++ b/configure
@@ -1823,7 +1823,7 @@ sub features
{
check => $gprof,
perl_pre => "
- \$cflags .= ' -pg -Og ';
+ \$cflags .= ' -pg -Og -g ';
",
},
diff --git a/src/binary_c_structures.h b/src/binary_c_structures.h
index 50fc4d185d3cd90659f29f4a6afcec658827f6aa..89e8cc0d64403c85777a627b504d7842c71c1001 100644
--- a/src/binary_c_structures.h
+++ b/src/binary_c_structures.h
@@ -1118,7 +1118,6 @@ struct disc_loss_t {
};
struct disc_t {
- struct stardata_t * stardata;
struct memoize_hash_t * memo; /* memoize data */
double M; // disc total mass
double alpha; // disc viscosity constant (dimensionless)
diff --git a/src/common_envelope/common_envelope_evolution_BSE.c b/src/common_envelope/common_envelope_evolution_BSE.c
index 25837166c8a84edc0b1fba16debd091f1edabdce..5bd2e5d511d4da6ab91c21d1e0dc5d7e8b5843da 100644
--- a/src/common_envelope/common_envelope_evolution_BSE.c
+++ b/src/common_envelope/common_envelope_evolution_BSE.c
@@ -2005,7 +2005,6 @@ end_section:
cbdisc = new_disc(DISC_CIRCUMBINARY,
stardata,
DISC_OUTER_EDGE);
- cbdisc->stardata = stardata;
double discM = mass_ejected * f_M;
double discJ = J_ejected * f_J;
diff --git a/src/disc/disc_evaporate_cbdiscs.c b/src/disc/disc_evaporate_cbdiscs.c
index 4efa4b22ae028d88e5a130e105a6b857c9445c88..cbaf4cd6b2859a0c16c3f350ae5ba6797b5613e8 100644
--- a/src/disc/disc_evaporate_cbdiscs.c
+++ b/src/disc/disc_evaporate_cbdiscs.c
@@ -24,6 +24,7 @@ void disc_evaporate_cbdiscs(struct stardata_t * const stardata)
disc->M = 0.0;
disc->J = 0.0;
disc->F = 0.0;
+ disc->end_count = disc->end_count==0 ? 1 : disc->end_count;
}
}
}
diff --git a/src/disc/disc_failure_mode.c b/src/disc/disc_failure_mode.c
index 83d60c8a77604dab0cd2fc667668c52503a684cf..a3e75f8a581a4313d6d25349b49ac373043c6ea9 100644
--- a/src/disc/disc_failure_mode.c
+++ b/src/disc/disc_failure_mode.c
@@ -27,7 +27,7 @@ char * disc_failure_mode(struct stardata_t * const stardata,
* Special case : flux is negative,
* disc should be evaporated
*/
- Boolean flux_negative =
+ const Boolean flux_negative =
status == CONVERGENCE_FLUX_NEGATIVE ?
TRUE : FALSE;
@@ -35,19 +35,19 @@ char * disc_failure_mode(struct stardata_t * const stardata,
* Reasons for failure...
*/
const double Ldisc = disc_total_luminosity(disc);
- Boolean too_dim =
+ const Boolean too_dim =
(stardata->preferences->cbdisc_minimum_luminosity >-TINY ?
(Ldisc < L_SUN * stardata->preferences->cbdisc_minimum_luminosity) :
(Ldisc/binary->L < stardata->preferences->cbdisc_minimum_luminosity));
- Boolean too_low_mass =
+ const Boolean too_low_mass =
disc->M < stardata->preferences->cbdisc_minimum_mass * M_SUN;
- Boolean too_ringlike =
+ const Boolean too_ringlike =
Is_not_zero(disc->fRing) &&
(disc->fRing < stardata->preferences->cbdisc_minimum_fRing);
- Boolean too_old =
+ const Boolean too_old =
disc->lifetime > DISC_MAX_LIFETIME ||
(disc->type == DISC_CIRCUMBINARY &&
Is_not_zero(stardata->preferences->cbdisc_max_lifetime) &&
@@ -64,22 +64,24 @@ char * disc_failure_mode(struct stardata_t * const stardata,
too_ringlike ||
disc_is_ring ||
too_old)
- {
-#define __STRING_N 155
- reason = Malloc(sizeof(char)*__STRING_N);
- snprintf(reason,(__STRING_N-1),"%s%s%s%s%s%s : disc->M=%6.3e Msun (%6.3e), disc->J=%6.3e (%6.3e), disc->F=%6.3e",
- (too_low_mass ? "too low mass" : ""), // 12
- (too_dim ? " too dim" : ""), // + 8 = 20
- (disc_is_ring ? " is ring" : ""), // + 8 = 28
- (flux_negative ? " flux negative" : ""), // + 14 = 42
- (too_ringlike ? " too ringlike" : ""), // + 13 = 55
- (too_old ? " too old " : ""), // + 8 = 63
- disc->M/M_SUN, //
- disc_total_mass(disc)/M_SUN,
- disc->J,
- disc_total_angular_momentum(disc,binary),
- disc->F
- );
+ {
+ if(asprintf(&reason,
+ "%s%s%s%s%s%s : disc->M=%6.3e Msun (%6.3e), disc->J=%6.3e (%6.3e), disc->F=%6.3e",
+ (too_low_mass ? "too low mass" : ""), // 12
+ (too_dim ? " too dim" : ""), // + 8 = 20
+ (disc_is_ring ? " is ring" : ""), // + 8 = 28
+ (flux_negative ? " flux negative" : ""), // + 14 = 42
+ (too_ringlike ? " too ringlike" : ""), // + 13 = 55
+ (too_old ? " too old " : ""), // + 8 = 63
+ disc->M/M_SUN, //
+ disc_total_mass(disc)/M_SUN,
+ disc->J,
+ disc_total_angular_momentum(disc,binary),
+ disc->F
+ )<0)
+ {
+ reason = NULL;
+ }
}
else
{
diff --git a/src/disc/disc_hottest_zone_at_radius.c b/src/disc/disc_hottest_zone_at_radius.c
index ea141b7698267b975053d4183db937dd1eb30282..f7e8745add8a92b4c3e6c3ef985fb72840edc499 100644
--- a/src/disc/disc_hottest_zone_at_radius.c
+++ b/src/disc/disc_hottest_zone_at_radius.c
@@ -2,35 +2,45 @@
#ifdef DISCS
-Disc_zone_counter disc_hottest_zone_at_radius(struct disc_thermal_zone_t * zones,
- const double radius,
- const Boolean all,
- double * Thottest_p)
+int disc_hottest_zone_at_radius(struct disc_thermal_zone_t * const zones,
+ const double radius,
+ const Boolean all,
+ double * Thottest_p)
{
/*
* Return the counter of the hottest zone at
* given radius, and set *Thottest_p to this temperature
* if Thottest_p != NULL.
*
- * Returns -1 on failure.
+ * Returns -1 on failure (note that the return type
+ * is int, not Disc_zone_counter which is an
+ * unsigned int).
*
* Invalid zones are skipped, unless all == TRUE.
*/
double Thottest = 0.0;
- Disc_zone_counter i,ihottest = -1;
-
+ Disc_zone_counter i;
+ int ihottest = -1;
+
for(i=0;i<DISCS_MAX_N_ZONES;i++)
{
if(all == TRUE ||
zones[i].valid == TRUE)
{
const double T = power_law(&zones[i].Tlaw,radius);
+ //printf("T[%u](%g) = %g\n",i,radius,T);
if(T > Thottest)
{
ihottest = i;
Thottest = T;
}
}
+ /*
+ else
+ {
+ printf("Zkip zone %u\n",i);
+ }
+ */
}
if(ihottest!=-1 && Thottest_p != NULL)
{
diff --git a/src/disc/disc_logging.c b/src/disc/disc_logging.c
index f6a3324fa84aec49c2c3662e8590cc055b029a24..7e42a1870192dc3862dfb5d328204ed07761899b 100644
--- a/src/disc/disc_logging.c
+++ b/src/disc/disc_logging.c
@@ -10,8 +10,9 @@
#include "valgrind/valgrind.h"
#endif
#include "disc_constraints.h"
-
-static void close_disc_logs(struct stardata_t * stardata);
+static void disc_log_end(struct stardata_t * const stardata,
+ struct disc_t * const disc);
+static void close_disc_logs(struct stardata_t * const stardata);
/*
* Perform logging for the disc.
@@ -32,9 +33,17 @@ void disc_logging(struct stardata_t * const stardata,
int i;
for(i=0; i<stardata->common.ndiscs; i++)
{
- struct disc_t * disc = & stardata->common.discs[i];
+ struct disc_t * const disc = & stardata->common.discs[i];
- if(disc->converged == TRUE)
+ if(Is_really_zero(disc->M) &&
+ disc->lifetime > 0.0)
+ {
+ /*
+ * Disc did exist, now doesn't. End its log.
+ */
+ disc_log_end(stardata,disc);
+ }
+ else if(disc->converged == TRUE)
{
Boolean do_log_output = FALSE;
@@ -246,38 +255,24 @@ void disc_logging(struct stardata_t * const stardata,
/*
* If the disc has evaporated, log its end.
*/
- if(Is_really_zero(disc->M)&&
+ if(Is_really_zero(disc->M) &&
disc->end_count==0)
{
- disc->end_count++;
- Printf("DISC_END %30.20g %g %g %g %d %d %g %g %g %g %g\n",
- stardata->model.time,
- stardata->model.probability,
- disc->lifetime/YEAR_LENGTH_IN_SECONDS,
- stardata->common.orbit.eccentricity,
- stardata->star[0].stellar_type,
- stardata->star[1].stellar_type,
- stardata->star[0].pms_mass,
- stardata->star[1].pms_mass,
- stardata->common.zams_period,
- stardata->common.zams_separation,
- stardata->common.zams_eccentricity
- );
- //Backtrace;
-
+ disc_log_end(stardata,disc);
+ }
+
#ifdef __HAVE_VALGRIND__
- /*
- * Discs can greatly slow evolution, and when
- * running with valgrind cause timeouts.
- * Reset the timeout here.
- */
- if(RUNNING_ON_VALGRIND)
- {
- reset_binary_c_timeout();
- }
+ /*
+ * Discs can greatly slow evolution, and when
+ * running with valgrind cause timeouts.
+ * Reset the timeout here.
+ */
+ if(RUNNING_ON_VALGRIND)
+ {
+ reset_binary_c_timeout();
+ }
#endif // __HAVE_VALGRIND__
- }
#endif // DISC_LOG_POPSYN
@@ -829,11 +824,19 @@ void disc_logging(struct stardata_t * const stardata,
}
#endif // DISC_LOG
+
+ if(disc->lifetime > 0.0 &&
+ Is_really_zero(disc->M) &&
+ disc->end_count==0)
+ {
+ disc_log_end(stardata,disc);
+ }
+
if(!Disc_is_disc(disc))
{
close_disc_logs(stardata);
}
-
+
/*
* Free memory if required
*/
@@ -843,14 +846,19 @@ void disc_logging(struct stardata_t * const stardata,
}
}
}
- else
+
+ /*
+ * Force a close of the disc log
+ */
+ if(disc->converged == FALSE ||
+ Is_really_zero(disc->M))
{
close_disc_logs(stardata);
}
}
}
-static void close_disc_logs(struct stardata_t * stardata)
+static void close_disc_logs(struct stardata_t * const stardata)
{
#ifdef DISC_LOG_2D
if(stardata->preferences->disc_log2d!= DISC_LOG_LEVEL_NONE &&
@@ -870,5 +878,25 @@ static void close_disc_logs(struct stardata_t * stardata)
#endif //DISC_LOG
}
+static void disc_log_end(struct stardata_t * const stardata,
+ struct disc_t * const disc)
+{
+ disc->end_count++;
+ Printf("DISC_END %30.20g %g %g %g %d %d %g %g %g %g %g\n",
+ stardata->model.time,
+ stardata->model.probability,
+ disc->lifetime/YEAR_LENGTH_IN_SECONDS,
+ stardata->common.orbit.eccentricity,
+ stardata->star[0].stellar_type,
+ stardata->star[1].stellar_type,
+ stardata->star[0].pms_mass,
+ stardata->star[1].pms_mass,
+ stardata->common.zams_period,
+ stardata->common.zams_separation,
+ stardata->common.zams_eccentricity
+ );
+ //Backtrace;
+
+}
#endif // DISCS and any DISC logging
diff --git a/src/disc/disc_macros.h b/src/disc/disc_macros.h
index 670ba75a37589acd3b0300f79cb5d8273d339b2f..644ae65dfa375c97f9b4a5887fe79be740b677c8 100644
--- a/src/disc/disc_macros.h
+++ b/src/disc/disc_macros.h
@@ -206,5 +206,16 @@
#define DISC_STRIPPING_TIMESCALE_ORBIT 3
+/*
+ * Macro to determine whether "const disc" should be used
+ * in function calls, or "disc". If we're using MEMOIZE,
+ * disc can be changed (because it contains disc->memo which
+ * is writable) so we don't want the const.
+ */
+#ifdef MEMOIZE
+#define DISC_CONST
+#else
+#define DISC_CONST const
+#endif
#endif // DISC_MACROS_H
diff --git a/src/disc/disc_nearest_zone_n.c b/src/disc/disc_nearest_zone_n.c
index 3704261bcea16c2d266a3f9fb2042d8da1b40d8c..df86b650cd219cfbf446f645a35ae489a263cd72 100644
--- a/src/disc/disc_nearest_zone_n.c
+++ b/src/disc/disc_nearest_zone_n.c
@@ -3,7 +3,7 @@
#ifdef DISCS
Disc_zone_counter disc_nearest_zone_n(const double radius,
- const struct disc_t * const disc)
+ struct disc_t * const disc)
{
/*
* Return the nearest zone to the given radius
diff --git a/src/disc/disc_new_zone_list.c b/src/disc/disc_new_zone_list.c
index a19c82e3c202836cb8bd75f54b06ba6e2070cdfa..482f0ad58d66e2f896ee80b1c4296058f33d9e4e 100644
--- a/src/disc/disc_new_zone_list.c
+++ b/src/disc/disc_new_zone_list.c
@@ -75,7 +75,7 @@ Boolean disc_new_zone_list(struct disc_t * const disc,
DISC_REFERENCE_LENGTH, /* R0=1 */
0.0 /* R1=0 is ignored */);
radiative_out->type = DISC_ZONE_TYPE_RADIATIVE_OUTER;
-
+
/*
* Set up the inner radiative zone
*/
@@ -84,13 +84,14 @@ Boolean disc_new_zone_list(struct disc_t * const disc,
Exponent_radiative_in,
DISC_REFERENCE_LENGTH, /* R0=1 */
0.0 /* R1=0 is ignored */);
- radiative_in->type = DISC_ZONE_TYPE_RADIATIVE_INNER;
-
+ radiative_in->type = DISC_ZONE_TYPE_RADIATIVE_INNER;
+
/*
* Set sigma0 from Tvisc0 (which is a free parameter)
*/
- disc->sigma0 = pow(disc->Tvisc0,5.0/2.0) / sqrt(Prefactor_viscous_fraction);
-
+ disc->sigma0 = pow(disc->Tvisc0,5.0/2.0)
+ / sqrt(Prefactor_viscous_fraction);
+
/*
* And hence update the disc angular momentum flux, which
* depends only on sigma0 and some constants of the disc (alpha, mu,
diff --git a/src/disc/disc_new_zone_radii.c b/src/disc/disc_new_zone_radii.c
index c6e000d10848029774659e6673af83063d7f63b2..a859ef1158f6a7edfe129e5dd2f83866373902b7 100644
--- a/src/disc/disc_new_zone_radii.c
+++ b/src/disc/disc_new_zone_radii.c
@@ -2,6 +2,9 @@
#ifdef DISCS
+/*
+ * Zprintf is the local debugging function
+ */
#undef Zprintf
#define Zprintf(...) /* */
//#define Zprintf(...) printf(__VA_ARGS__);fflush(NULL);
@@ -86,8 +89,9 @@ Boolean disc_new_zone_radii(struct disc_t * const disc)
Disc_zone_counter zones_remaining = DISCS_MAX_N_ZONES;
Disc_zone_counter zones_allocated = 0;
+//#define ___DEBUGGING
#ifdef ___DEBUGGING
- if(0)
+ if(1)
{
/*
* Debugging:
@@ -98,13 +102,16 @@ Boolean disc_new_zone_radii(struct disc_t * const disc)
{
for(j=0;j<=i;j++)
{
- const double Rx = disc_Rcross(&P[i],
- &P[j]);
- printf("%d x %d : %g %s\n",
- i,
- j,
- Rx/R_SUN,
- R_SOLAR);
+ if(i!=j)
+ {
+ const double Rx = disc_Rcross(&P[i],
+ &P[j]);
+ printf("%u x %u : %g %s\n",
+ i,
+ j,
+ Rx/R_SUN,
+ R_SOLAR);
+ }
}
}
}
@@ -120,27 +127,27 @@ Boolean disc_new_zone_radii(struct disc_t * const disc)
*
* Zones that have already been used are ignored.
*/
- Disc_zone_counter ihottest = disc_hottest_zone_at_radius(P,r,FALSE,NULL);
+ int ihottest =
+ disc_hottest_zone_at_radius(P,r,FALSE,NULL);
- Zprintf("Hottest at r = %30.15e %s is zone %u\n",
+ Zprintf("Hottest at r = %30.15e %s is zone %d (disc Rin %g, Rout %g)\n",
r/R_SUN,
R_SOLAR,
- ihottest);
+ ihottest,
+ disc->Rin,
+ disc->Rout);
/*
* If ihottest is -1 we're in trouble:
- * this should not happen!
+ * this should not happen, but it can when
+ * one of the power laws has A=0 or something
+ * like that. Such solutions should be rejected,
+ * so return FALSE.
*/
if(ihottest == -1)
{
Safe_free(P);
Safe_free(Q);
- struct stardata_t * const stardata = disc->stardata;
- char * s = cmd_line_argstring(stardata);
- Exit_binary_c(2,
- "args were : %s\n",
- cmd_line_argstring(stardata));
- Safe_free(s);
return FALSE;
}
@@ -151,12 +158,12 @@ Boolean disc_new_zone_radii(struct disc_t * const disc)
* its next position, and set P[ihottest] to be invalid
* so it cannot be set into Q twice.
*/
- Zprintf("Set zone P[%u] in Q[%u]\n",ihottest,zones_allocated);
+ Zprintf("Set zone P[%d] in Q[%u]\n",ihottest,zones_allocated);
Copy_zone(&P[ihottest],
&Q[zones_allocated]);
P[ihottest].valid = FALSE;
zones_remaining--;
-
+
/*
* Set the end of this new zone, and r,
* to the next zone boundary. This is where
@@ -177,7 +184,7 @@ Boolean disc_new_zone_radii(struct disc_t * const disc)
if(i != ihottest &&
P[i].valid == TRUE)
{
- Zprintf("Call RCross %u (%g %g %g %s) %u (%g %g %g %s)\n",
+ Zprintf("Call RCross %d (%g %g %g %s) %u (%g %g %g %s)\n",
ihottest,
P[ihottest].Tlaw.A0,
P[ihottest].Tlaw.exponent,
diff --git a/src/disc/disc_parameters.h b/src/disc/disc_parameters.h
index 128fa9ccaf3b589624f1dda60f03750d8d96df58..f53dc57c62447300392021fded1de5361c0a502f 100644
--- a/src/disc/disc_parameters.h
+++ b/src/disc/disc_parameters.h
@@ -9,7 +9,7 @@
/*
* Logging for individual systems with discs
*/
-//#define DISC_LOG
+#define DISC_LOG
//#define DISC_LOG_2D
/*
diff --git a/src/disc/disc_prototypes.h b/src/disc/disc_prototypes.h
index 1c7b06c81d3cd9cc52f787d9a2f0249a2383a9e5..9a7ee6f335ddf992f762d6661c2f6c5dc8fa4907 100644
--- a/src/disc/disc_prototypes.h
+++ b/src/disc/disc_prototypes.h
@@ -65,7 +65,7 @@ double disc_temperature(const double radius,
double Pure_function disc_inner_edge_temperature(const struct disc_t * const disc);
double Pure_function disc_outer_edge_temperature(const struct disc_t * const disc);
double disc_column_density(const double radius,
- const struct disc_t * const disc);
+ DISC_CONST struct disc_t * const disc);
double zone_mass_numerical(struct disc_t * disc,
const struct disc_thermal_zone_t * zone);
double disc_specific_angular_momentum(const double radius,
@@ -77,7 +77,7 @@ double disc_injection_radius(const double h,
double disc_velocity(const double radius,
const struct binary_system_t * const binary);
double disc_gravitational_pressure(const double radius,
- const struct disc_t * const disc,
+ DISC_CONST struct disc_t * const disc,
const struct binary_system_t * const binary);
double disc_pressure_radius(struct disc_t * const disc,
const struct binary_system_t * const binary,
@@ -154,13 +154,15 @@ double disc_viscous_timescale(const double radius,
double disc_mass_weighted_viscous_timescale(struct disc_t * const disc,
const struct binary_system_t * const binary);
Disc_zone_counter disc_zone_n(const double radius,
- const struct disc_t * const disc);
+ struct disc_t * RESTRICT const disc);
+Disc_zone_counter disc_zone_n_memoize(const double radius,
+ struct disc_t * RESTRICT const disc);
double generic_power_law(double radius,
- const struct disc_t * RESTRICT const disc,
+ struct disc_t * RESTRICT const disc,
const int n);
Disc_zone_counter disc_nearest_zone_n(const double radius,
- const struct disc_t * const disc);
+ struct disc_t * const disc);
void cbdisc_eccentricity_pumping_rate(struct stardata_t * stardata,
struct disc_t * disc,
@@ -494,10 +496,10 @@ Boolean disc_monotonic_check(struct stardata_t * const stardata,
);
-Disc_zone_counter disc_hottest_zone_at_radius(struct disc_thermal_zone_t * zones,
- const double radius,
- const Boolean all,
- double * Thottest_p);
+int disc_hottest_zone_at_radius(struct disc_thermal_zone_t * const zones,
+ const double radius,
+ const Boolean all,
+ double * Thottest_p);
#endif // DISCS
#endif // DISC_PROTOTYPES_H
diff --git a/src/disc/disc_structure_functions.c b/src/disc/disc_structure_functions.c
index 003af166ac58604c7d160e62244111ae7d811d5a..5fe13ba10c823d8fa40fe975391f27674a134d0c 100644
--- a/src/disc/disc_structure_functions.c
+++ b/src/disc/disc_structure_functions.c
@@ -49,7 +49,7 @@ double Pure_function disc_outer_edge_temperature(const struct disc_t * const dis
}
double generic_power_law(double radius,
- const struct disc_t * RESTRICT const disc,
+ struct disc_t * RESTRICT const disc,
const int n)
{
/*
@@ -84,7 +84,7 @@ double generic_power_law(double radius,
double disc_column_density(const double radius,
- const struct disc_t * const disc)
+ DISC_CONST struct disc_t * const disc)
{
/*
* Column density as a function of radius
@@ -134,7 +134,7 @@ double Pure_function disc_velocity(const double radius,
double disc_gravitational_pressure(const double radius,
- const struct disc_t * const disc,
+ DISC_CONST struct disc_t * const disc,
const struct binary_system_t * const binary)
{
/*
diff --git a/src/disc/disc_thermal_zones.h b/src/disc/disc_thermal_zones.h
index d6a4207a0e96c4f83f5c5187b8867d8f19bd29a3..18e664b153600a5dbdf1481085f9356020c48f26 100644
--- a/src/disc/disc_thermal_zones.h
+++ b/src/disc/disc_thermal_zones.h
@@ -17,8 +17,11 @@
*/
//#define Disc_zone_n(R,D) ((D)->n_thermal_zones==1 ? 0 : (R) < (D)->thermal_zones[1].rstart ? 0 : 1)
-// was 2.38s
+#ifdef MEMOIZE
+#define Disc_zone_n(R,D) (disc_zone_n_memoize((R),(D)))
+#else//MEMOIZE
#define Disc_zone_n(R,D) (disc_zone_n((R),(D)))
+#endif//MEMOIZE
#define Disc_nearest_zone_n(R,D) (disc_nearest_zone_n((R),(D)))
diff --git a/src/disc/disc_zone_n.c b/src/disc/disc_zone_n.c
index 8bf29f19d64632c49cb18ecc67fe18b9faa95b28..47b49e00bb0806ef642b21339c4f110c186fdd2f 100644
--- a/src/disc/disc_zone_n.c
+++ b/src/disc/disc_zone_n.c
@@ -2,8 +2,27 @@
#ifdef DISCS
+Disc_zone_counter disc_zone_n_memoize(const double radius,
+ struct disc_t * RESTRICT const disc)
+{
+ return
+ (Disc_zone_counter)
+ Memoize(
+ disc->memo,
+ "disc_zone_n",
+ 1,
+ scalar,
+ Disc_zone_counter,
+ 1,
+ radius,
+ const double,
+ 1,
+ disc_zone_n(radius,disc)
+ );
+}
+
Disc_zone_counter disc_zone_n(const double radius,
- const struct disc_t * const disc)
+ struct disc_t * RESTRICT const disc)
{
/*
* Given a disc structure and a radial coordinate (radius)
@@ -14,6 +33,9 @@ Disc_zone_counter disc_zone_n(const double radius,
* This prevents floating point errors.
*/
+
+
+
#if DISC_DEBUG >= 3
double d = radius-Disc_inner_zone(disc)->rstart;
Discdebug(3,
@@ -38,10 +60,10 @@ Disc_zone_counter disc_zone_n(const double radius,
*/
#define ZONE_FLUFF (TINY)
-#define Is_not_our_zone(RADIUS,ZONE) ( \
- (ZONE)->valid == FALSE || \
- (RADIUS) < (ZONE)->rstart*(1.0 - ZONE_FLUFF) || \
- (RADIUS) > (ZONE)->rend*(1.0 + ZONE_FLUFF) \
+#define Is_not_our_zone(RADIUS,ZONE) ( \
+ unlikely((ZONE)->valid == FALSE) || \
+ (RADIUS) < (ZONE)->rstart*(1.0 - ZONE_FLUFF) || \
+ (RADIUS) > (ZONE)->rend*(1.0 + ZONE_FLUFF) \
)
#define Is_our_zone(RADIUS,ZONE) \
@@ -55,14 +77,14 @@ Disc_zone_counter disc_zone_n(const double radius,
/* compare to inner and outer edges */
else if(unlikely(Float_same_within_eps(radius/Disc_outer_zone(disc)->rend,
- 1.0,
- DISC_EDGE_EPS)))
+ 1.0,
+ DISC_EDGE_EPS)))
{
return Disc_outer_zone_n(disc);
}
else if(unlikely(Float_same_within_eps(radius/Disc_inner_zone(disc)->rstart,
- 1.0,
- DISC_EDGE_EPS)))
+ 1.0,
+ DISC_EDGE_EPS)))
{
return Disc_inner_zone_n(disc);
}
@@ -108,7 +130,7 @@ Disc_zone_counter disc_zone_n(const double radius,
for(;i<DISCS_MAX_N_ZONES;i++,z++)
{
- if(! Is_not_our_zone(radius,z))
+ if(Is_our_zone(radius,z))
{
/* all conditions ok : this is our zone! */
return i;
diff --git a/src/perl/scripts2/cbdiscs.pl b/src/perl/scripts2/cbdiscs.pl
index 1dd1230bc5ecb74956f206bd0ab736586b24195d..5e683ffbdbaafa4330e0fcc4081a6625b25eb2c9 100755
--- a/src/perl/scripts2/cbdiscs.pl
+++ b/src/perl/scripts2/cbdiscs.pl
@@ -431,7 +431,7 @@ sub parse_data
# we want orbital period in days
my $orbital_period_binned = $population->rebin(safelog10($orbital_period*365.25),0.1);
my $eccentricity_binned = $population->rebin($eccentricity,0.025);
- my $logLstar_binned = $population->rebin(safelog10($Lstar),0.1);
+ my $logLstar_binned = $population->rebin(safelog10($Lstar),0.025);
# 1D distributions
$r->{Mdot_viscous}->{$population->rebin(safelog10($Mdot_viscous),0.1)} += $dtp;
@@ -441,34 +441,34 @@ sub parse_data
$r->{Mdot_FUV}->{$population->rebin(safelog10($Mdot_FUV),0.1)} += $dtp;
$r->{Mdot_Xray}->{$population->rebin(safelog10($Mdot_Xray),0.1)} += $dtp;
$r->{logMdisc}->{$population->rebin(safelog10($Mdisc),0.1)} += $dtp;
- $r->{logMbinary}->{$population->rebin(safelog10($Mbinary),0.1)} += $dtp;
- $r->{logLdisc}->{$population->rebin(safelog10($Ldisc),0.1)} += $dtp;
+ $r->{logMbinary}->{$population->rebin(safelog10($Mbinary),0.025)} += $dtp;
+ $r->{logLdisc}->{$population->rebin(safelog10($Ldisc),0.025)} += $dtp;
$r->{logJdisc}->{$population->rebin(safelog10($Jdisc),0.1)} += $dtp;
- $r->{logJbinary}->{$population->rebin(safelog10($Jbinary),0.1)} += $dtp;
+ $r->{logJbinary}->{$population->rebin(safelog10($Jbinary),0.025)} += $dtp;
$r->{age}->{$population->rebin(safelog10($t),0.1)} += $dtp;
- $r->{separation}->{$population->rebin(safelog10($separation),0.1)} += $dtp;
+ $r->{separation}->{$population->rebin(safelog10($separation),0.025)} += $dtp;
$r->{orbital_period}->{$orbital_period_binned} += $dtp;
$r->{eccentricity}->{$eccentricity_binned} += $dtp;
$r->{edot_resonance}->{$population->rebin($edot_resonance,0.1)} += $dtp;
$r->{stellar_type0}->{$stellar_type0} += $dtp;
$r->{stellar_type1}->{$stellar_type1} += $dtp;
- $r->{logRin}->{$population->rebin(safelog10($Rin),0.1)} += $dtp;
- $r->{logRout}->{$population->rebin(safelog10($Rout),0.1)} += $dtp;
- $r->{logRout_over_Rin}->{$population->rebin(safelog10($Rout/MAX(1e-50,$Rin)),0.1)} += $dtp;
- $r->{logRin_over_separation}->{$population->rebin(safelog10($Rin/MAX(1e-50,$separation)),0.1)} += $dtp;
- $r->{logRout_over_separation}->{$population->rebin(safelog10($Rout/MAX(1e-50,$separation)),0.1)} += $dtp;
- $r->{logHRin}->{$population->rebin(safelog10($HRin),0.1)} += $dtp;
- $r->{logHRout}->{$population->rebin(safelog10($HRout),0.1)} += $dtp;
+ $r->{logRin}->{$population->rebin(safelog10($Rin),0.025)} += $dtp;
+ $r->{logRout}->{$population->rebin(safelog10($Rout),0.025)} += $dtp;
+ $r->{logRout_over_Rin}->{$population->rebin(safelog10($Rout/MAX(1e-50,$Rin)),0.025)} += $dtp;
+ $r->{logRin_over_separation}->{$population->rebin(safelog10($Rin/MAX(1e-50,$separation)),0.025)} += $dtp;
+ $r->{logRout_over_separation}->{$population->rebin(safelog10($Rout/MAX(1e-50,$separation)),0.025)} += $dtp;
+ $r->{logHRin}->{$population->rebin(safelog10($HRin),0.025)} += $dtp;
+ $r->{logHRout}->{$population->rebin(safelog10($HRout),0.025)} += $dtp;
$r->{logLdisc}->{$population->rebin(safelog10($Ldisc),0.1)} += $dtp;
$r->{logLstar}->{$logLstar_binned} += $dtp;
$r->{logLXstar}->{$population->rebin(safelog10($LXstar),0.1)} += $dtp;
$r->{logdensityin}->{$population->rebin(safelog10($density_in),0.1)} += $dtp;
$r->{logdensityout}->{$population->rebin(safelog10($density_out),0.1)} += $dtp;
- $r->{logTin}->{$population->rebin(safelog10($Tin),0.1)} += $dtp;
- $r->{logTout}->{$population->rebin(safelog10($Tout),0.1)} += $dtp;
+ $r->{logTin}->{$population->rebin(safelog10($Tin),0.025)} += $dtp;
+ $r->{logTout}->{$population->rebin(safelog10($Tout),0.025)} += $dtp;
$r->{logPgrav_in}->{$population->rebin(safelog10($Pgrav_in),0.1)} += $dtp;
$r->{logPgrav_out}->{$population->rebin(safelog10($Pgrav_out),0.1)} += $dtp;
- $r->{logTeffmax}->{$population->rebin(safelog10($Teffmax),0.1)} += $dtp;
+ $r->{logTeffmax}->{$population->rebin(safelog10($Teffmax),0.025)} += $dtp;
$r->{logtvisc_in}->{$population->rebin(safelog10($tvisc_in),0.1)} += $dtp;
# 2D plots
@@ -504,16 +504,22 @@ sub parse_data
$systems{$x}++;
}
+
+ printf "DISC_END Set prog type %s loglifetime = %g from %g += %g\n",
+ $progenitor_stellar_type,
+ $population->rebin(safelog10($lifetime),0.1),
+ $lifetime,
+ $p;
foreach my $set ('all',
$progenitor_stellar_type)
{
my $r = $results->{$set};
- $r->{logendtime}->{$population->rebin(safelog10($endtime),0.1)} += $p;
- $r->{loglifetime}->{$population->rebin(safelog10($lifetime),0.1)} += $p;
+ $r->{logendtime}->{$population->rebin(safelog10($endtime),0.025)} += $p;
+ $r->{loglifetime}->{$population->rebin(safelog10($lifetime),0.025)} += $p;
$r->{final_eccentricity}->{$population->rebin($final_eccentricity,0.025)} += $p;
}
-
+
if(defined $initparamspace_fp)
{
@@ -550,6 +556,7 @@ sub parse_data
$init_a,
$init_e
) = @$la;
+ print "Start donor type $donor_stellar_type (== $progenitor_stellar_type ?)\n";
foreach my $set ('all',
$progenitor_stellar_type)
@@ -626,9 +633,14 @@ sub output
my $results = $population->results;
my $nhash = rob_misc::hashcount($population->results);
- if(($population->{_grid_options}{slurm_command} eq 'join' ||
- !defined $population->{_grid_options}{slurm_command})&&
- $nhash)
+ if(
+ (
+ $population->{_grid_options}{slurm} == 0 ||
+ $population->{_grid_options}{slurm_command} eq 'join' ||
+ !defined $population->{_grid_options}{slurm_command}
+ )
+ &&
+ $nhash)
{
printf "OUTPUT_ALL : hash size %d\n",$nhash;
diff --git a/src/stellar_structure/stellar_structure_HeStar.c b/src/stellar_structure/stellar_structure_HeStar.c
index 8f6e6c66eedab612e8cfc669e9430ce34ba49c26..630b2f4b838d3eb6b106afdcc31358f1767e188c 100644
--- a/src/stellar_structure/stellar_structure_HeStar.c
+++ b/src/stellar_structure/stellar_structure_HeStar.c
@@ -370,7 +370,7 @@ void stellar_structure_HeStar(struct star_t * newstar,
/*
* Mass is too low to ignite neon, hence
* the supernova is an electron capture,
- * not a canonical core collapse (which
+ * not a canonical core collapse which
* ignites neon, ..., silicon and forms
* an iron core.
*/
diff --git a/tbse b/tbse
index a1a4ab17c1a7f57b97cdf812db8c15de4b1e4db6..bd56573ae56bf79130bfbd2ac0d55a2c410de343 100755
--- a/tbse
+++ b/tbse
@@ -168,8 +168,12 @@ RANDOM_SKIP=0
# Core mass limits
MAX_TPAGB_CORE_MASS=1.38
CHANDRASEKHAR_MASS=1.44
+
+# mass above which neutron stars collapse to black holes
MAX_NEUTRON_STAR_MASS=1.8
+# stellar mass below which carbon fails to ignite
MINIMUM_MASS_FOR_CARBON_IGNITION=1.6
+# helium-star core mass lower limit for neon ignition
MINIMUM_MASS_FOR_NEON_IGNITION=2.85
# When the core of a red giant is exposed, its mass
@@ -575,7 +579,7 @@ DISC_TIMESTEP_FACTOR=1.0
# DISC_LOG_LEVEL_SUBTIMESTEP 2 = on with sub-timestep logging
# DISC_LOG_LEVEL_NORMAL_FIRST_DISC_ONLY = as 1 for first disc only
# DISC_LOG_LEVEL_SUBTIMESTEP_FIRST_DISC_ONLY = as 2 for first disc only
-DISC_LOG=DISC_LOG_LEVEL_NORMAL_FIRST_DISC_ONLY
+DISC_LOG=DISC_LOG_LEVEL_NONE
# soft 2D disc log switch, same options as DISC_LOG
DISC_LOG2D=DISC_LOG_LEVEL_NONE
@@ -1721,6 +1725,7 @@ then
shift
fi
+
# check for RANDOM_SEED
if [ -n "$RANDOM_SEED" ]
then
@@ -2232,7 +2237,6 @@ then
elif [ "$RUNMODE" == "gprof" ];
then
-
# run through gprof, suppressing stdout
# (requires CFLAGS="-pg" during configure step)
POSTCOMMANDS=("gprof $BIN_C")