From 5de62327cf64774ec2775e9317b644d54a34e254 Mon Sep 17 00:00:00 2001
From: Robert Izzard <r.izzard@surrey.ac.uk>
Date: Thu, 18 May 2017 16:19:00 +0000
Subject: [PATCH] many more bug fixes

---
 runall.pl                                     |  4 +-
 src/RLOF/RLOF_NSNS_merger.c                   |  2 +-
 src/RLOF/RLOF_black_hole_merger.c             |  2 +-
 src/RLOF/RLOF_dynamical_comenv.c              |  2 +-
 ...OF_dynamical_transfer_from_WD_to_remnant.c |  4 +-
 src/RLOF/RLOF_orbital_angular_momentum_loss.c |  2 +-
 src/binary_macros.h                           | 69 ++++++++++++-------
 src/binary_parameters.h                       |  2 +-
 ...ar_momentum_and_eccentricity_derivatives.c |  4 +-
 .../calc_wind_loss_and_gain_rates.c           | 21 ++++--
 .../calculate_orbital_period.c                |  4 +-
 .../calculate_orbital_separation.c            |  4 +-
 .../check_for_collision_at_periastron.c       | 10 +--
 src/binary_star_functions/contact_system.c    |  2 +-
 .../gravitational_radiation_edot.c            |  8 +--
 .../gravitational_radiation_jdot.c            |  9 +--
 .../limit_accretion_rates.c                   |  2 +-
 src/binary_star_functions/mix_stars.c         |  8 +--
 .../non_conservative_angular_momentum_loss.c  |  2 +-
 src/binary_star_functions/radial_velocity_K.c |  4 +-
 .../test_for_roche_lobe_overflow.c            |  2 +-
 .../update_orbital_variables.c                |  2 +-
 src/binary_structures.h                       |  1 +
 src/code_options.h                            |  2 +-
 src/debug/debug_fprintf.c                     |  6 +-
 src/debug_macros.h                            |  2 +-
 src/evolution/evolution_prototypes.h          |  2 +-
 src/evolution/initialize_parameters.c         |  2 +-
 src/evolution/loop_evolution.c                | 22 ++++--
 .../system_is_observationally_single.c        | 28 ++++++++
 src/evolution/system_is_single.c              | 32 ---------
 src/logging/log_every_timestep.c              |  4 +-
 src/nucsyn/nucsyn_gce_log.c                   |  4 +-
 src/nucsyn/nucsyn_sn_yield.h                  |  6 +-
 src/setup/version.c                           |  2 +-
 src/stellar_structure/stellar_structure_AGB.c |  4 +-
 .../stellar_structure_HeStar.c                |  6 +-
 src/stellar_structure/stellar_structure_NS.c  |  2 +-
 src/stellar_structure/stellar_structure_WD.c  |  8 +--
 src/stellar_types.h                           |  7 +-
 src/supernovae/new_supernova.c                | 16 +++++
 src/supernovae/supernova.c                    |  1 +
 src/supernovae/supernovae_prototypes.h        |  1 +
 tbse                                          |  6 +-
 44 files changed, 194 insertions(+), 139 deletions(-)
 create mode 100644 src/evolution/system_is_observationally_single.c
 delete mode 100644 src/evolution/system_is_single.c

diff --git a/runall.pl b/runall.pl
index 2ea8f5055..9d9c10b98 100755
--- a/runall.pl
+++ b/runall.pl
@@ -2,8 +2,8 @@
 
 use strict;
 use rob_misc;
-
-foreach my $f (`ls /tmp/c_log-*.dat`)
+my $d = $ARGV[0] // '/tmp/';
+foreach my $f (`ls $d/c_log-*.dat`) 
 {
     chomp $f;
     my $x = slurp($f);
diff --git a/src/RLOF/RLOF_NSNS_merger.c b/src/RLOF/RLOF_NSNS_merger.c
index 77f6fa006..172b7f227 100644
--- a/src/RLOF/RLOF_NSNS_merger.c
+++ b/src/RLOF/RLOF_NSNS_merger.c
@@ -40,7 +40,7 @@ void RLOF_NSNS_merger(struct stardata_t *stardata,
     }
 
     struct new_stellar_structure_t * news = 
-        new_supernova(stardata,accretor,NULL);
+        new_supernova(stardata,accretor,donor,NULL);
     Dprint("new stellar structure at %p\n",news);
 
     /* 
diff --git a/src/RLOF/RLOF_black_hole_merger.c b/src/RLOF/RLOF_black_hole_merger.c
index 1a8478f6a..0869d30df 100644
--- a/src/RLOF/RLOF_black_hole_merger.c
+++ b/src/RLOF/RLOF_black_hole_merger.c
@@ -15,7 +15,7 @@ void RLOF_black_hole_merger(struct stardata_t * RESTRICT stardata,
 
     accretor->SN_type = SN_BH_BH;
     struct new_stellar_structure_t * news = 
-        new_supernova(stardata,accretor,NULL);
+        new_supernova(stardata,accretor,donor,NULL);
     Dprint("new stellar structure at %p\n",news);
     
     news->mass = accretor->mass + RLOF_orbit->dM_RLOF_transfer;
diff --git a/src/RLOF/RLOF_dynamical_comenv.c b/src/RLOF/RLOF_dynamical_comenv.c
index 770f06c0b..429daec14 100644
--- a/src/RLOF/RLOF_dynamical_comenv.c
+++ b/src/RLOF/RLOF_dynamical_comenv.c
@@ -88,7 +88,7 @@ Boolean RLOF_dynamical_comenv(struct stardata_t * RESTRICT stardata,
 	     * Reset orbital parameters as separation may have changed.
 	     */
 
-	    stardata->common.orbital_period = (stardata->common.separation/AU_IN_SOLAR_RADII)*sqrt(stardata->common.separation/(AU_IN_SOLAR_RADII*TOTAL_MASS));
+	    stardata->common.orbital_period = (stardata->common.separation/AU_IN_SOLAR_RADII)*sqrt(stardata->common.separation/(AU_IN_SOLAR_RADII*Total_mass));
 	    stardata->common.orbital_angular_frequency = TWOPI/stardata->common.orbital_period;
 
 	    action = CONTINUE_RLOF;
diff --git a/src/RLOF/RLOF_dynamical_transfer_from_WD_to_remnant.c b/src/RLOF/RLOF_dynamical_transfer_from_WD_to_remnant.c
index b8c6d8aca..196572449 100644
--- a/src/RLOF/RLOF_dynamical_transfer_from_WD_to_remnant.c
+++ b/src/RLOF/RLOF_dynamical_transfer_from_WD_to_remnant.c
@@ -57,7 +57,7 @@ void RLOF_dynamical_transfer_from_WD_to_remnant(struct stardata_t * RESTRICT sta
          */
         accretor->SN_type = SN_IA_He_Coel;
         struct new_stellar_structure_t * news = 
-            new_supernova(stardata,accretor,NULL);
+            new_supernova(stardata,accretor,donor,NULL);
         news->stellar_type = MASSLESS_REMNANT;
         stellar_structure_make_massless_remnant(news);
 #else
@@ -133,7 +133,7 @@ void RLOF_dynamical_transfer_from_WD_to_remnant(struct stardata_t * RESTRICT sta
     { 
         accretor->SN_type = SN_IA_CHAND_Coel;
         struct new_stellar_structure_t * news = 
-            new_supernova(stardata,accretor,NULL);
+            new_supernova(stardata,accretor,donor,NULL);
         stellar_structure_make_massless_remnant(news);
     }
 
diff --git a/src/RLOF/RLOF_orbital_angular_momentum_loss.c b/src/RLOF/RLOF_orbital_angular_momentum_loss.c
index 5e0d3d9fb..649f3c994 100644
--- a/src/RLOF/RLOF_orbital_angular_momentum_loss.c
+++ b/src/RLOF/RLOF_orbital_angular_momentum_loss.c
@@ -44,7 +44,7 @@ void RLOF_orbital_angular_momentum_loss(struct stardata_t * RESTRICT  stardata)
         /*
          * Specific orbital angular momentum
          */
-        double Lorb = Jorb / TOTAL_MASS;
+        double Lorb = Jorb / Total_mass;
 
         /*
          * gamma is the factor by which dm_lose_dt is multiplied to calculate
diff --git a/src/binary_macros.h b/src/binary_macros.h
index 8c0a53a76..9752ac834 100644
--- a/src/binary_macros.h
+++ b/src/binary_macros.h
@@ -70,10 +70,6 @@
 #define MEMMAP_STRING_LENGTH 82
 #endif /* MEMMAP */
 
-
-
-
-
 #define HRDIAG_FILENAME "/tmp/c_binary_hrdiag.dat"
 #include "stellar_structure/stellar_structure_macros.h"
 
@@ -84,16 +80,28 @@
 
 #include "binary_c_physical_constants.h"
 
-#define SYSTEM_IS_SINGLE (stardata->model.sgl==TRUE)
-#define SYSTEM_IS_BINARY (stardata->model.sgl==FALSE)
+/*
+ * Macro to define whether a system is single (or not)
+ */
+#define System_is_single (Boolean_(                                     \
+        (stardata->model.sgl==TRUE ||                                   \
+         stardata->star[0].stellar_type == MASSLESS_REMNANT ||          \
+         stardata->star[1].stellar_type == MASSLESS_REMNANT ||          \
+         stardata->common.eccentricity < TINY ||                        \
+         stardata->common.separation < TINY ||                          \
+         MORE_OR_EQUAL(stardata->common.eccentricity,1.0))))
+
+#define System_is_binary (((System_is_single)==TRUE) ? FALSE : TRUE)
 
+#define MIN_OBSERVATIONAL_STELLAR_MASS (0.01)
+#define MAX_OBSERVATIONAL_BINARY_ORBITAL_PERIOD (1e7*YEAR_LENGTH_IN_DAYS)
 
 /* 
  * Define an observable binary as one in which 
  * either star has a radial velocity above 
- * the threshold, and that SYSTEM_IS_BINARY is TRUE.
+ * the threshold, and that System_is_binary is TRUE.
  */
-#define OBSERVABLE_BINARY (SYSTEM_IS_BINARY && \
+#define Observable_binary (System_is_binary && \
                            ((radial_velocity_K(stardata,90,1)>(stardata->preferences->observable_radial_velocity_minimum)) || \
                             (radial_velocity_K(stardata,90,2)>(stardata->preferences->observable_radial_velocity_minimum))))
 #define OBSERVABLE_RADIAL_VELOCITY_MIN_DEFAULT 1.0
@@ -112,7 +120,7 @@
 #define NGTV -1.0
 #define NGTV2 -2.0
 
-#define TOTAL_MASS ((double)(stardata->star[0].mass+stardata->star[1].mass))
+#define Total_mass ((double)(stardata->star[0].mass+stardata->star[1].mass))
 #define K2 0.1
 #define K3 0.21
 #define MR23YR 0.4311
@@ -320,7 +328,23 @@
 #define T_TPAGB_FIRST_PULSE 13
 #define T_TMCMAX 14
 
-#define TIMESCALES_STRINGS {"MS","BGB","HE_IGNITION","HE_BURNING","GIANT_TINF1","GIANT_TINF2","GIANT_T","EAGB_TINF_1","EAGB_TINF_2","EAGB_T","TPAGB_TINF_1","TPAGB_TINF_2","TPAGB_T","TPAGB_FIRST_PULSE","T_TMCMAX"}
+#define TIMESCALES_STRINGS {                            \
+                            "MS",                       \
+                            "BGB",                      \
+                            "HE_IGNITION",              \
+                            "HE_BURNING",               \
+                            "GIANT_TINF1",              \
+                            "GIANT_TINF2",              \
+                            "GIANT_T",                  \
+                            "EAGB_TINF_1",              \
+                            "EAGB_TINF_2",              \
+                            "EAGB_T",                   \
+                            "TPAGB_TINF_1",             \
+                            "TPAGB_TINF_2",             \
+                            "TPAGB_T",                  \
+                            "TPAGB_FIRST_PULSE",        \
+                            "T_TMCMAX"                  \
+    }
 
 /* ditto for luminosity macros */
 #define L_ZAMS 1
@@ -404,9 +428,6 @@
 #define COMENV_MERGER_SPIN_METHOD_CONSERVE 1
 #define COMENV_MERGER_SPIN_METHOD_OMEGA 2
 
-
-
-
 #define REMNANT_LUMINOSITY 1e-10
 #define REMNANT_RADIUS 1e-10
 
@@ -730,10 +751,10 @@
 #define NUMBER_OF_ACCRETION_LIMITS 4
 #define ACCRETION_LIMIT_STRINGS {               \
         "Eddington",                            \
-            "Thermal",                          \
-            "Dynamical",                        \
-            "Hachisu"                           \
-}
+        "Thermal",                              \
+        "Dynamical",                            \
+        "Hachisu"                               \
+    }
 
 /*
  * Donor limitors
@@ -745,23 +766,21 @@
 #define NUMBER_OF_DONOR_LIMITS 4
 #define DONOR_LIMIT_STRINGS {                   \
         "Thermal",                              \
-            "Dynamical",                        \
-            "Envelope",                         \
-            "HG fractional envelope",           \
-            }
+        "Dynamical",                            \
+        "Envelope",                             \
+        "HG fractional envelope",               \
+    }
 
 #define NOVA_STATE_NONE 0
 #define NOVA_STATE_TRIGGER 1
 #define NOVA_STATE_POST 2
 #define NOVA_STATE_QUIESCENT 3
 #define NOVA_STATE_STRINGS {                    \
-    "None",                                     \
+        "None",                                 \
         "Triggered",                            \
         "Post-nova",                            \
         "Quiescent"                             \
-        }
-
-
+    }
 
 /* Time of tpagb start */
 #define TPAGB_start_time (star->time_first_pulse)
diff --git a/src/binary_parameters.h b/src/binary_parameters.h
index e0e661ac6..bc583a91e 100644
--- a/src/binary_parameters.h
+++ b/src/binary_parameters.h
@@ -686,7 +686,7 @@
 /*
  * Unit conversion factor in Bondi-Hoyle wind accretion 
  */
-#define ACC1 1.9067e+5
+#define BONDI_HOYLE_CONVERSION_FACTOR 1.9067e+5
 
 /*
  * acc2 is the Bondi-Hoyle wind accretion factor (3/2) = alpha_w 
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 694728231..b80df7150 100644
--- a/src/binary_star_functions/angular_momentum_and_eccentricity_derivatives.c
+++ b/src/binary_star_functions/angular_momentum_and_eccentricity_derivatives.c
@@ -33,7 +33,7 @@ void angular_momentum_and_eccentricity_derivatives(
         Zero_stellar_angmom_derivatives(star);
     }
 
-    if(SYSTEM_IS_BINARY)
+    if(System_is_binary)
     {
         /*
          * Eccentricity changes because of wind loss/gain:
@@ -160,7 +160,7 @@ void angular_momentum_and_eccentricity_derivatives(
             );
     }
 
-    if(SYSTEM_IS_BINARY)
+    if(System_is_binary)
     {
         /*
          * Orbital angular momentum loss because of 
diff --git a/src/binary_star_functions/calc_wind_loss_and_gain_rates.c b/src/binary_star_functions/calc_wind_loss_and_gain_rates.c
index e9d1fb13e..7aa2ad8ce 100644
--- a/src/binary_star_functions/calc_wind_loss_and_gain_rates.c
+++ b/src/binary_star_functions/calc_wind_loss_and_gain_rates.c
@@ -58,8 +58,11 @@ void calc_wind_loss_and_gain_rates(struct stardata_t * RESTRICT stardata,
                 Dprint("from calc_stellar_wind_loss dmr(%d)=%12.12g\n",
                        k,star->derivative[DERIVATIVE_STELLAR_MASS_WIND_LOSS]);
 
-                
-                if(POSITIVE_NONZERO(stardata->common.separation))
+
+                /*
+                 * Binary systems have wind accretion
+                 */
+                if(System_is_binary)
                 {
                     /*
                      * Calculate how much of wind mass loss from companion will be
@@ -67,7 +70,7 @@ void calc_wind_loss_and_gain_rates(struct stardata_t * RESTRICT stardata,
                      */
 
                     /* V1.10 different version of the above */
-                    vwind2 = 2.0*BETA*ACC1*star->mass/radius;
+                    vwind2 = 2.0*BETA*BONDI_HOYLE_CONVERSION_FACTOR*star->mass/radius;
 
 #ifdef USE_WIND_VELOCITY_OF_VW93_ON_TPAGB
 
@@ -109,7 +112,7 @@ void calc_wind_loss_and_gain_rates(struct stardata_t * RESTRICT stardata,
 
                     stardata->star[Other_star(k)].derivative[DERIVATIVE_STELLAR_MASS_WIND_GAIN] =
                         ivsqm*stardata->preferences->Bondi_Hoyle_accretion_factor*-stardata->star[k].derivative[DERIVATIVE_STELLAR_MASS_WIND_LOSS]*
-                        POW2((ACC1*stardata->star[Other_star(k)].mass/vwind2))/
+                        POW2((BONDI_HOYLE_CONVERSION_FACTOR*stardata->star[Other_star(k)].mass/vwind2))/
                         (2.0*Separation_squared*omv2);
                     
 #ifdef WRLOF_MASS_TRANSFER
@@ -125,10 +128,14 @@ void calc_wind_loss_and_gain_rates(struct stardata_t * RESTRICT stardata,
 
 
                     Dprint("wind accretion: acc2=%g dmt(%d)=%g [%3.2f %% of companion wind loss] (from ivsqm=%g dmr=%12.12g ac1=%12.12g mass=%g vwind2=%g sep=%g omv2=%g\n",
-                           stardata->preferences->Bondi_Hoyle_accretion_factor,Other_star(k),
+                           stardata->preferences->Bondi_Hoyle_accretion_factor,
+                           Other_star(k),
                            stardata->star[Other_star(k)].derivative[DERIVATIVE_STELLAR_MASS_WIND_GAIN],
-                           100.0*stardata->star[Other_star(k)].derivative[DERIVATIVE_STELLAR_MASS_WIND_GAIN]/MAX(1e-20,-stardata->star[k].derivative[DERIVATIVE_STELLAR_MASS_WIND_LOSS]),
-                           ivsqm,-star->derivative[DERIVATIVE_STELLAR_MASS_WIND_LOSS],ACC1,
+                           100.0*stardata->star[Other_star(k)].derivative[DERIVATIVE_STELLAR_MASS_WIND_GAIN]/
+                           MAX(1e-20,-fabs(stardata->star[k].derivative[DERIVATIVE_STELLAR_MASS_WIND_LOSS])),
+                           ivsqm,
+                           -star->derivative[DERIVATIVE_STELLAR_MASS_WIND_LOSS],
+                           BONDI_HOYLE_CONVERSION_FACTOR,
                            stardata->star[Other_star(k)].mass,
                            vwind2,stardata->common.separation,omv2);
                 }
diff --git a/src/binary_star_functions/calculate_orbital_period.c b/src/binary_star_functions/calculate_orbital_period.c
index a41610e55..60c7b7e9d 100644
--- a/src/binary_star_functions/calculate_orbital_period.c
+++ b/src/binary_star_functions/calculate_orbital_period.c
@@ -5,11 +5,11 @@ double calculate_orbital_period (struct stardata_t * RESTRICT stardata)
     /* calculate the orbital period using Kepler's law, return it in days */
 
     double per=sqrt(POW3(stardata->common.separation/
-			 (double)AU_IN_SOLAR_RADII)/(TOTAL_MASS))
+			 (double)AU_IN_SOLAR_RADII)/(Total_mass))
 	*YEAR_LENGTH_IN_DAYS;
     
     Dprint("Calculating orbital period: %g (M1+M2=%g a=%g)\n",per,
-           TOTAL_MASS,stardata->common.separation);
+           Total_mass,stardata->common.separation);
 
     return per;
 }
diff --git a/src/binary_star_functions/calculate_orbital_separation.c b/src/binary_star_functions/calculate_orbital_separation.c
index fd267fcf0..2c1f600ce 100644
--- a/src/binary_star_functions/calculate_orbital_separation.c
+++ b/src/binary_star_functions/calculate_orbital_separation.c
@@ -6,10 +6,10 @@ double calculate_orbital_separation (struct stardata_t * RESTRICT stardata)
 {
     /* calculate the orbital separation using Kepler's law */
     double sep = (double)AU_IN_SOLAR_RADII*
-        cbrt(POW2(stardata->common.orbital_period)*TOTAL_MASS);
+        cbrt(POW2(stardata->common.orbital_period)*Total_mass);
 
     Dprint("Calculating stellar separation: %g (M1+M2=%g P=%g)\n",sep,
-           TOTAL_MASS,stardata->common.orbital_period);
+           Total_mass,stardata->common.orbital_period);
 
     return sep;
 }
diff --git a/src/binary_star_functions/check_for_collision_at_periastron.c b/src/binary_star_functions/check_for_collision_at_periastron.c
index bc6288cb9..35af98216 100644
--- a/src/binary_star_functions/check_for_collision_at_periastron.c
+++ b/src/binary_star_functions/check_for_collision_at_periastron.c
@@ -22,10 +22,10 @@ Boolean check_for_collision_at_periastron(struct stardata_t * RESTRICT stardata)
            stardata->star[0].mass,
            stardata->star[1].stellar_type,
            stardata->star[1].mass,
-           SYSTEM_IS_BINARY
+           System_is_binary
         );
 
-    Boolean x = Boolean_(SYSTEM_IS_BINARY && 
+    Boolean x = Boolean_(System_is_binary && 
                          PERIASTRON_DISTANCE<TOTAL_RADIUS && 
                          stardata->model.intpol==0 );
 
@@ -40,13 +40,15 @@ Boolean check_for_collision_at_periastron(struct stardata_t * RESTRICT stardata)
     if(x==TRUE)
     {
         Append_logstring(LOG_PERIASTRON_COLLISION,
-                         "R1+R2=%g+%g=%g > aperi=%g*(1-%g)=%g",
+                         "R1+R2=%g+%g=%g > aperi=%g*(1-%g)=%g (binary? %d)",
                          stardata->star[0].radius,
                          stardata->star[1].radius,
                          TOTAL_RADIUS,
                          stardata->common.separation,
                          stardata->common.eccentricity,
-                         PERIASTRON_DISTANCE);
+                         PERIASTRON_DISTANCE,
+                         System_is_binary
+            );
     }
 #endif //FILE_LOG
 
diff --git a/src/binary_star_functions/contact_system.c b/src/binary_star_functions/contact_system.c
index b00e99625..75136c96b 100644
--- a/src/binary_star_functions/contact_system.c
+++ b/src/binary_star_functions/contact_system.c
@@ -89,7 +89,7 @@ int contact_system(struct stardata_t * RESTRICT stardata)
                 /*
                  * Reset orbital parameters as separation may have changed.
                  */
-                total_mass=TOTAL_MASS;
+                total_mass=Total_mass;
                 stardata->common.orbital_period = (stardata->common.separation/AU_IN_SOLAR_RADII)*
                     sqrt(stardata->common.separation/(AU_IN_SOLAR_RADII*total_mass));
                 stardata->common.orbital_angular_frequency = TWOPI/stardata->common.orbital_period;
diff --git a/src/binary_star_functions/gravitational_radiation_edot.c b/src/binary_star_functions/gravitational_radiation_edot.c
index 80212effc..59f6f3f97 100644
--- a/src/binary_star_functions/gravitational_radiation_edot.c
+++ b/src/binary_star_functions/gravitational_radiation_edot.c
@@ -12,11 +12,7 @@ double gravitational_radiation_edot(struct stardata_t * RESTRICT stardata)
     if(stardata->preferences->gravitational_radiation_model!=
        GRAVITATIONAL_RADIATION_NONE
        &&
-       IS_NOT_ZERO(stardata->common.separation)
-       &&
-       MORE_OR_EQUAL(stardata->common.eccentricity,0.0)
-       &&
-       stardata->common.eccentricity<1.0
+       System_is_binary
         )
     {
         if(
@@ -32,7 +28,7 @@ double gravitational_radiation_edot(struct stardata_t * RESTRICT stardata)
             /*
              * Eggleton's weak field formula
              */
-            edot = (double) -8.315e-10*MASS_PRODUCT*TOTAL_MASS/(Separation_squared*Separation_squared)
+            edot = (double) -8.315e-10*MASS_PRODUCT*Total_mass/(Separation_squared*Separation_squared)
                 *stardata->common.eccentricity*f1/POW5(sqrt(1.0-ecc2));
         }
         /*
diff --git a/src/binary_star_functions/gravitational_radiation_jdot.c b/src/binary_star_functions/gravitational_radiation_jdot.c
index 84b09aad3..89a506b04 100644
--- a/src/binary_star_functions/gravitational_radiation_jdot.c
+++ b/src/binary_star_functions/gravitational_radiation_jdot.c
@@ -15,12 +15,7 @@ double gravitational_radiation_jdot(struct stardata_t * RESTRICT stardata)
     if(stardata->preferences->gravitational_radiation_model!=
        GRAVITATIONAL_RADIATION_NONE
        &&
-       IS_NOT_ZERO(stardata->common.separation)
-       &&
-       MORE_OR_EQUAL(stardata->common.eccentricity,0.0)
-       &&
-       stardata->common.eccentricity<1.0
-        )
+       System_is_binary)
     {
         if(
             stardata->preferences->gravitational_radiation_model==GRAVITATIONAL_RADIATION_H02 
@@ -33,7 +28,7 @@ double gravitational_radiation_jdot(struct stardata_t * RESTRICT stardata)
              * Eggleton's weak field approximation
              */
             const double ecc2 = ECC_SQUARED;
-            jdot = -8.315e-10*MASS_PRODUCT*TOTAL_MASS/(Separation_squared*Separation_squared)*
+            jdot = -8.315e-10*MASS_PRODUCT*Total_mass/(Separation_squared*Separation_squared)*
                 stardata->common.orbital_angular_momentum*
                 (1.0 + 7.0/8.0 * ecc2)/POW5( sqrt(1.0 - ecc2) );
         }
diff --git a/src/binary_star_functions/limit_accretion_rates.c b/src/binary_star_functions/limit_accretion_rates.c
index 96492bb86..dd759251a 100644
--- a/src/binary_star_functions/limit_accretion_rates.c
+++ b/src/binary_star_functions/limit_accretion_rates.c
@@ -130,7 +130,7 @@ void limit_accretion_rates(struct stardata_t * RESTRICT stardata)
             }
         
             /* impose accretion limit */
-            double Lorb = JORB / TOTAL_MASS;
+            double Lorb = JORB / Total_mass;
             
             for(i=0;i<NUMBER_OF_ACCRETION_LIMITS;i++)
             {
diff --git a/src/binary_star_functions/mix_stars.c b/src/binary_star_functions/mix_stars.c
index 201ec23c9..d2b6d6430 100644
--- a/src/binary_star_functions/mix_stars.c
+++ b/src/binary_star_functions/mix_stars.c
@@ -599,7 +599,7 @@ void mix_stars(struct star_t *stars,
             Dprint("Yield SNIa m=%g\n",m1+m2);
             stars[0].SN_type = SN_IA_CHAND_Coel;
             struct new_stellar_structure_t * news = 
-                new_supernova(stardata,&stars[0],NULL);
+                new_supernova(stardata,&stars[0],NULL,NULL);
             stellar_structure_make_massless_remnant(news);
         }
     }
@@ -650,7 +650,7 @@ void mix_stars(struct star_t *stars,
 
             stars[0].SN_type = SN_TZ;
             struct new_stellar_structure_t * news = 
-                new_supernova(stardata,&stars[0],NULL);
+                new_supernova(stardata,&stars[0],NULL,NULL);
             Dprint("new stellar structure at %p\n",news);
             news->stellar_type = merged_stellar_type;
         }
@@ -682,7 +682,7 @@ void mix_stars(struct star_t *stars,
                    merged_stellar_type);
 
             struct new_stellar_structure_t * news = 
-                new_supernova(stardata,&stars[0],NULL);
+                new_supernova(stardata,&stars[0],NULL,NULL);
 
             Dprint("new stellar structure at %p\n",news);
             news->stellar_type = merged_stellar_type;
@@ -708,7 +708,7 @@ void mix_stars(struct star_t *stars,
 
         stars[0].SN_type = SN_IA_He;
         struct new_stellar_structure_t * news = 
-            new_supernova(stardata,&stars[0],NULL);
+            new_supernova(stardata,&stars[0],NULL,NULL);
         stellar_structure_make_massless_remnant(news);
         
 #ifdef NUCSYN
diff --git a/src/binary_star_functions/non_conservative_angular_momentum_loss.c b/src/binary_star_functions/non_conservative_angular_momentum_loss.c
index 460889dc1..cbdc409e4 100644
--- a/src/binary_star_functions/non_conservative_angular_momentum_loss.c
+++ b/src/binary_star_functions/non_conservative_angular_momentum_loss.c
@@ -11,7 +11,7 @@ void non_conservative_angular_momentum_loss(struct stardata_t * stardata,
      * hit or because of non-conservative RLOF.
      */
 
-    double Lorb = JORB / TOTAL_MASS;
+    double Lorb = JORB / Total_mass;
     Star_number k;
     STARLOOP(k)
     {
diff --git a/src/binary_star_functions/radial_velocity_K.c b/src/binary_star_functions/radial_velocity_K.c
index 8665e90ad..30f32a949 100644
--- a/src/binary_star_functions/radial_velocity_K.c
+++ b/src/binary_star_functions/radial_velocity_K.c
@@ -16,7 +16,7 @@ double radial_velocity_K(struct stardata_t * RESTRICT stardata,
     double K;
     Star_number nstar = n-1;
 
-    if(SYSTEM_IS_BINARY)
+    if(System_is_binary)
     {
         double aa = stardata->common.separation * R_SUN;
 
@@ -47,7 +47,7 @@ double radial_velocity_K(struct stardata_t * RESTRICT stardata,
         1.0/sqrt(stardata->common.separation * R_SUN);
 
     printf("K (binary=%d this star %d, other star %d) = sqrt(%g * %g * (%g+%g)*(1+%g)/((1-%g) * %g*%g)) = sqrt(%g) = %g (Kalt = %g)\n",
-           SYSTEM_IS_BINARY,
+           System_is_binary,
            nstar,
            Other_star(nstar),
            GRAVITATIONAL_CONSTANT,
diff --git a/src/binary_star_functions/test_for_roche_lobe_overflow.c b/src/binary_star_functions/test_for_roche_lobe_overflow.c
index fb1cf4b08..e779765a0 100644
--- a/src/binary_star_functions/test_for_roche_lobe_overflow.c
+++ b/src/binary_star_functions/test_for_roche_lobe_overflow.c
@@ -16,7 +16,7 @@ unsigned int test_for_roche_lobe_overflow(struct stardata_t * RESTRICT stardata)
     /*
      * If we're a single-star system, we cannot have RLOF
      */
-    if(!SYSTEM_IS_BINARY) return NO_ROCHE_OVERFLOW;
+    if(!System_is_binary) return NO_ROCHE_OVERFLOW;
     
     /*
      * Set model->ndonor    to the donor, 
diff --git a/src/binary_star_functions/update_orbital_variables.c b/src/binary_star_functions/update_orbital_variables.c
index 21d7835e2..c4daba24a 100644
--- a/src/binary_star_functions/update_orbital_variables.c
+++ b/src/binary_star_functions/update_orbital_variables.c
@@ -8,7 +8,7 @@ void update_orbital_variables(struct stardata_t * RESTRICT stardata)
      * and update Roche radii
      */
     
-    double total_mass=TOTAL_MASS;
+    double total_mass=Total_mass;
     double awas=stardata->common.separation;
 
     if(MORE_OR_EQUAL(stardata->common.eccentricity,1.0)||
diff --git a/src/binary_structures.h b/src/binary_structures.h
index c316aa44a..0e75ae3e9 100644
--- a/src/binary_structures.h
+++ b/src/binary_structures.h
@@ -215,6 +215,7 @@ struct tmpstore_t {
 /* new supernova struct */
 struct new_supernova_t {
     struct new_stellar_structure_t * new_stellar_structure;
+    struct new_stellar_structure_t * new_companion_structure;
 #ifdef NUCSYN
     Abundance X[ISOTOPE_ARRAY_SIZE];
 #endif
diff --git a/src/code_options.h b/src/code_options.h
index 6bd5e87b9..7d6ed1933 100644
--- a/src/code_options.h
+++ b/src/code_options.h
@@ -51,7 +51,7 @@
  *
  * See debug_macros.h for further options.
  */
-#define DEBUG 1
+#define DEBUG 0
 #define FLUSH_LOG
 
 /*
diff --git a/src/debug/debug_fprintf.c b/src/debug/debug_fprintf.c
index c5f45e00e..f119a869c 100644
--- a/src/debug/debug_fprintf.c
+++ b/src/debug/debug_fprintf.c
@@ -63,12 +63,14 @@ void debug_fprintf(struct stardata_t * RESTRICT stardata,
     {
         snprintf(sid,
                  99,
-                 "(ID %d) SN %g %g maxt=%g ",
+                 "(ID %d) SN %g %g dt=%g %g ",
                  stardata->model.id_number,
                  stardata->star[0].mass,
                  stardata->star[1].mass,
-                 stardata->model.max_evolution_time
+                 stardata->model.dt,
+                 stardata->model.dtm
             );
+
     }
     else
 #endif // BINARY_C_API
diff --git a/src/debug_macros.h b/src/debug_macros.h
index c04c2af03..3bb56e322 100644
--- a/src/debug_macros.h
+++ b/src/debug_macros.h
@@ -40,7 +40,7 @@
  * variable.
  */
 #define Debug_expression (1)
-//#define Debug_expression (stardata!=NULL && stardata->model.time>7.2)
+//#define Debug_expression (stardata!=NULL && stardata->model.time>0.0)
 
 /* enable line number output from Dprint */
 #define DEBUG_LINENUMBERS
diff --git a/src/evolution/evolution_prototypes.h b/src/evolution/evolution_prototypes.h
index e0989e925..cf3309d2b 100644
--- a/src/evolution/evolution_prototypes.h
+++ b/src/evolution/evolution_prototypes.h
@@ -35,7 +35,7 @@ void init_model(struct model_t * RESTRICT model);
 void init_common(struct stardata_t * RESTRICT stardata);
 
 void evolution_nanchecks(struct stardata_t * RESTRICT stardata);
-Boolean system_is_single(struct stardata_t * RESTRICT stardata);
+Boolean system_is_observationally_single(struct stardata_t * RESTRICT stardata);
 void zero_derivatives(struct stardata_t * RESTRICT stardata);
 void evolution_cleanup(struct stardata_t * RESTRICT stardata,
                        const Boolean free_pointers,
diff --git a/src/evolution/initialize_parameters.c b/src/evolution/initialize_parameters.c
index 0128a67a5..ec458d41c 100644
--- a/src/evolution/initialize_parameters.c
+++ b/src/evolution/initialize_parameters.c
@@ -164,7 +164,7 @@ void initialize_parameters(struct stardata_t * RESTRICT stardata)
 #endif
     {
         common->orbital_angular_frequency = TWOPI/common->orbital_period;
-        common->orbital_angular_momentum = MASS_PRODUCT/TOTAL_MASS*sqrt(1.0-ECC_SQUARED)*Separation_squared*(common->orbital_angular_frequency);
+        common->orbital_angular_momentum = MASS_PRODUCT/Total_mass*sqrt(1.0-ECC_SQUARED)*Separation_squared*(common->orbital_angular_frequency);
 
         if(star1->omega<0.0) star1->omega=common->orbital_angular_frequency;
         if(star2->omega<0.0) star2->omega=common->orbital_angular_frequency;
diff --git a/src/evolution/loop_evolution.c b/src/evolution/loop_evolution.c
index d59a7f07c..2c4cbb4fd 100644
--- a/src/evolution/loop_evolution.c
+++ b/src/evolution/loop_evolution.c
@@ -41,10 +41,24 @@ Boolean loop_evolution(struct stardata_t * RESTRICT stardata)
                 star->stellar_timestep = stardata->model.max_evolution_time;
             }
         }
-        stardata->common.eccentricity = -1.0;
-        stardata->common.separation = 0.0;
-        stardata->model.dtm = 0.0;
+        Dprint("Binary check : sep=%g ecc=%g kw=%d,%d : %d\n",
+               stardata->common.separation,
+               stardata->common.eccentricity,
+               stardata->star[0].stellar_type,
+               stardata->star[1].stellar_type,
+               System_is_binary);
 
+        /*
+         * No longer a binary
+         */
+        if(MORE_OR_EQUAL(stardata->common.eccentricity,0.0) &&
+           stardata->common.separation > TINY)
+        {
+            Dprint("system no longer binary, set sep=0, ecc=-1, dtm=0");
+            stardata->common.eccentricity = -1.0;
+            stardata->common.separation = 0.0;
+            stardata->model.dtm = 0.0;
+        }
         /* if a merger happened, set the appropriate flag */
         if(stardata->model.coel==TRUE)
         {
@@ -52,7 +66,7 @@ Boolean loop_evolution(struct stardata_t * RESTRICT stardata)
             stardata->model.merged = TRUE;
             stardata->model.coel   = FALSE;
         }
-
+            
         Dprint("m1=%g m10=%g: return TRUE, sgl=%d\n",
                stardata->star[0].mass,
                stardata->star[0].phase_start_mass,
diff --git a/src/evolution/system_is_observationally_single.c b/src/evolution/system_is_observationally_single.c
new file mode 100644
index 000000000..d32c7f321
--- /dev/null
+++ b/src/evolution/system_is_observationally_single.c
@@ -0,0 +1,28 @@
+#include "../binary_c.h"
+
+Boolean system_is_observationally_single(struct stardata_t * RESTRICT stardata)
+{
+    /*
+     * Is a system single (TRUE) or binary (FALSE) in 
+     * an observational sense.
+     */
+    Boolean single = 
+        System_is_single ||
+        // very low mass -> not a star
+        stardata->star[1].mass < MIN_OBSERVATIONAL_STELLAR_MASS ||
+        // at very long period, star will be 'apparently single' 
+        stardata->common.orbital_period >  MAX_OBSERVATIONAL_BINARY_ORBITAL_PERIOD 
+        ;
+
+    Dprint("system_is_observationally_single sgl %d : MLR1 %d : MLR2 %d  : M2<0.01 %d : P>long %d : a<0 %d : e<0 %d : single ? %d\n",
+           stardata->model.sgl,
+           stardata->star[0].stellar_type==MASSLESS_REMNANT ,
+           stardata->star[1].stellar_type==MASSLESS_REMNANT ,
+           stardata->star[1].mass < 0.010001 , 
+           stardata->common.orbital_period > 1e7*YEAR_LENGTH_IN_DAYS ,
+           stardata->common.separation < 0.0 ,
+           stardata->common.eccentricity < 0.0,
+           single);
+    
+    return single;
+}
diff --git a/src/evolution/system_is_single.c b/src/evolution/system_is_single.c
deleted file mode 100644
index 576d9c1fc..000000000
--- a/src/evolution/system_is_single.c
+++ /dev/null
@@ -1,32 +0,0 @@
-#include "../binary_c.h"
-
-Boolean system_is_single(struct stardata_t * RESTRICT stardata)
-{
-    /*
-     * Is a system single (TRUE) or binary (FALSE) 
-     */
-    Boolean single = 
-        // sgl flag 
-        stardata->model.sgl ||
-        stardata->star[0].stellar_type==MASSLESS_REMNANT ||
-        stardata->star[1].stellar_type==MASSLESS_REMNANT ||
-        stardata->star[1].mass < 0.010001 || 
-        // at very long period, star will be 'apparently single' 
-        stardata->common.orbital_period > 1e7*YEAR_LENGTH_IN_DAYS ||
-        // obviously single stars (broken binaries)
-        stardata->common.separation < 0.0 ||
-        stardata->common.eccentricity < 0.0
-        ;
-
-    Dprint("system_is_single sgl %d : MLR1 %d : MLR2 %d  : M2<0.01 %d : P>long %d : a<0 %d : e<0 %d : single ? %d\n",
-           stardata->model.sgl,
-           stardata->star[0].stellar_type==MASSLESS_REMNANT ,
-           stardata->star[1].stellar_type==MASSLESS_REMNANT ,
-           stardata->star[1].mass < 0.010001 , 
-           stardata->common.orbital_period > 1e7*YEAR_LENGTH_IN_DAYS ,
-           stardata->common.separation < 0.0 ,
-           stardata->common.eccentricity < 0.0,
-           single);
-    
-    return single;
-}
diff --git a/src/logging/log_every_timestep.c b/src/logging/log_every_timestep.c
index 287396e7f..a7408593b 100644
--- a/src/logging/log_every_timestep.c
+++ b/src/logging/log_every_timestep.c
@@ -41,7 +41,7 @@ void log_every_timestep(struct stardata_t * RESTRICT stardata)
     stardata->model.prev_log_t = stardata->model.time;
 
     /* single star condition */
-    Boolean single MAYBE_UNUSED = system_is_single(stardata); 
+    Boolean single MAYBE_UNUSED = system_is_observationally_single(stardata); 
     Boolean binary MAYBE_UNUSED = NOT(single);
 
     if(0)PRINTF("TESTLOG %30.12e %g %g %d %g %g %g\n",
@@ -2009,7 +2009,7 @@ void log_every_timestep(struct stardata_t * RESTRICT stardata)
                    stardata->star[0].pms_mass,
                    stardata->star[1].pms_mass);
 
-            if(SYSTEM_IS_BINARY)
+            if(System_is_binary)
             {
                 /* current binary conditions */
                 PRINTF("%g %g %g %g %g ",
diff --git a/src/nucsyn/nucsyn_gce_log.c b/src/nucsyn/nucsyn_gce_log.c
index ac9d743dd..7393e3dc1 100644
--- a/src/nucsyn/nucsyn_gce_log.c
+++ b/src/nucsyn/nucsyn_gce_log.c
@@ -208,7 +208,7 @@ void nucsyn_gce_log(struct stardata_t * RESTRICT stardata,
             { 
                 Dprint("Binary statistics");
 
-                ensemble_weight[OBSERVABLE_BINARY ?
+                ensemble_weight[Observable_binary ?
                                 NUCSYN_ENSEMBLE_NUMBER_OF_OBSERVABLE_BINARY_STARS :
                                 NUCSYN_ENSEMBLE_NUMBER_OF_OBSERVABLE_SINGLE_STARS]=1.0;
 
@@ -979,7 +979,7 @@ void nucsyn_gce_log(struct stardata_t * RESTRICT stardata,
             {
                 ensemble_weight[NUCSYN_ENSEMBLE_BLUE_STRAGGLER]=1.0;
             
-                if(OBSERVABLE_BINARY)
+                if(Observable_binary)
                 {
                     ensemble_weight[NUCSYN_ENSEMBLE_BINARY_BLUE_STRAGGLER]=1.0;
                 }
diff --git a/src/nucsyn/nucsyn_sn_yield.h b/src/nucsyn/nucsyn_sn_yield.h
index f451a5e54..885ed360a 100644
--- a/src/nucsyn/nucsyn_sn_yield.h
+++ b/src/nucsyn/nucsyn_sn_yield.h
@@ -18,6 +18,10 @@ static void nucsyn_sn_test(struct stardata_t *stardata);
  */
 #define Other_exploder(A) Exit_binary_c(        \
         OUT_OF_RANGE,                           \
-        "Other exploder error %d\n",A);
+        "Other exploder error %d : stellar types pre_explosion_star %d, exploder %d, companion %d\n", \
+        A,                                                              \
+        pre_explosion_star->stellar_type,                               \
+        exploder->stellar_type,                                         \
+        companion->stellar_type);
 
 #endif // NUCSYN_SN_YIELD_H
diff --git a/src/setup/version.c b/src/setup/version.c
index 820a9038c..5ddd85fe0 100644
--- a/src/setup/version.c
+++ b/src/setup/version.c
@@ -525,7 +525,7 @@ void version(struct stardata_t * RESTRICT stardata)
     PRINTF("LUM0 %g\nKAP %g\n",LUM0,KAP);
     PRINTF("AHE %g\nACO %g\nTAUMIN %g\nMLP %g\n",
            AHE,ACO,TAUMIN,MLP);
-    PRINTF("ACC1 %g\nBONDI_HOYLE_ACCRETION_FACTOR_DEFAULT %g\n",ACC1,BONDI_HOYLE_ACCRETION_FACTOR_DEFAULT);
+    PRINTF("BONDI_HOYLE_CONVERSION_FACTOR %g\nBONDI_HOYLE_ACCRETION_FACTOR_DEFAULT %g\n",BONDI_HOYLE_CONVERSION_FACTOR,BONDI_HOYLE_ACCRETION_FACTOR_DEFAULT);
 
     PRINTF("DTFAC_DEFAULT = %g\n",
            DTFAC_DEFAULT);
diff --git a/src/stellar_structure/stellar_structure_AGB.c b/src/stellar_structure/stellar_structure_AGB.c
index 1f73bb28d..a327eecca 100644
--- a/src/stellar_structure/stellar_structure_AGB.c
+++ b/src/stellar_structure/stellar_structure_AGB.c
@@ -204,7 +204,7 @@ Stellar_type stellar_structure_AGB(struct new_stellar_structure_t * s,
                     Dprint("STELLAR_STRUCTURE SN IIa\n");
                     s->SN_type = SN_IIa;
                     struct new_stellar_structure_t * news = 
-                        new_supernova(stardata,star,s);
+                        new_supernova(stardata,star,NULL,s);
                     news->stellar_type = MASSLESS_REMNANT;
                     news->age = 0.0;
                     news->mass = 0.0;
@@ -219,7 +219,7 @@ Stellar_type stellar_structure_AGB(struct new_stellar_structure_t * s,
                     Dprint("STELLAR_STRUCTURE core collapse SN\n");
                     s->SN_type = SN_II;
                     struct new_stellar_structure_t * news = 
-                        new_supernova(stardata,star,s);
+                        new_supernova(stardata,star,NULL,s);
                     news->mass = bh_mass(s->mass,
                                          s->core_mass,
                                          stardata,star);
diff --git a/src/stellar_structure/stellar_structure_HeStar.c b/src/stellar_structure/stellar_structure_HeStar.c
index c066db26b..c924a7281 100644
--- a/src/stellar_structure/stellar_structure_HeStar.c
+++ b/src/stellar_structure/stellar_structure_HeStar.c
@@ -252,7 +252,7 @@ void stellar_structure_HeStar(struct new_stellar_structure_t * s,
 
                     s->SN_type = SN_HeStarIa;
                     struct new_stellar_structure_t * news = 
-                        new_supernova(stardata,star,s);
+                        new_supernova(stardata,star,NULL,s);
                     stellar_structure_make_massless_remnant(news);
                 }
                 else
@@ -270,7 +270,7 @@ void stellar_structure_HeStar(struct new_stellar_structure_t * s,
                          * an iron core. 
                          */
                         s->SN_type = SN_ECAP;
-                        news = new_supernova(stardata,star,s);
+                        news = new_supernova(stardata,star,NULL,s);
                         news->mass = bh_mass(s->mass,s->core_mass,stardata,star);
                         news->core_mass = s->mass;
                     }
@@ -336,7 +336,7 @@ void stellar_structure_HeStar(struct new_stellar_structure_t * s,
                         {
                             s->SN_type = SN_GRB_COLLAPSAR;
                         }
-                        news = new_supernova(stardata,star,s);
+                        news = new_supernova(stardata,star,NULL,s);
                         news->mass = bh_mass(s->mass,
                                              s->core_mass, 
                                              stardata,
diff --git a/src/stellar_structure/stellar_structure_NS.c b/src/stellar_structure/stellar_structure_NS.c
index 2881514cf..342c56563 100644
--- a/src/stellar_structure/stellar_structure_NS.c
+++ b/src/stellar_structure/stellar_structure_NS.c
@@ -21,7 +21,7 @@ void stellar_structure_NS(struct new_stellar_structure_t * s,
         {
             s->SN_type = SN_AIC_BH;
             struct new_stellar_structure_t * news =
-                new_supernova(stardata,star,s);
+                new_supernova(stardata,star,NULL,s);
             news->stellar_type = BLACK_HOLE;
             news->age = 0.0;
         }
diff --git a/src/stellar_structure/stellar_structure_WD.c b/src/stellar_structure/stellar_structure_WD.c
index 752793cd0..71f7bf32b 100644
--- a/src/stellar_structure/stellar_structure_WD.c
+++ b/src/stellar_structure/stellar_structure_WD.c
@@ -34,7 +34,7 @@ void stellar_structure_WD(struct new_stellar_structure_t * s,
         if(s->stellar_type==ONeWD)
         {
             s->SN_type = SN_AIC;
-            news = new_supernova(stardata,star,s);
+            news = new_supernova(stardata,star,NULL,s);
             news->stellar_type = NEUTRON_STAR;
             Dprint("accretion induced sn with core mass %g\n",s->mass);
             news->mass = bh_mass(news->mass,news->mass,stardata,star);
@@ -49,7 +49,7 @@ void stellar_structure_WD(struct new_stellar_structure_t * s,
              */
             s->SN_type = SN_IA_CHAND;
             Dprint("Make new supernova\n");
-            news = new_supernova(stardata,star,s);
+            news = new_supernova(stardata,star,NULL,s);
             stellar_structure_make_massless_remnant(news);
             Append_logstring(LOG_ROCHE,
                              "SNIa M=%g>MCh=%g (%g %g) companion type %d",
@@ -79,7 +79,7 @@ void stellar_structure_WD(struct new_stellar_structure_t * s,
         {
             Dprint("COWD ELD");
             s->SN_type=SN_IA_ELD;
-            struct new_stellar_structure_t * news = new_supernova(stardata,star,s);
+            struct new_stellar_structure_t * news = new_supernova(stardata,star,NULL,s);
             
             Append_logstring(LOG_ROCHE,
                              "SNIa (ELD) DM=%g > %g (phase_start_mass = %g, companion type %d)",
@@ -96,7 +96,7 @@ void stellar_structure_WD(struct new_stellar_structure_t * s,
         {
             Dprint("HeWD SN");
             s->SN_type = SN_IA_He;
-            struct new_stellar_structure_t * news = new_supernova(stardata,star,s);
+            struct new_stellar_structure_t * news = new_supernova(stardata,star,NULL,s);
             Append_logstring(LOG_ROCHE,
                              "SNI? (Helium) M=%g > %g (companion type %d)\n",
                              news->mass,
diff --git a/src/stellar_types.h b/src/stellar_types.h
index 27f34dde6..941ef6143 100644
--- a/src/stellar_types.h
+++ b/src/stellar_types.h
@@ -130,11 +130,12 @@
 #define IS_A_PLANET(A) (stardata->star[(A)].stellar_type==LMMS && stardata->star[(A)].mass<MAXIMUM_PLANET_MASS))
 #define WAS_A_STAR(A) (stardata->star[(A)].pms_mass>MAXIMUM_PLANET_MASS)
 #define MASSLESS(A) (stardata->star[(A)].stellar_type == MASSLESS_REMNANT)
-#define EVOLVING(A) (!MASSLESS(A))
-#define IS_A_STAR(A) (EVOLVING(A) && stardata->star[(A)].mass>MAXIMUM_PLANET_MASS)
+#define EVOLVING(A) (!(MASSLESS(A)))
+#define IS_A_STAR(A) ((EVOLVING(A))\
+                      && (stardata->star[(A)].mass>MAXIMUM_PLANET_MASS))
 
 
-#define NUCLEAR_BURNING_BINARY (NUCLEAR_BURNING(stardata->star[0].stellar_type) || NUCLEAR_BURNING(stardata->star[1].stellar_type<HeWD))
+#define NUCLEAR_BURNING_BINARY ((NUCLEAR_BURNING(stardata->star[0].stellar_type)) || (NUCLEAR_BURNING(stardata->star[1].stellar_type<HeWD)))
 
 
 /* Giant filter for WR candidates */
diff --git a/src/supernovae/new_supernova.c b/src/supernovae/new_supernova.c
index bcaa34294..7d7630722 100644
--- a/src/supernovae/new_supernova.c
+++ b/src/supernovae/new_supernova.c
@@ -17,10 +17,14 @@
  * This routine allocates memory to store the appropriate structures,
  * and sets up the new structure ready to be modified by the
  * stellar structure routine.
+ *
+ * companion stores the companion stellar structure in cases where
+ * this is required. It is ignored if it is NULL.
  */
 
 struct new_stellar_structure_t * new_supernova(struct stardata_t * stardata,
                                                struct star_t * star,
+                                               struct star_t * companion,
                                                struct new_stellar_structure_t * s)
 {
     Dprint("NEW SUPERNOVA star %d t=%g dt=%g : new_supernova was allocated as %p, s=%p\n",
@@ -35,6 +39,7 @@ struct new_stellar_structure_t * new_supernova(struct stardata_t * stardata,
         Dprint("malloc new stellar structure");
         star->new_supernova = MALLOC(sizeof(struct new_supernova_t));
         star->new_supernova->new_stellar_structure = MALLOC(sizeof(struct new_stellar_structure_t));
+        star->new_supernova->new_companion_structure = NULL;
     }
     else
     {
@@ -64,6 +69,17 @@ struct new_stellar_structure_t * new_supernova(struct stardata_t * stardata,
                sizeof(struct new_stellar_structure_t));
     }
 
+    if(companion != NULL)
+    {
+        /*
+         * Add details about the companion
+         */
+        star->new_supernova->new_companion_structure = MALLOC(sizeof(struct new_stellar_structure_t));
+        set_stellar_structure_struct_from_star(
+            star->new_supernova->new_stellar_structure,
+            companion);
+    }
+    
     Dprint("Made new SN %p -> new structure %p\n",
            star->new_supernova,
            star->new_supernova->new_stellar_structure);
diff --git a/src/supernovae/supernova.c b/src/supernovae/supernova.c
index 201d95298..aed007cc8 100644
--- a/src/supernovae/supernova.c
+++ b/src/supernovae/supernova.c
@@ -200,6 +200,7 @@ void supernova(struct stardata_t * stardata,
      * Free the memory associated with the new_supernova
      */
     Safe_free(star->new_supernova->new_stellar_structure);
+    Safe_free(star->new_supernova->new_companion_structure);
     Safe_free(star->new_supernova);
     Safe_free(pre_explosion_star);
 }
diff --git a/src/supernovae/supernovae_prototypes.h b/src/supernovae/supernovae_prototypes.h
index 24218300e..588634125 100644
--- a/src/supernovae/supernovae_prototypes.h
+++ b/src/supernovae/supernovae_prototypes.h
@@ -61,6 +61,7 @@ double ejection_speed(struct star_t * star);
 
 struct new_stellar_structure_t * new_supernova(struct stardata_t * stardata,
                                                struct star_t * star,
+                                               struct star_t * companion,
                                                struct new_stellar_structure_t * s);
     
 #endif // SUPERNOVA_PROTOTYPES_H
diff --git a/tbse b/tbse
index 4f0ef93cf..a6b03441c 100755
--- a/tbse
+++ b/tbse
@@ -58,13 +58,13 @@ STELLAR_TYPE2=1
 
 # stellar masses (solar units)
 M1=7
-M2=.1
+M2=6
 
 # orbit: 
 # If the period (days) is given, use it. 
 # If the period is zero, use the separation (in Rsun).
 # Eccentricity should be in the range 0.0-1.0.
-PER=20000
+PER=2
 SEP=0
 ECC=0.0
 
@@ -510,7 +510,7 @@ DONOR_LIMIT_DYNAMICAL_MULTIPLIER=1.0
 #
 #
 # Eddington rate for accretion (was 'eddfac')
-ACCRETION_LIMIT_EDDINGTON_MULTIPLIER=1.0
+ACCRETION_LIMIT_EDDINGTON_MULTIPLIER=1.0e6
 
 # Thermal limit for accretion onto a MS, HG or CHeB star
 ACCRETION_LIMIT_THERMAL_MULTIPLIER=1.0
-- 
GitLab