diff --git a/src/logging/log_every_timestep.c b/src/logging/log_every_timestep.c
index e783f1f6dc1f1d5ecff83a1a2f312f5b1f545976..b91c49d0358ae2b6f167ccb1b878e9ad11c65068 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 %d %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 %g %g\n",
stardata->model.time,
stardata->star[0].mass,
dt,
@@ -2204,7 +2204,9 @@ void log_every_timestep(struct stardata_t * RESTRICT const stardata)
(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
+ stardata->common.orbit.separation, //13
+ stardata->star[0].dmacc,
+ stardata->star[1].dmacc
);
}
#endif//INDIVIDUAL_NOVAE
diff --git a/src/novae/nova_event_handler.c b/src/novae/nova_event_handler.c
index 2b266494b64afac1ac932f5e14f9d06392c96d24..c59fda166797e7f5b821d9279501f00ab00ae84c 100644
--- a/src/novae/nova_event_handler.c
+++ b/src/novae/nova_event_handler.c
@@ -4,6 +4,11 @@
/*
* Event handling function for novae
+ *
+ * When you choose individual_novae to be on, nova
+ * derivatives (in star->derivative and stardata->model.derivative)
+ * are zero. Instead, novae are treated as instanteous events
+ * here, much like a supernova would be.
*/
Event_handler_function nova_event_handler(void * const eventp,
@@ -33,27 +38,50 @@ Event_handler_function nova_event_handler(void * const eventp,
const double dm_lost = (1.0 - nova_data->f) * exploder->dm_novaH;
const double dm_loss_rate = - dm_lost / stardata->model.dt;
+ /*
+ * Re-accretion of a fraction
+ */
+ const double dm_gain = exploder->nova_beta * dm_lost;
+ const double dm_gain_rate = +dm_gain / stardata->model.dt;
+
+
/*
* Simulate a nova in the derivative
* (although this is only really useful for logging)
*/
exploder->derivative[DERIVATIVE_STELLAR_MASS_NOVA] = dm_loss_rate;
+ exploder->derivative[DERIVATIVE_STELLAR_MASS_NOVA] = dm_gain_rate;
/*
- * No hydrogen is left on the surface of the white dwarf
+ * Remove material from the white dwarf:
+ * no hydrogen is left on the surface
+ *
+ * Here I just remove dm_lost from dm_novaH... but you
+ * could argue this material is mixed in and becomes part
+ * of the white dwarf after the nova.
*/
exploder->dm_novaH = 0.0;
+ exploder->dmacc = 0.0;
exploder->mass -= dm_lost;
exploder->novastate = NOVA_STATE_POST;
+
+ /*
+ * Accrete some on the companion
+ */
+ donor->mass += dm_gain;
/*
- * yield nova material
+ * calculate nova yield abundances
*/
double Xnova[ISOTOPE_ARRAY_SIZE];
nucsyn_set_nova_abunds(stardata,
exploder,
donor->Xenv,//THIS IS WRONG
Xnova);
+
+ /*
+ * yield material in nova explosion
+ */
nucsyn_calc_yields(stardata,
exploder,
dm_lost,
@@ -64,11 +92,35 @@ Event_handler_function nova_event_handler(void * const eventp,
YIELD_NOT_FINAL,
NUCSYN_SOURCE_NOVAE);
- /*
- * Re-accretion of a fraction
- */
- const double dm_gain = exploder->nova_beta * dm_lost;
- donor->mass += dm_gain;
+ /* should we call nucsyn remove dm here? */
+
+ if(Is_not_zero(dm_gain))
+ {
+
+ /*
+ * Put material on companion star: we put it
+ * in the accretion layer, let it get mixed
+ * on the next timestep
+ */
+ nucsyn_dilute_shell(
+ donor->dmacc,donor->Xacc,
+ dm_gain,Xnova
+ );
+ donor->dmacc += dm_gain;
+
+ /*
+ * "negative-yield" accreted material
+ */
+ nucsyn_calc_yields(stardata,
+ donor,
+ 0.0,
+ NULL,
+ dm_gain,
+ Xnova,
+ exploder->starnum,
+ YIELD_NOT_FINAL,
+ NUCSYN_SOURCE_NOVAE);
+ }
/*
* update stellar and orbital angular momenta
diff --git a/src/nucsyn/nucsyn_add_mass_to_surface.c b/src/nucsyn/nucsyn_add_mass_to_surface.c
new file mode 100644
index 0000000000000000000000000000000000000000..546450f46e00dfce9a2ca2c3afd9b9dede657620
--- /dev/null
+++ b/src/nucsyn/nucsyn_add_mass_to_surface.c
@@ -0,0 +1,30 @@
+#include "../binary_c.h"
+
+#ifdef NUCSYN
+void nucsyn_add_mass_to_surface(const Star_number k,
+ const double dm,
+ const Abundance * RESTRICT const X,
+ double * RESTRICT const ndmacc,
+ const Abundance * const * const RESTRICT Xacc,
+ Abundance ** RESTRICT newXacc)
+{
+ /*
+ * Add material to the surface of a star by putting it
+ * in its accretion layer
+ */
+ if(ndmacc[k]>VERY_TINY)
+ {
+ /* Existing accretion layer : have to mix */
+ nucsyn_dilute_shell(
+ ndmacc[k],newXacc[k],
+ dm,X
+ );
+ }
+ else
+ {
+ /* new accretion layer */
+ Copy_abundances(X,newXacc[k]);
+ }
+ ndmacc[k] += dm;
+}
+#endif//NUCSYN
diff --git a/src/nucsyn/nucsyn_prototypes.h b/src/nucsyn/nucsyn_prototypes.h
index a1c6b888f550d93c422bae72ebab142d8de2d171..a1a9f82a44916e8a5632b03091dc05167c20081f 100644
--- a/src/nucsyn/nucsyn_prototypes.h
+++ b/src/nucsyn/nucsyn_prototypes.h
@@ -617,6 +617,23 @@ Abundance nucsyn_elemental_abundance_by_number(const char * RESTRICT const eleme
void nucsyn_make_icache(struct store_t * const store);
void nucsyn_free_icache(Isotope ** icache);
+void nucsyn_add_mass_to_surface(const Star_number k,
+ const double dm,
+ const Abundance * RESTRICT const X,
+ double * RESTRICT const ndmacc,
+ const Abundance * const * RESTRICT const Xacc,
+ Abundance ** RESTRICT nXacc);
+
+void nucsyn_remove_mass_from_surface(
+#if defined DEBUG && DEBUG > 0
+ struct stardata_t * stardata,
+#endif
+ struct star_t * RESTRICT star,
+ const double dm,
+ Abundance ** RESTRICT X,
+ double * RESTRICT ndmacc,
+ Boolean *RESTRICT allocated);
+
#endif /* NUCSYN */
#endif /* NUCSYN_PROTOYPES */
diff --git a/src/nucsyn/nucsyn_remove_mass_from_surface.c b/src/nucsyn/nucsyn_remove_mass_from_surface.c
new file mode 100644
index 0000000000000000000000000000000000000000..1d61c9312d75934301c2c472175515c92dcdd300
--- /dev/null
+++ b/src/nucsyn/nucsyn_remove_mass_from_surface.c
@@ -0,0 +1,37 @@
+#include "../binary_c.h"
+
+void nucsyn_remove_mass_from_surface(
+#if defined DEBUG && DEBUG > 0
+ struct stardata_t * stardata,
+#endif
+ struct star_t * RESTRICT star,
+ const double dm,
+ Abundance ** RESTRICT Xp,
+ double * RESTRICT ndmacc,
+ Boolean * RESTRICT allocated)
+{
+ /*
+ * Remove material dm (>0.0) from the star, and
+ * set X, the abundance of the material removed,
+ * if X is non-NULL.
+ *
+ * Allocated can be NULL but only if X is NULL.
+ */
+
+ ndmacc[star->starnum] = Max(0.0, ndmacc[star->starnum] - dm);
+ if(likely(Xp!=NULL))
+ {
+ nucsyn_remove_dm_from_surface(
+#if defined DEBUG && DEBUG > 0
+ stardata,
+#endif
+ star,
+ dm,
+ Xp,
+ allocated);
+ }
+ else if(allocated != NULL)
+ {
+ *allocated = FALSE;
+ }
+}
diff --git a/src/nucsyn/nucsyn_update_abundances.c b/src/nucsyn/nucsyn_update_abundances.c
index a859686358f45a2a4db174b834c688f43dc204cb..717a590adbaedf5c9313119408fcf518edd3db4f 100644
--- a/src/nucsyn/nucsyn_update_abundances.c
+++ b/src/nucsyn/nucsyn_update_abundances.c
@@ -227,12 +227,12 @@ void Hot_function nucsyn_update_abundances(struct stardata_t * RESTRICT const st
/*
* Wind accretion, abundance calculated above
*/
- add_material(k,
- dm_wind_gain,
- Xwind_gain,
- ndmacc,
- Xacc,
- nXacc);
+ nucsyn_add_mass_to_surface(k,
+ dm_wind_gain,
+ Xwind_gain,
+ ndmacc,
+ Xacc,
+ nXacc);
/*
* RLOF accretion and non-conservative loss:
@@ -249,12 +249,12 @@ void Hot_function nucsyn_update_abundances(struct stardata_t * RESTRICT const st
&allocated_XRLOF_gain);
/* only the RLOF gain is accreted */
- add_material(k,
- dm_RLOF_gain,
- XRLOF_gain,
- ndmacc,
- Xacc,
- nXacc);
+ nucsyn_add_mass_to_surface(k,
+ dm_RLOF_gain,
+ XRLOF_gain,
+ ndmacc,
+ Xacc,
+ nXacc);
#ifdef DISCS
{
@@ -265,12 +265,12 @@ void Hot_function nucsyn_update_abundances(struct stardata_t * RESTRICT const st
if(stardata->common.ndiscs>=1 &&
stardata->common.discs[0].type==DISC_CIRCUMBINARY)
{
- add_material(k,
- dm_cbdisc_gain,
- stardata->common.discs[0].X,
- ndmacc,
- Xacc,
- nXacc);
+ nucsyn_add_mass_to_surface(k,
+ dm_cbdisc_gain,
+ stardata->common.discs[0].X,
+ ndmacc,
+ Xacc,
+ nXacc);
}
}
#endif // DISCS
@@ -278,7 +278,7 @@ void Hot_function nucsyn_update_abundances(struct stardata_t * RESTRICT const st
/*
* Wind loss
*/
- remove_material(
+ nucsyn_remove_mass_from_surface(
#if defined DEBUG && DEBUG > 0
stardata,
#endif
@@ -292,7 +292,7 @@ void Hot_function nucsyn_update_abundances(struct stardata_t * RESTRICT const st
/*
* RLOF and excretion disc loss
*/
- remove_material(
+ nucsyn_remove_mass_from_surface(
#if defined DEBUG && DEBUG > 0
stardata,
#endif
@@ -357,7 +357,7 @@ void Hot_function nucsyn_update_abundances(struct stardata_t * RESTRICT const st
/*
* Nova loss from the white dwarf
*/
- remove_material(
+ nucsyn_remove_mass_from_surface(
#if defined DEBUG && DEBUG > 0
stardata,
#endif
@@ -376,12 +376,12 @@ void Hot_function nucsyn_update_abundances(struct stardata_t * RESTRICT const st
/*
* Nova "reaccretion"
*/
- add_material(k,
- dm_novae,
- Xnovae,
- ndmacc,
- Xacc,
- nXacc);
+ nucsyn_add_mass_to_surface(k,
+ dm_novae,
+ Xnovae,
+ ndmacc,
+ Xacc,
+ nXacc);
Add_yield(0.0,NULL,
dm_novae,Xnovae,
NUCSYN_SOURCE_NOVAE);
@@ -627,67 +627,6 @@ static int id_wind_source(const struct star_t * const star,
}
#endif //NUCSYN_ID_SOURCES
-static void add_material(const Star_number k,
- const double dm,
- const Abundance * RESTRICT const X,
- double * RESTRICT const ndmacc,
- const Abundance * const * const RESTRICT Xacc,
- Abundance ** RESTRICT newXacc)
-{
- /*
- * Add material to the surface of a star by putting it
- * in its accretion layer
- */
- if(ndmacc[k]>VERY_TINY)
- {
- /* Existing accretion layer : have to mix */
- nucsyn_dilute_shell(
- ndmacc[k],newXacc[k],
- dm,X
- );
- }
- else
- {
- /* new accretion layer */
- Copy_abundances(X,newXacc[k]);
- }
- ndmacc[k] += dm;
-}
-
-static void remove_material(
-#if defined DEBUG && DEBUG > 0
- struct stardata_t * stardata,
-#endif
- struct star_t * RESTRICT star,
- const double dm,
- Abundance ** RESTRICT Xp,
- double * RESTRICT ndmacc,
- Boolean * RESTRICT allocated)
-{
- /*
- * Remove material dm (>0.0) from the star, and
- * set X, the abundance of the material removed,
- * if X is non-NULL.
- *
- * Allocated can be NULL but only if X is NULL.
- */
- ndmacc[star->starnum] = Max(0.0, ndmacc[star->starnum] - dm);
- if(likely(Xp!=NULL))
- {
- nucsyn_remove_dm_from_surface(
-#if defined DEBUG && DEBUG > 0
- stardata,
-#endif
- star,
- dm,
- Xp,
- allocated);
- }
- else if(allocated != NULL)
- {
- *allocated = FALSE;
- }
-}
#endif /* NUCSYN */
diff --git a/src/nucsyn/nucsyn_update_abundances.h b/src/nucsyn/nucsyn_update_abundances.h
index 5fa8a4db5c3d8b257b32fc95e4e7473a5178637d..1fe48b669965c1ff81f907c5fbbd2819f240be40 100644
--- a/src/nucsyn/nucsyn_update_abundances.h
+++ b/src/nucsyn/nucsyn_update_abundances.h
@@ -103,21 +103,6 @@ double coratio(double c,double o)
#endif
-static void add_material(const Star_number k,
- const double dm,
- const Abundance * RESTRICT const X,
- double * RESTRICT const ndmacc,
- const Abundance * const * RESTRICT const Xacc,
- Abundance ** RESTRICT nXacc);
-static void remove_material(
-#if defined DEBUG && DEBUG > 0
- struct stardata_t * stardata,
-#endif
- struct star_t * RESTRICT star,
- const double dm,
- Abundance ** RESTRICT X,
- double * RESTRICT ndmacc,
- Boolean *RESTRICT allocated);
#ifdef NUCSYN_ID_SOURCES