diff --git a/src/binary_c_parameters.h b/src/binary_c_parameters.h
index 253a2b10932f68e81c648a653ebd7e7b6dcff9cb..ca7d451f6be052600cb1f502eebd5d171d10097c 100644
--- a/src/binary_c_parameters.h
+++ b/src/binary_c_parameters.h
@@ -718,7 +718,7 @@
* NOVAE_TIMESTEP_FACTOR of the recurrance time
* (time between nova explosions).
*/
-//#define TIME_RESOLVE_NOVAE
+#define TIME_RESOLVE_NOVAE
//#define NOVAE_TIMESTEP_FACTOR 0.1
/*
diff --git a/src/binary_c_structures.h b/src/binary_c_structures.h
index 93289a383f8e7a6345016bf53a31131e34dcefa6..cf032a6f19b17dd4c923aa0ec2c5072f83c805a3 100644
--- a/src/binary_c_structures.h
+++ b/src/binary_c_structures.h
@@ -1368,7 +1368,7 @@ struct star_t {
#endif// TIME_RESOLVE_NOVAE
#ifdef INDIVIDUAL_NOVAE
- double dm_novaH;
+ double dm_novaH,nova_timeout;
#endif //INDIVIDUAL_NOVAE
double nova_beta,nova_faml;
@@ -1554,7 +1554,7 @@ struct star_t {
#ifdef LOG_SUPERNOVAE
Boolean supernova;
#endif//LOG_SUPERNOAVE
-
+ Boolean nova;
#ifdef NUCSYN
Boolean rstar; /* true if star is R type */
Boolean first_dredge_up; /* stores whether we've had first dredge up */
diff --git a/src/binary_star_functions/angular_momentum_and_eccentricity_derivatives.c b/src/binary_star_functions/angular_momentum_and_eccentricity_derivatives.c
index 075770d689fab8351b506cb9a39bf9a3ba04900f..908fa54af12d70893313842490ac0837710e488f 100644
--- a/src/binary_star_functions/angular_momentum_and_eccentricity_derivatives.c
+++ b/src/binary_star_functions/angular_momentum_and_eccentricity_derivatives.c
@@ -201,32 +201,11 @@ void angular_momentum_and_eccentricity_derivatives(
SETstars(k);
if(Is_not_zero(star->derivative[DERIVATIVE_STELLAR_MASS_NOVA]))
{
- double q = companion->mass / star->mass;
-
- /*
- * Assume novae lose mass spherically,
- * and the companion gains mass in a Keplerian disc.
- */
- double k = star->derivative[DERIVATIVE_STELLAR_MASS_NOVA] < 0.0
- ? (2.0/3.0) : 1.0;
-
- star->derivative[DERIVATIVE_STELLAR_ANGMOM_NOVA] =
- star->derivative[DERIVATIVE_STELLAR_MASS_NOVA] *
- k * Pow2(star->radius) * star->omega;
-
- stardata->model.derivative[DERIVATIVE_ORBIT_ANGMOM_NOVA] +=
- q / (1.0 + q) *
- star->derivative[DERIVATIVE_STELLAR_MASS_NOVA] / star->mass
- * stardata->common.orbit.angular_momentum;
-
- if(star->derivative[DERIVATIVE_STELLAR_MASS_NOVA] < 0.0)
- {
- /* this is the noving WD, only let this contribute to FAML */
- stardata->model.derivative[DERIVATIVE_ORBIT_ANGMOM_NOVA] +=
- star->nova_faml *
- star->derivative[DERIVATIVE_STELLAR_MASS_NOVA] *
- stardata->common.orbit.angular_momentum;
- }
+ nova_angular_momentum_changes(stardata,
+ star,
+ star->derivative[DERIVATIVE_STELLAR_MASS_NOVA],
+ &star->derivative[DERIVATIVE_STELLAR_ANGMOM_NOVA],
+ &stardata->model.derivative[DERIVATIVE_ORBIT_ANGMOM_NOVA]);
}
}
diff --git a/src/binary_star_functions/limit_accretion_rates.c b/src/binary_star_functions/limit_accretion_rates.c
index 322c2311f3c9a1fe0dfcd7bd56de921b9f2a17d2..a7c1261cd79a4da66e3f06e0cf19b72e65e9fc1e 100644
--- a/src/binary_star_functions/limit_accretion_rates.c
+++ b/src/binary_star_functions/limit_accretion_rates.c
@@ -20,6 +20,7 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
{
struct star_t * accretor = &stardata->star[k];
struct star_t * donor = &stardata->star[Other_star(k)];
+ accretor->nova = FALSE;
Dprint("limit acc rate star %d : is currently %g\n",
k,
@@ -50,7 +51,7 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
mdot *= YEAR_LENGTH_IN_SECONDS / M_SUN;
accretor->derivative[DERIVATIVE_STELLAR_MASS_RLOF_GAIN] += mdot;
*/
-
+
/*
* Limit mass accretion rate.
* Only applies if the accretion rate > 0
@@ -167,7 +168,7 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
* (e.g. novae)
*/
double accretion_rate = Mdot_net(accretor);
-
+
if(accretion_rate > 0.0 &&
COMPACT_OBJECT(accretor->stellar_type) &&
!COMPACT_OBJECT(donor->stellar_type))
@@ -239,6 +240,7 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
{
steady_burn_rate = steady_burn_rate_He;
nova_type = NOVA_TYPE_HELIUM;
+ accretor->nova = TRUE;
}
/* restore stellar type */
@@ -280,6 +282,7 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
* Angular momentum is dealth with in
* angular_momentum_and_eccentricity_derivatives()
*/
+
const double f = nova_retention_fraction(stardata,
accretor,
donor,
@@ -287,6 +290,7 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
steady_burn_rate);
Dprint("nova_retention_factor f = %g\n",f);
+
/*
* Accretion-induced irradiation
@@ -364,13 +368,7 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
*/
accretor->derivative[DERIVATIVE_STELLAR_H_LAYER_MASS] = accretion_rate;
- printf("dm NOVA %g\n",
- accretor->dm_novaH);
-
#ifdef INDIVIDUAL_NOVAE_EVENTS
- printf("cf DM %g vs dMCRIT %g\n",
- accretor->dm_novaH,
- dMcrit);
/*
* Check if we have sufficient mass for explosion
*/
@@ -386,7 +384,7 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
new_nova->accretion_rate = accretion_rate;
new_nova->steady_burn_rate = steady_burn_rate;
new_nova->f = f;
- printf("novae event\n");
+ accretor->nova = TRUE;
if(Add_new_event(stardata,
BINARY_C_EVENT_NOVA,
&nova_event_handler,
@@ -423,17 +421,19 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
* in the layer. We burn it all to maintain mass
* conservation.
*/
+ accretor->nova = TRUE;
if(accretor->novastate == NOVA_STATE_TRIGGER)
{
/*
* A nova has been triggered.
*
- * The timestep is (at most) 1 year.
+ * The timestep is short.
*
* Calculate dm_lost, the mass lost from
* the white dwarf in the nova explosion,
* hence the nova mass loss rate.
*/
+
const double dm_lost = (1.0-f) * accretor->dm_novaH;
const double dm_loss_rate = - dm_lost / stardata->model.dt;
accretor->derivative[DERIVATIVE_STELLAR_MASS_NOVA] +=
@@ -487,6 +487,7 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
if(accretor->novastate != NOVA_STATE_NONE)
{
+ accretor->nova = TRUE;
accretor->novastate = NOVA_STATE_QUIESCENT;
}
}
@@ -501,6 +502,7 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
*/
accretor->derivative[DERIVATIVE_STELLAR_MASS_NOVA] +=
- (1.0-f) * accretion_rate;
+ accretor->nova = TRUE;
Dprint("accretion rate %g -> dM/dt nova = %g\n",
accretion_rate,
accretor->derivative[DERIVATIVE_STELLAR_MASS_NOVA]);
@@ -590,6 +592,13 @@ void limit_accretion_rates(struct stardata_t * RESTRICT const stardata)
accretor->novastate =
accretor->novastate == NOVA_STATE_NONE ?
NOVA_STATE_NONE : NOVA_STATE_QUIESCENT;
+
+ if((accretor->novastate == NOVA_STATE_QUIESCENT ||
+ Is_not_zero(accretor->derivative[DERIVATIVE_STELLAR_MASS_NOVA])) &&
+ (Is_not_zero(accretor->nova_timeout) && stardata->model.time > accretor->nova_timeout))
+ {
+ accretor->nova = TRUE;
+ }
#endif
}
}
diff --git a/src/evolution/evolution_difflog.c b/src/evolution/evolution_difflog.c
index 9274b67818e31cded2ebcb75f124d403178b2cd4..a62434d91d2f17b3cd15206637c71d4083cac25a 100644
--- a/src/evolution/evolution_difflog.c
+++ b/src/evolution/evolution_difflog.c
@@ -386,7 +386,7 @@ static void set_diffstats(struct stardata_t * stardata,
now->luminosity[k] = star->luminosity;
now->novarate[k] = star->derivative[DERIVATIVE_STELLAR_MASS_NOVA];
#ifdef INDIVIDUAL_NOVAE
- now->novastate[k] = star->novastate;
+ now->novastate[k] = star->nova;
#endif
now->v_eq_ratio[k] = star->v_eq_ratio;
now->SN_type[k] = star->SN_type;
@@ -480,13 +480,13 @@ static void cf_diffstats(struct stardata_t * RESTRICT const stardata,
#else
Boolean novae = Boolean_(now->novarate[k] < 0.0);
#endif
- if(novae == TRUE &&
+ if((novae == TRUE) &&
lockflags[k][LOCK_NOVAE] == FALSE)
{
lockflags[k][LOCK_NOVAE] = TRUE;
flags[BEG_NOVAE_LABEL] = TRUE;
}
- else if(novae == FALSE &&
+ else if((novae == FALSE) &&
lockflags[k][LOCK_NOVAE] == TRUE)
{
lockflags[k][LOCK_NOVAE] = FALSE;
diff --git a/src/logging/log_every_timestep.c b/src/logging/log_every_timestep.c
index 3a34272029604375de69b80f679f2af0f0deab45..e783f1f6dc1f1d5ecff83a1a2f312f5b1f545976 100644
--- a/src/logging/log_every_timestep.c
+++ b/src/logging/log_every_timestep.c
@@ -2191,7 +2191,7 @@ void log_every_timestep(struct stardata_t * RESTRICT const stardata)
#ifdef INDIVIDUAL_NOVAE
if(WHITE_DWARF(stardata->star[0].stellar_type))
{
- printf("NOVAE %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e\n",
+ printf("NOVAE %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %d %20.12e %20.12e %20.12e\n",
stardata->model.time,
stardata->star[0].mass,
dt,
@@ -2200,7 +2200,11 @@ void log_every_timestep(struct stardata_t * RESTRICT const stardata)
stardata->star[1].radius,
Mdot_net(&stardata->star[1]),
stardata->star[0].dm_novaH, //8
- stardata->star[0].derivative[DERIVATIVE_STELLAR_MASS_NOVA]
+ stardata->star[0].derivative[DERIVATIVE_STELLAR_MASS_NOVA],
+ (stardata->star[0].nova || stardata->star[1].nova) ? 1 : 0, //10
+ Mdot_net(&stardata->star[0]), //11 accretor rate (+ve when accreting, -ve when exploding)
+ Mdot_net(&stardata->star[1]), //12 donor rate (-ve)
+ stardata->common.orbit.separation //13
);
}
#endif//INDIVIDUAL_NOVAE
diff --git a/src/novae/nova_angular_momentum_changes.c b/src/novae/nova_angular_momentum_changes.c
new file mode 100644
index 0000000000000000000000000000000000000000..c669ff91ef1dc648fe83ec9b14f553bf839d312c
--- /dev/null
+++ b/src/novae/nova_angular_momentum_changes.c
@@ -0,0 +1,51 @@
+#include "../binary_c.h"
+
+void nova_angular_momentum_changes(
+ struct stardata_t * const stardata,
+ struct star_t * const star,
+ const double dm,
+ double * const dJstar,
+ double * const dJorb
+ )
+{
+ /*
+ * Angular momentum changes that happen because of novae.
+ *
+ * dm is either a mass lost (e.g. in a nova event) or the
+ * time deriative. If it is a derivative, then similarly
+ * dJstar and dJorb, the changes of the stellar and orbital
+ * angular momenta, should point to the appropriate derivatives.
+ *
+ * Usually dm is negative, because the exploding star loses mass.
+ *
+ * If dm is positive, this function should deal with nova re-accretion.
+ */
+
+ struct star_t * const companion =
+ Other_star_struct(star);
+
+ const double q = companion->mass / star->mass;
+
+ /*
+ * Assume novae lose mass spherically,
+ * and the companion gains mass in a Keplerian disc.
+ */
+ const double k = dm < 0.0 ? (2.0/3.0) : 1.0;
+
+ *dJstar =
+ dm * k * Pow2(star->radius) * star->omega;
+
+ *dJorb =
+ dm / star->mass *
+ q / (1.0 + q) *
+ stardata->common.orbit.angular_momentum;
+
+ if(dm < 0.0)
+ {
+ /* this is the noving WD, only let this contribute to FAML */
+ *dJorb +=
+ star->nova_faml *
+ dm *
+ stardata->common.orbit.angular_momentum;
+ }
+}
diff --git a/src/novae/nova_event_handler.c b/src/novae/nova_event_handler.c
index e1e1c8a9a4d23bc82c1c765c0444105a6e7e28f4..2b266494b64afac1ac932f5e14f9d06392c96d24 100644
--- a/src/novae/nova_event_handler.c
+++ b/src/novae/nova_event_handler.c
@@ -16,50 +16,100 @@ Event_handler_function nova_event_handler(void * const eventp,
(struct binary_c_nova_event_t *) data;
struct star_t * const donor = nova_data->donor;
- struct star_t * const accretor = nova_data->accretor;
- const double f = nova_data->f;
+ struct star_t * const exploder = nova_data->accretor;
- printf("caught nova event at %20.12g : accretor = %d with H-layer mass %g, donor = %d : f = %g\n",
+ printf("caught nova event at %20.12g : exploder = %d with H-layer mass %g, donor = %d : f = %g\n",
stardata->model.time,
- accretor->starnum,
- accretor->dm_novaH,
+ exploder->starnum,
+ exploder->dm_novaH,
donor->starnum,
- f
+ nova_data->f
);
-
- const double dm_lost = (1.0-f) * accretor->dm_novaH;
+
+ /*
+ * Calculate mass lost: note that dm_lost is positive,
+ * while dm_loss_rate is negative.
+ */
+ const double dm_lost = (1.0 - nova_data->f) * exploder->dm_novaH;
const double dm_loss_rate = - dm_lost / stardata->model.dt;
- accretor->derivative[DERIVATIVE_STELLAR_MASS_NOVA] +=
- dm_loss_rate;
+ /*
+ * Simulate a nova in the derivative
+ * (although this is only really useful for logging)
+ */
+ exploder->derivative[DERIVATIVE_STELLAR_MASS_NOVA] = dm_loss_rate;
+
/*
* No hydrogen is left on the surface of the white dwarf
*/
- accretor->dm_novaH = 0.0;
- accretor->mass -= dm_lost;
- accretor->novastate = NOVA_STATE_POST;
-
+ exploder->dm_novaH = 0.0;
+ exploder->mass -= dm_lost;
+ exploder->novastate = NOVA_STATE_POST;
+
/*
* yield nova material
*/
double Xnova[ISOTOPE_ARRAY_SIZE];
nucsyn_set_nova_abunds(stardata,
- accretor,
+ exploder,
donor->Xenv,//THIS IS WRONG
Xnova);
nucsyn_calc_yields(stardata,
- accretor,
+ exploder,
dm_lost,
Xnova,
0.0,
NULL,
- accretor->starnum,
+ exploder->starnum,
YIELD_NOT_FINAL,
NUCSYN_SOURCE_NOVAE);
/*
- * Todo : re-accretion of a fraction
+ * Re-accretion of a fraction
+ */
+ const double dm_gain = exploder->nova_beta * dm_lost;
+ donor->mass += dm_gain;
+
+ /*
+ * update stellar and orbital angular momenta
+ * caused by mass loss from the exploder
+ */
+ double dJstar,dJorb;
+ nova_angular_momentum_changes(
+ stardata,
+ exploder,
+ -dm_lost,
+ &dJstar,
+ &dJorb
+ );
+ exploder->angular_momentum += dJstar;
+ stardata->common.orbit.angular_momentum += dJorb;
+ fprintf(stderr,"NOVA exploder %g %g\n",dJstar,dJorb);
+
+ /*
+ * update stellar and orbital angular momenta
+ * caused by mass gained onto the donor
+ */
+ nova_angular_momentum_changes(
+ stardata,
+ donor,
+ dm_gain,
+ &dJstar,
+ &dJorb
+ );
+ donor->angular_momentum += dJstar;
+ stardata->common.orbit.angular_momentum += dJorb;
+ fprintf(stderr,"NOVA gainer (reaccretion) %g %g\n",dJstar,dJorb);
+
+ /*
+ * update orbit
*/
+ fprintf(stderr,"sep was %20.12e .. ",stardata->common.orbit.separation);
+ update_orbital_variables(stardata,
+ &stardata->common.orbit,
+ &stardata->star[0],
+ &stardata->star[1]);
+ fprintf(stderr,"now %20.12e\n",stardata->common.orbit.separation);
return NULL;
}
diff --git a/src/novae/novae_prototypes.h b/src/novae/novae_prototypes.h
index 6c4142de1fc0bf458fdec5d79feff5d949b63482..d4f14ac3263794bd26d7eaac2c79da89f4418878 100644
--- a/src/novae/novae_prototypes.h
+++ b/src/novae/novae_prototypes.h
@@ -29,5 +29,13 @@ Event_handler_function nova_erase_event_handler(void * eventp,
void * data);
#endif // INDIVIDUAL_NOVAE
+void nova_angular_momentum_changes(
+ struct stardata_t * const stardata,
+ struct star_t * const exploder,
+ const double dm,
+ double * const dJstar,
+ double * const dJorb
+ );
+
#endif // NOVAE_PROTOTYPES_H
diff --git a/src/timestep/setup_fixed_timesteps.c b/src/timestep/setup_fixed_timesteps.c
index 102ce228084d7b68e41605a3d81bf2b7923ab03c..78147ccde401d3028512140a4a1e0eacc3d08f58 100644
--- a/src/timestep/setup_fixed_timesteps.c
+++ b/src/timestep/setup_fixed_timesteps.c
@@ -13,7 +13,7 @@ void setup_fixed_timesteps(struct stardata_t * RESTRICT const stardata)
* use this by default, although allow it
* to be turned off on the command line.
*/
- if(stardata->preferences->use_fixed_timestep[FIXED_TIMESTEP_YIELDS])
+ if(stardata->preferences->use_fixed_timestep[FIXED_TIMESTEP_YIELDS] == TRUE)
{
t->enabled = TRUE;
t->begin =
@@ -47,7 +47,7 @@ void setup_fixed_timesteps(struct stardata_t * RESTRICT const stardata)
* The test timestep outputs every 1Myr. It is disabled by default.
*/
t = &stardata->model.fixed_timesteps[FIXED_TIMESTEP_TEST];
- if(stardata->preferences->use_fixed_timestep[FIXED_TIMESTEP_TEST])
+ if(stardata->preferences->use_fixed_timestep[FIXED_TIMESTEP_TEST] == TRUE)
{
t->enabled = FALSE;
t->begin = 0.0;
@@ -70,5 +70,4 @@ void setup_fixed_timesteps(struct stardata_t * RESTRICT const stardata)
* for it in timestep.h and increment the NUMBER_OF_FIXED_TIMESTEPS
* appropriately.
*/
-
}
diff --git a/src/timestep/timestep_limits.c b/src/timestep/timestep_limits.c
index 8606f7951ace0950171f2de64288a9a39019e4a2..244fd698ad004daf085700fa58ad3d8cb3199364 100644
--- a/src/timestep/timestep_limits.c
+++ b/src/timestep/timestep_limits.c
@@ -269,7 +269,8 @@ void timestep_limits(Timestep_prototype_args)
}
#ifdef INDIVIDUAL_NOVAE
- if(stardata->preferences->individual_novae == TRUE)
+ if(star->nova == TRUE &&
+ stardata->preferences->individual_novae == TRUE)
{
/*
* When novae are switched on and off because the system
@@ -284,35 +285,49 @@ void timestep_limits(Timestep_prototype_args)
if(Is_zero(star->dtnova))
{
- star->dtnova = nova_recurrence_time(stardata,
- 1e-9,
+ star->dtnova = nova_recurrence_time(1e-9,
star->mass);
}
if(star->novastate == NOVA_STATE_TRIGGER)
{
- dtlim = 1.0;
/*
- * We know we're in a state of mass transfer
- * so we can update the recurrance time.
+ * Nova triggered: short timestep
+ * required
*/
- star->dtnova = nova_recurrence_time(stardata,
- accretion_rate,
+ dtlim = stardata->preferences->minimum_timestep;
+ star->dtnova = nova_recurrence_time(accretion_rate,
star->mass);
+
+ /*
+ * Keep the nova state either when accreting or
+ * until after nova_timeout: this keeps the timestep
+ * short enough to resolve novae for a while, just in case
+ * they restart.
+ */
+ star->nova_timeout = stardata->model.time + 3.0 * star->dtnova;
}
else if(star->novastate == NOVA_STATE_POST)
{
- dtlim = 1.0;
+ /*
+ * Just after nova: also a short timsetep
+ */
+ dtlim = stardata->preferences->minimum_timestep;
}
else
{
- dtlim = NOVAE_TIMESTEP_FACTOR
- * star->dtnova;
+ /*
+ * We have between nova explosions, try
+ * to resolve them
+ */
+ dtlim = stardata->preferences->timestep_multipliers[DT_LIMIT_NOVAE] * star->dtnova;
}
dtlim *= 1e-6; // convert to Myr
- /* limit timestep because of novae */
+ /*
+ * limit timestep because of novae
+ */
Limit_timestep(*dt,
dtlim,
star,
@@ -322,8 +337,7 @@ void timestep_limits(Timestep_prototype_args)
#endif // INDIVIDUAL_NOVAE
{
/* dtnova is in YEARS, */
- double dtnova = nova_recurrence_time(stardata,
- accretion_rate,
+ double dtnova = nova_recurrence_time(accretion_rate,
star->mass);
if(0)
diff --git a/src/timestep/timestep_set_default_multipliers.c b/src/timestep/timestep_set_default_multipliers.c
index 3a179f28afde8dbc9208ed5ceed1b1231e6eee7b..8bc5d219776bb9c0701ac8568a77cc52b0cf5314 100644
--- a/src/timestep/timestep_set_default_multipliers.c
+++ b/src/timestep/timestep_set_default_multipliers.c
@@ -57,7 +57,7 @@ void timestep_set_default_multipliers(struct preferences_t * preferences)
preferences->timestep_multipliers[DT_LIMIT_YVT] = 1.0; /* 36 */
preferences->timestep_multipliers[DT_LIMIT_MINIMUM_TIMESTEP] = 1.0; /* 37 */
preferences->timestep_multipliers[DT_LIMIT_MAXIMUM_TIMESTEP] = 1.0; /* 38 */
- preferences->timestep_multipliers[DT_LIMIT_NOVAE] = 1.0; /* 39 */
+ preferences->timestep_multipliers[DT_LIMIT_NOVAE] = 1e-2; /* 39 */
preferences->timestep_multipliers[DT_LIMIT_ARTIFICIAL_ACCRETION] = 1.0; /* 40 */
preferences->timestep_multipliers[DT_LIMIT_SN] = 1.0; /* 41 */
preferences->timestep_multipliers[DT_LIMIT_MASS_GAIN] = 1.0; /* 42 */