From ae506200ae6e63a265fcf4450a023c468ebc323a Mon Sep 17 00:00:00 2001
From: Robert Izzard <r.izzard@surrey.ac.uk>
Date: Mon, 21 Jan 2019 20:56:47 +0000
Subject: [PATCH] fixed some issues with supernova:
post-SN orbital angular velocity was infinite (it was 1/0!) oops.
tides now act on everything except black holes and massless remnants.
also TWOPI*PI in wind-fed discs should have been TWOPI (this wasn't used, however)
---
src/binary_c_debug.h | 15 ++++----
src/binary_c_parameters.h | 2 +-
...ar_momentum_and_eccentricity_derivatives.c | 2 +-
src/binary_star_functions/tides.c | 19 +++++++---
src/binary_star_functions/tides.h | 22 +++++++++---
.../update_masses_and_angular_momenta.c | 2 +-
.../update_orbital_variables.c | 2 +-
src/disc/disc_stellar_wind_to_cbdisc.c | 2 +-
src/evolution/catch_events.c | 6 ++--
src/evolution/erase_events.c | 1 +
.../evolution_restore_from_previous.c | 1 +
src/evolution/evolution_split.c | 31 +++++++++++++---
src/stellar_structure/stellar_structure.c | 3 +-
src/stellar_structure/stellar_structure_BSE.c | 3 +-
src/stellar_structure/stellar_structure_NS.c | 9 ++---
src/supernovae/set_new_orbit.c | 17 +++++----
src/supernovae/supernova.c | 36 ++++++++++++-------
src/supernovae/supernova_kick.c | 2 ++
src/supernovae/supernovae.h | 2 +-
19 files changed, 124 insertions(+), 53 deletions(-)
diff --git a/src/binary_c_debug.h b/src/binary_c_debug.h
index ca1590134..a5ad663ff 100644
--- a/src/binary_c_debug.h
+++ b/src/binary_c_debug.h
@@ -20,7 +20,7 @@
*
* See below for further options.
*/
-#define DEBUG 1
+#define DEBUG 0
/*
* To remove debugging possibility from the code, globally, hence
@@ -56,8 +56,8 @@
* the same, but stardata is defined everywhere as a global
* variable.
*/
-//#define Debug_expression (1)
-#define Debug_expression (stardata!=NULL && stardata->model.time>3170 && stardata->model.time < 3190)
+#define Debug_expression (1)
+//#define Debug_expression (stardata!=NULL && stardata->model.time>3170 && stardata->model.time < 3190)
/*
* If you define Debug_stop_expression, and it is at any time TRUE,
@@ -72,8 +72,8 @@
*/
-//#define Debug_stop_expression (stardata->star[0].mass < stardata->star[0].core_mass)
-#undef Debug_stop_expression
+//#define Debug_stop_expression (isinf(stardata->common.orbital_angular_frequency))
+#undef Debug_stop_expression
/*
* Debug_show_expression, if defined, is shown at the beginning
@@ -85,13 +85,14 @@
* this is not shown.
*/
-#define Debug_show_expression " %g st=[%d,%d] omega=[%g,%g] Jorb=%g ", \
+#define Debug_show_expression " %g st=[%d,%d] omega=[%g,%g] Jorb=%g omega_orb=%g ", \
stardata->model.time, \
stardata->star[0].stellar_type, \
stardata->star[1].stellar_type, \
stardata->star[0].omega, \
stardata->star[1].omega, \
- stardata->common.orbital_angular_momentum
+ stardata->common.orbital_angular_momentum, \
+ stardata->common.orbital_angular_frequency
//#undef Debug_show_expression
diff --git a/src/binary_c_parameters.h b/src/binary_c_parameters.h
index 1b4333b8f..0fa04f0e0 100644
--- a/src/binary_c_parameters.h
+++ b/src/binary_c_parameters.h
@@ -1651,7 +1651,7 @@
//#define TPAGBTRACKS
/* log for Window To The Stars 2 */
-#define WTTS_LOG
+//#define WTTS_LOG
/* anti TZ project */
//#define ANTI_TZ
diff --git a/src/binary_star_functions/apply_orbital_angular_momentum_and_eccentricity_derivatives.c b/src/binary_star_functions/apply_orbital_angular_momentum_and_eccentricity_derivatives.c
index 0b2d5c17d..95e6c1ac0 100644
--- a/src/binary_star_functions/apply_orbital_angular_momentum_and_eccentricity_derivatives.c
+++ b/src/binary_star_functions/apply_orbital_angular_momentum_and_eccentricity_derivatives.c
@@ -47,7 +47,7 @@ Boolean apply_orbital_angular_momentum_and_eccentricity_derivatives(struct stard
/*
* Angular momentum
*/
- show_derivatives(stardata);
+ //show_derivatives(stardata);
common->orbital_angular_momentum +=
stardata->model.derivative[DERIVATIVE_ORBIT_ANGMOM] *
dt;
diff --git a/src/binary_star_functions/tides.c b/src/binary_star_functions/tides.c
index 7186a7d20..1394b54af 100644
--- a/src/binary_star_functions/tides.c
+++ b/src/binary_star_functions/tides.c
@@ -1,8 +1,6 @@
#include "../binary_c.h"
#include "tides.h"
-#undef DEBUG
-#define DEBUG 0
/*
* This function calculates the derivatives of
@@ -26,7 +24,7 @@ void tides(struct stardata_t *stardata,
double sq6 = POW6(sqome2);
double isqome2_13 = 1.0/(sq6*sq6*sqome2);
double starI = moment_of_inertia(star,radius);
- const int ndonor = stardata->model.ndonor;
+ const int ndonor MAYBE_UNUSED = stardata->model.ndonor;
/*
* Calculate tidal effects : circularization, orbital changes and spin up.
@@ -163,8 +161,19 @@ void tides(struct stardata_t *stardata,
(3.0*cq*tcqr/(rg2*POW6(omecc2)))*
(f2*stardata->common.orbital_angular_frequency - sqome3f5*star->omega);
- Dprint("Spin up star->dspint = %g\n",star->derivative[DERIVATIVE_STELLAR_ANGULAR_VELOCITY]);
-
+ Dprint("Spin up star->dspint = %g from cq=%g tcqr=%g rg2=%g omecc2=%g f2=%g sqome3f5=%g omega_orb=%g omega_star=%g\n",
+ star->derivative[DERIVATIVE_STELLAR_ANGULAR_VELOCITY],
+ cq,
+ tcqr,
+ rg2,
+ omecc2,
+ f2,
+ sqome3f5,
+ stardata->common.orbital_angular_frequency,
+ star->omega
+ );
+ NANCHECK(star->derivative[DERIVATIVE_STELLAR_ANGULAR_VELOCITY]);
+
/*
* Calculate the equilibrium spin at which no angular momentum
* can be transferred. The stellar angular velocity derivative
diff --git a/src/binary_star_functions/tides.h b/src/binary_star_functions/tides.h
index bfa922418..9be2a9bb3 100644
--- a/src/binary_star_functions/tides.h
+++ b/src/binary_star_functions/tides.h
@@ -1,10 +1,24 @@
#ifndef TIDES_H
#define TIDES_H
-#define DO_TIDES (stardata->preferences->tidal_strength_factor>TINY && \
- ((star->stellar_type<=HeGB && \
- MORE_OR_EQUAL(star->radius,0.01*star->roche_radius))|| \
- (star->stellar_type>=HeWD && k==ndonor)))
+/* BSE definition of when to allow tides */
+#define DO_TIDES ( \
+ /* tides must be on */ \
+ stardata->preferences->tidal_strength_factor>TINY \
+ && \
+ ( \
+ /* star is pre-WD */ \
+ ( \
+ star->stellar_type<HeWD && \
+ MORE_OR_EQUAL(star->radius,1e-2*star->roche_radius) \
+ ) \
+ || \
+ /* or star is a WD donor */ \
+ (star->stellar_type>=HeWD && k==ndonor)))
+
+/* better: always do tides except for black holes and massless remnants */
+#undef DO_TIDES
+#define DO_TIDES (star->stellar_type < BLACK_HOLE)
#ifdef PRE_MAIN_SEQUENCE
diff --git a/src/binary_star_functions/update_masses_and_angular_momenta.c b/src/binary_star_functions/update_masses_and_angular_momenta.c
index b0ca058d8..689d517b8 100644
--- a/src/binary_star_functions/update_masses_and_angular_momenta.c
+++ b/src/binary_star_functions/update_masses_and_angular_momenta.c
@@ -114,7 +114,7 @@ int update_masses_and_angular_momenta(struct stardata_t * stardata,
NANCHECK(stardata->common.orbital_angular_momentum);
Boolean retval = apply_orbital_angular_momentum_and_eccentricity_derivatives(stardata,RLOF);
- Dprint("post-apply orb and ecc : J=%g nan %d\n",stardata->common.orbital_angular_momentum,isnan(stardata->common.orbital_angular_momentum));
+ Dprint("post-apply orb and ecc : J=%g NaN? %d\n",stardata->common.orbital_angular_momentum,isnan(stardata->common.orbital_angular_momentum));
NANCHECK(stardata->common.orbital_angular_momentum);
return(retval);
diff --git a/src/binary_star_functions/update_orbital_variables.c b/src/binary_star_functions/update_orbital_variables.c
index 68291a27b..a92dde6d6 100644
--- a/src/binary_star_functions/update_orbital_variables.c
+++ b/src/binary_star_functions/update_orbital_variables.c
@@ -25,7 +25,7 @@ void update_orbital_variables(struct stardata_t * RESTRICT stardata)
else
{
double denom = POW2(MASS_PRODUCT*TWOPI)*(POW3(AU_IN_SOLAR_RADII)*(1.0-ECC_SQUARED));
- stardata->common.separation = denom > TINY ? total_mass*JORB_SQUARED/denom : 0.0;
+ stardata->common.separation = denom > TINY ? (total_mass*JORB_SQUARED/denom) : 0.0;
}
stardata->model.derivative[DERIVATIVE_ORBIT_SEMI_MAJOR_AXIS] =
diff --git a/src/disc/disc_stellar_wind_to_cbdisc.c b/src/disc/disc_stellar_wind_to_cbdisc.c
index b80b52340..48d2df4c8 100644
--- a/src/disc/disc_stellar_wind_to_cbdisc.c
+++ b/src/disc/disc_stellar_wind_to_cbdisc.c
@@ -102,7 +102,7 @@ void disc_stellar_wind_to_cbdisc(struct stardata_t * stardata)
/*
* Angular frequency in radians per year
*/
- double orbital_angular_frequency = TWOPI*PI/orbital_period_years;
+ double orbital_angular_frequency = TWOPI/orbital_period_years;
/*
* And in radians per second
diff --git a/src/evolution/catch_events.c b/src/evolution/catch_events.c
index 5c0b1edd7..a2c7c0535 100644
--- a/src/evolution/catch_events.c
+++ b/src/evolution/catch_events.c
@@ -38,6 +38,7 @@ void catch_events(struct stardata_t * stardata)
* kicks, mass change, yields, etc.
*/
supernova(stardata,star);
+ Dprint("after supernova\n");
/*
* Check for broken system
@@ -60,10 +61,11 @@ void catch_events(struct stardata_t * stardata)
IS_ZERO(stardata->common.separation) ? 0.0 :
((stardata->common.separation/AU_IN_SOLAR_RADII)
*sqrt(stardata->common.separation/(AU_IN_SOLAR_RADII*(stardata->star[0].mass + stardata->star[1].mass))));
-
+
+ /* and the orbital angular frequency */
stardata->common.orbital_angular_frequency =
IS_ZERO(stardata->common.separation) ? 0.0 :
- (TWOPI/stardata->common.orbital_period)/0.0;
+ (TWOPI/stardata->common.orbital_period);
}
}
}
diff --git a/src/evolution/erase_events.c b/src/evolution/erase_events.c
index def0a98f5..716947bdb 100644
--- a/src/evolution/erase_events.c
+++ b/src/evolution/erase_events.c
@@ -12,6 +12,7 @@ void erase_events(struct stardata_t * stardata)
STARLOOP(k)
{
SETstar(k);
+ Dprint("call free supernova\n");
free_supernova(stardata,star);
}
}
diff --git a/src/evolution/evolution_restore_from_previous.c b/src/evolution/evolution_restore_from_previous.c
index 038dab6da..60e88715d 100644
--- a/src/evolution/evolution_restore_from_previous.c
+++ b/src/evolution/evolution_restore_from_previous.c
@@ -1,6 +1,7 @@
#include "../binary_c.h"
+
void evolution_restore_from_previous(struct stardata_t * RESTRICT stardata)
{
/*
diff --git a/src/evolution/evolution_split.c b/src/evolution/evolution_split.c
index 43aecd73d..327b526d8 100644
--- a/src/evolution/evolution_split.c
+++ b/src/evolution/evolution_split.c
@@ -67,6 +67,18 @@ int testindex = NUCSYN_ENSEMBLE_RSG_LUM;
#endif // EVOLUTION_SPLITTING
+#define DDD \
+ SDEBUG("stardata new_supernovae %p %p st %d %d omega %g %g J %g %g Jorb %g\n", \
+ p->stardata->star[0].new_supernova, \
+ p->stardata->star[1].new_supernova, \
+ p->stardata->star[0].stellar_type, \
+ p->stardata->star[1].stellar_type, \
+ p->stardata->star[0].omega, \
+ p->stardata->star[1].omega, \
+ p->stardata->star[0].angular_momentum, \
+ p->stardata->star[1].angular_momentum, \
+ p->stardata->common.orbital_angular_momentum \
+ );
int evolution_split(struct stardata_t * RESTRICT stardata,
@@ -100,7 +112,10 @@ int evolution_split(struct stardata_t * RESTRICT stardata,
* Check if we need to rerun from a previous
* splitpoint
*/
- SDEBUG("SPLIT status is evolution_break : currently working on splitinfo struct %d at %p, of which count = %d\n",stardata->preferences->current_splitdepth-1,p,p->count);
+ SDEBUG("SPLIT status is evolution_break : currently working on splitinfo struct %d at %p, of which count = %d\n",
+ stardata->preferences->current_splitdepth-1,
+ p,
+ p->count);
if(p->count > 0)
{
@@ -116,6 +131,8 @@ int evolution_split(struct stardata_t * RESTRICT stardata,
Random_buffer random_buffer = stardata->common.random_buffer;
Boolean random_buffer_set = stardata->common.random_buffer_set;
+ DDD;
+
p->count--;
Copy_stardata(p->stardata,stardata);
@@ -128,7 +145,10 @@ int evolution_split(struct stardata_t * RESTRICT stardata,
stardata->preferences->current_splitdepth - 1,
p->count,
stardata->model.time);
- Append_logstring(LOG_SPLIT_EVOLUTION,"Evolution restart @ t=%g, %d.%d",stardata->model.time,stardata->preferences->current_splitdepth-1,p->count);
+ Append_logstring(LOG_SPLIT_EVOLUTION,"Evolution restart @ t=%g, %d.%d",
+ stardata->model.time,
+ stardata->preferences->current_splitdepth-1,
+ p->count);
#ifdef LINUX
reset_binary_c_timeout();
#endif
@@ -137,7 +157,9 @@ int evolution_split(struct stardata_t * RESTRICT stardata,
else
{
stardata->preferences->current_splitdepth--;
- SDEBUG("SPLIT last done (count=%d) : next depth is %d\n",p->count,stardata->preferences->current_splitdepth-1);
+ SDEBUG("SPLIT last done (count=%d) : next depth is %d\n",
+ p->count,
+ stardata->preferences->current_splitdepth-1);
/* clean up memory */
Safe_free(stardata->preferences->splitinfo[stardata->preferences->current_splitdepth]->stardata);
@@ -158,7 +180,7 @@ int evolution_split(struct stardata_t * RESTRICT stardata,
SDEBUG("SPLIT AT %g : split depth = %d\n",
stardata->model.time,
stardata->preferences->current_splitdepth);
-
+
/* add new struct */
stardata->preferences->splitinfo[stardata->preferences->current_splitdepth] = malloc(sizeof(struct splitinfo_t));
@@ -174,6 +196,7 @@ int evolution_split(struct stardata_t * RESTRICT stardata,
/* allocate space and copy previous stardata for restore */
p->stardata = New_stardata_from(stardata->previous_stardata);
SDEBUG("starting post-split at t=%g\n",p->stardata->model.time);
+ DDD;
/*
* modulate probability, yields, ensemble etc. which
diff --git a/src/stellar_structure/stellar_structure.c b/src/stellar_structure/stellar_structure.c
index 32c115e04..16f84c4c5 100644
--- a/src/stellar_structure/stellar_structure.c
+++ b/src/stellar_structure/stellar_structure.c
@@ -80,7 +80,7 @@ void stellar_structure(struct star_t *star,
* We have a SN, but are not allowed to have one.
* Set the stellar structure from the new supernova
* struct and free the associated memory.
- */
+ */
if(newstar->new_supernova != NULL)
{
Dprint("Remove new supernova : SNe not allowed\n");
@@ -121,6 +121,7 @@ void stellar_structure(struct star_t *star,
Safe_free(newstar);
+ Dprint("return \n");
return;
#ifdef OLDCODE
diff --git a/src/stellar_structure/stellar_structure_BSE.c b/src/stellar_structure/stellar_structure_BSE.c
index f96f4456b..b0606b5d0 100644
--- a/src/stellar_structure/stellar_structure_BSE.c
+++ b/src/stellar_structure/stellar_structure_BSE.c
@@ -692,7 +692,8 @@ static Boolean helium_ignition_test(struct star_t * newstar,
ignite_on_HG
);
}
-
+
+ Dprint("return %d\n",helium_ignition);
return helium_ignition;
}
#endif
diff --git a/src/stellar_structure/stellar_structure_NS.c b/src/stellar_structure/stellar_structure_NS.c
index 1a9123845..07b047181 100644
--- a/src/stellar_structure/stellar_structure_NS.c
+++ b/src/stellar_structure/stellar_structure_NS.c
@@ -8,7 +8,7 @@ void stellar_structure_NS(struct star_t * newstar,
newstar->core_mass = newstar->mass;
newstar->GB_core_mass = 0.0;
newstar->CO_core_mass = newstar->mass;
- newstar->luminosity=lns(newstar->mass,newstar->age);
+ newstar->luminosity = lns(newstar->mass,newstar->age);
newstar->radius = rns(newstar->mass);
if(newstar->core_mass > stardata->preferences->max_neutron_star_mass)
@@ -21,13 +21,14 @@ void stellar_structure_NS(struct star_t * newstar,
if(newstar->SN_type == SN_NONE)
{
newstar->SN_type = SN_AIC_BH;
- struct star_t * newsSN =
+ struct star_t * news =
new_supernova(stardata,
newstar,
Other_star_struct(newstar),
newstar);
- newsSN->stellar_type = BLACK_HOLE;
- newsSN->age = 0.0;
+ news->stellar_type = BLACK_HOLE;
+ news->age = 0.0;
+ news->mass = newstar->mass;
}
}
diff --git a/src/supernovae/set_new_orbit.c b/src/supernovae/set_new_orbit.c
index 202e167a2..f93f71a5e 100644
--- a/src/supernovae/set_new_orbit.c
+++ b/src/supernovae/set_new_orbit.c
@@ -97,7 +97,7 @@ void set_new_orbit(struct stardata_t * stardata,
kick_dispersion,
kick_companion);
}
-
+
Dprint("kick_speed = %g, orbital_speed = %g, new_orbital_speed_squared = %g\n",
kick_system->kick_speed,
kick_system->orbital_speed,
@@ -123,14 +123,18 @@ void set_new_orbit(struct stardata_t * stardata,
*separation);
}
- SNprint("SN new separation %g from r=%g vn2=%g new_system_mass=%g\n",
- *separation,
- kick_system->separation,
- kick_system->new_orbital_speed_squared,
- new_system_mass
+ SNprint("SN new separation %g from r=%g vn2=%g new_system_mass=%g : old omega_orb=%g %g\n",
+ *separation,
+ kick_system->separation,
+ kick_system->new_orbital_speed_squared,
+ new_system_mass,
+ pre_explosion_stardata->common.orbital_angular_frequency,
+ stardata->common.orbital_angular_frequency
);
+
Dprint("Pre-sep calc : sgl = %d\n",stardata->model.sgl);
+
if(*separation < MINIMUM_SEPARATION_TO_BE_CALLED_BINARY ||
*separation > MAXIMUM_SEPARATION_TO_BE_CALLED_BINARY ||
stardata->model.sgl)
@@ -204,5 +208,6 @@ void set_new_orbit(struct stardata_t * stardata,
Dprint("New COM velocity (%g, %g, %g)\n",vsx,vsy,vsz);
}
+
}
diff --git a/src/supernovae/supernova.c b/src/supernovae/supernova.c
index a065c627f..0523f3eda 100644
--- a/src/supernovae/supernova.c
+++ b/src/supernovae/supernova.c
@@ -186,26 +186,35 @@ void supernova(struct stardata_t * stardata,
*/
star->omega = pre_explosion_star->omega;
//star->omega=2.0e8;
- printf("SNomega %g from pre explosion M=%g Mc=%g R=%g kw=%d omega=%g\n",
+
+ const double I = K3*POW2(star->radius)*star->mass;
+
+ star->angular_momentum = I * star->omega;
+
+ const double Jmax = GRAVITATIONAL_CONSTANT * POW2(star->mass * M_SUN) /
+ SPEED_OF_LIGHT / ANGULAR_MOMENTUM_CGS;
+
+ /*
+ printf("SNomega %g from pre explosion M=%g Mc=%g R=%g kw=%d omega=%g : remNaNt J = %g, Jmax = %g\n",
star->omega,
pre_explosion_star->mass,
pre_explosion_star->core_mass,
pre_explosion_star->radius,
pre_explosion_star->stellar_type,
- pre_explosion_star->omega
+ pre_explosion_star->omega,
+ star->angular_momentum,
+ Jmax
);
fflush(NULL);
-
- star->angular_momentum = K3*POW2(star->radius)*
- star->mass*star->omega;
-
+ */
+
/*
* Object cannot rotate beyond GM^2/c
*/
star->angular_momentum =
- MIN(GRAVITATIONAL_CONSTANT * POW2(star->mass * M_SUN) /
- SPEED_OF_LIGHT / ANGULAR_MOMENTUM_CGS,
- star->angular_momentum);
+ MIN(Jmax, star->angular_momentum);
+
+ star->omega = star->angular_momentum / I;
}
/*
@@ -245,13 +254,12 @@ void supernova(struct stardata_t * stardata,
0.0,
star->new_supernova->new_companion_structure->stellar_type
);
+#endif // NUCSYN && NUCSYN_SUPERNOVAE
-
-#endif
-
/*
* Logging
*/
+ Dprint("Logging\n");
#ifdef LOG_SUPERNOVAE
log_sn(stardata,star->starnum,star->SN_type,POST_SN);
#endif
@@ -271,7 +279,7 @@ void supernova(struct stardata_t * stardata,
#endif
#if defined DISCS && defined DISC_LOG_POPSYN
-
+ Dprint("Logging (popsyn)\n");
int i;
for(i=0; i<stardata->common.ndiscs; i++)
{
@@ -303,7 +311,9 @@ void supernova(struct stardata_t * stardata,
/*
* Free the memory associated with the new_supernova
*/
+ Dprint("call free supernova\n");
free_supernova(stardata,star);
Safe_free(pre_explosion_star);
Safe_free(pre_explosion_stardata);
+ Dprint("Done\n");
}
diff --git a/src/supernovae/supernova_kick.c b/src/supernovae/supernova_kick.c
index d3d71c5f3..018ef4869 100644
--- a/src/supernovae/supernova_kick.c
+++ b/src/supernovae/supernova_kick.c
@@ -26,4 +26,6 @@ void supernova_kick(struct stardata_t * stardata,
stardata,
pre_explosion_stardata);
}
+
+
}
diff --git a/src/supernovae/supernovae.h b/src/supernovae/supernovae.h
index d4bbca911..67eb208a1 100644
--- a/src/supernovae/supernovae.h
+++ b/src/supernovae/supernovae.h
@@ -7,7 +7,7 @@
FEQUAL((kick_companion),-(A)))
#if defined DEBUG && DEBUG==1
-#define SNprint(...) printf(__VA_ARGS__)
+#define SNprint(...) printf(__VA_ARGS__);fflush(NULL);
#else
#define SNprint(...)
#endif
--
GitLab