diff --git a/runall.pl b/runall.pl
index 2ea8f5055f2f42fd122fa26987a03a7f511a709c..9d9c10b986875597a938074b59097c8805e4eabf 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 77f6fa006139994b77ae853ec7ca2c9008c58431..172b7f227333b129b4ff39e56e92dba0f3f22f0a 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 1a8478f6a07b3f55d69c44b64fd56ba83eb9e32e..0869d30dfd890a07577226c1b775cfb27320561d 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 770f06c0bc8b7c22d8ba2480dbb33bdcb88109a2..429daec148bf183a83fc79e66040db03b24db2e9 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 b8c6d8aca19c9a9e08a556783be5e037b9949817..196572449c14e50705a255fcc63c95dfcdcb0777 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 5e0d3d9fb4481faa993e8aff8938a474f773aae0..649f3c994c874277e346da20c6736aa94b578968 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 8c0a53a7682d1a0e384edeacdba43f4f0fd0ae99..9752ac83486fe6df80f71588e21a6dd32e19170d 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 e0e661ac6f5ab6aaad426479ee73ab07534b5bfe..bc583a91e3ea413236ab19d11305f53c227f39b8 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 6947282314b1050056e6184a739628059de3a154..b80df7150b9db54d319b3bfac3026bc2cf52a5ca 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 e9d1fb13efbabaf8e9d7e135691986940462e91d..7aa2ad8ce4f8d6624f181c6d4dea13aa227e0036 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 a41610e55749a65ddb568f1f68a1a63e054d31e6..60c7b7e9d54c0933e0b0c62c80bee1fa0c0731b3 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 fd267fcf0d3bbf435140dfebce00cd0fafb1a4bc..2c1f600ce44d9cfd208b3c3b6c1911324c7e1280 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 bc6288cb9e41a45f6fb1789ad25daa32fd0cfc62..35af98216204e5cd293eaeb4e99b72ddc9a4a692 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 b00e99625f26f119eeb31c22364b4b310de28c66..75136c96bb40f43cdd9ed01ca3798afaa4505ae2 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 80212effcff9114aa0413751695528f25ce4429d..59f6f3f97a930b450f0e89c4d58010d637af7e5d 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 84b09aad35adc453c304123bc5eaaf7240de2334..89a506b04b1abca7ee41a2a8c38908c3540cd66a 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 96492bb86985f6535ca36f737b861410245bac53..dd759251aa10101614e2da8d761ade8f5b95c1b8 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 201ec23c907a7983afa3813a47839698ed5621d9..d2b6d6430fd1cad29f047dded5946c53ad293f68 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 460889dc1be7208419f27a52f9a7ffb7769bc661..cbdc409e496f8fdba7973c8ac0399a5061aaf3fe 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 8665e90ad861912661937cbb6b6fa9554cec8374..30f32a949e243e9db303173cb10a23f82a43be94 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 fb1cf4b08fbefb96b9c21e77615ef349b8507888..e779765a0af9c855ffee5cff44a8985a6d57c50d 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 21d7835e254ea9311e70208043bfbd99f7503df2..c4daba24a5d92c6d910e7350b7349b37f0cc6f6e 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 c316aa44a465f391bb91df2801d598606e1913a5..0e75ae3e933beebb210a6a49269baf90b0d7617e 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 6bd5e87b93d98a15a8ba37950060639eacbbb071..7d6ed19336550b172849b42193975b4359616f2f 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 c5f45e00eb742ab57fb49badfdb48a872f33137a..f119a869c933537012aee941497e82d6ec97b78e 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 c04c2af0396f7e05bda2ff095f919d1043f75b0e..3bb56e32262bed14c4c8f5afc007269a245e1af8 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 e0989e925575d53cca1bf80f95230900216bf785..cf3309d2b37662c5e44227fec7f5ebedfaef94c6 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 0128a67a5668247456c57b89dbb47e64158cb7bb..ec458d41cf1ac0bf91792b48cab6b59a2341a553 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 d59a7f07c07cfff0c3289529b79eb9460e7f6ad7..2c4cbb4fdc2ab76db91ce8f4f243f202abe50f14 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 0000000000000000000000000000000000000000..d32c7f321dc06d5f328d54326ebfee99b8b0ea90
--- /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 576d9c1fcbbd8fc78d0227d12d7b20f05684a9a9..0000000000000000000000000000000000000000
--- 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 287396e7f813de4c6e95f0672b7aac7d69db764b..a7408593bed27ae4e4965910b18a4dcdc7370a78 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 ac9d743ddbfa1ed90815eb1c64109e786bd02ad3..7393e3dc1a6b25585089c4cd6c292c6bbbc5c28f 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 f451a5e540c9240fdb7ca51abc9dc389349cc773..885ed360aede31a01bbd3f09ced80233cbb83e84 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 820a9038c819772589128f4d57cd01f0f889fcfa..5ddd85fe0e5c8d87a090e4e5c5a7027c2a05b6cc 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 1f73bb28d16499ddd24ee6fc98963eed55a17d88..a327eecca1562bb23c31c58bf506e61aab60b758 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 c066db26bee28fb685e0c07596c1c6a9dea1368b..c924a7281dc8d620e14e164e05fc294bc25cc40e 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 2881514cf5c0ae1ba2cfc687a3889234745bb1c3..342c5656313a681c7835ba90a7871a275ef0bd78 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 752793cd0cee91cc617bad4d9a38ffdba3a07d13..71f7bf32b86188d8bd060865c673ef9fdc79124b 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 27f34dde68d197f10795027c489fc8afc12aae79..941ef614370510dc1fdaae6f6d5f5be2713082f9 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 bcaa342943b7e268466170af2c2d26ea06344fcd..7d76307226dc81515311f449a43aef3c6d8f858b 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 201d95298e65ec75f66ffb9aa765e621f140f7f6..aed007cc8d46ae51e69183ede53f716bfdeeba22 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 24218300e984a13ef38e5976d26093a451bc0f8d..588634125a433e92a81cf5eb1798320178562ac2 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 4f0ef93cf11e4ea0e11e21bb6af78e51c859c0ba..a6b03441c723f0d9f1e0a547dc74d6f8ec62be48 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