diff --git a/apitest/Make_fortran b/apitest/Make_fortran
index 77a00ea003d37b63ec11b5887d7ee4b6861ac903..1c82131091e68cac4a17003af7d859b78779dc40 100755
--- a/apitest/Make_fortran
+++ b/apitest/Make_fortran
@@ -16,7 +16,7 @@ F77="gfortran"
 
 LIBS=`../binary_c-config libs`
 LIBDIRS=`../binary_c-config libdirs`
-CFLAGS="-O3"
+CFLAGS=`../binary_c-config cflags`
 INCDIRS=`../binary_c-config incdirs`
 
 CMD="gfortran -lbinary_c $INCDIRS -I../src/ -L../src/ $LIBDIRS $LIBS $CFLAGS -ffree-line-length-0 apitest.f90 ../src/libbinary_c.so  -o apitest-f90"
diff --git a/apitest/apitest.c b/apitest/apitest.c
index de13fcd8c87b9939ec955af615904e193b05c827..92a7976b67d502d1f8b0d8020ec50ea2e62aa80d 100644
--- a/apitest/apitest.c
+++ b/apitest/apitest.c
@@ -54,7 +54,6 @@
 #define SHOW_BUFFER TRUE
 #define MAX_EVOLUTION_TIME 15000.0
 #define TIMESTEP 10.0
-#define ARGSTRING_LENGTH 4096
 
 /* local function prototypes */
 static double randd(double min, double max);
@@ -92,13 +91,12 @@ int main ( int argc, char * *  argv )
     /* memory for N binary systems */
     struct libbinary_c_stardata_t *stardata[N]; 
     struct libbinary_c_store_t * store = NULL;
-    char * argstring = malloc(sizeof(char)*ARGSTRING_LENGTH);
+    char * argstring;
     for(i=0;i<N;i++)
     {
         double m1 = exp(randd(log(0.1),log(100.0)));
         
-        int n_writ = snprintf(argstring,
-                              ARGSTRING_LENGTH-1,
+        int n_writ = asprintf(&argstring,
                               "binary_c M_1 %g M_2 %g eccentricity %g metallicity %g max_evolution_time %g orbital_period %g",
                               m1, // m2
                               randd(0.1,m1), // 0.1 < m2 < m1
@@ -114,13 +112,6 @@ int main ( int argc, char * *  argv )
                     n_writ);
             exit(BINARY_C_WRITE_FAILED);
         }
-        else if(n_writ >= ARGSTRING_LENGTH-1)
-        {
-            fprintf(stderr,
-                    "snprintf into argstring failed because buffer is full (n_writ = %d)\n",
-                    n_writ);
-            exit(BINARY_C_WRITE_FAILED);
-        }
 
         printf("argstring : %s\n",argstring);
         fflush(stdout);
@@ -151,8 +142,9 @@ int main ( int argc, char * *  argv )
 
         stardata[i]->preferences->internal_buffering = INTERNAL_BUFFERING_STORE;
         stardata[i]->preferences->batchmode = BATCHMODE_LIBRARY;
+        free(argstring);
     }
-    free(argstring);
+
 
     /* find available system RAM (bytes) */
     long int RAM_available=(long int)getTotalSystemMemory();
@@ -290,12 +282,13 @@ static void APIprintf(char * format,...)
     va_start(args,format);
 
     /* s contains the message */
-    static char s[10000];
-    vsnprintf(s,10000,format,args);
+    char * s;
+    asprintf(&s,format,args);
     chomp(s); // from binary_c
     fprintf(stderr,"APITEST : %s\n",s);
     fflush(stderr);
     va_end(args);
+    free(s);
 #endif
 }
 
diff --git a/apitest/apitest.f90 b/apitest/apitest.f90
index 333ead645db44d0b6c080a207a1703e0d2eba7c0..16d89c2497a6a21b4eabd0e4bbf18e1b4c4b4016 100644
--- a/apitest/apitest.f90
+++ b/apitest/apitest.f90
@@ -22,11 +22,12 @@ program apitestf
   ! These are passed to binary_c to set up the stars
   ! and/or returned by binary_c
   double precision m1,m2,mc1,mc2,per,ecc,metallicity
-  double precision t,maxt,dt,binary_c_time
+  double precision t,newt,maxt,dt,binary_c_time
   double precision omega1,omega2,r1,r2,vrot1,vrot2
   double precision logg1,logg2,lum1,lum2,teff1,teff2
   double precision xh_1,xhe_1,xc_1,xn_1,xo_1,xfe_1
   double precision xh_2,xhe_2,xc_2,xn_2,xo_2,xfe_2
+  double precision timediff
   integer stellar_type1, stellar_type2, av,ac
   character (len=1024) :: argstring
   character (len=1024) :: format_string
@@ -62,7 +63,7 @@ program apitestf
   ecc = 0.d0
   metallicity = 0.02
   maxt=15000d0
-  format_string = '( "binary_c M_1 ",E10.3," M_2 ",E10.3," metallicity ",E10.3," orbital_period ",E10.3," eccentricity ",E10.3," max_evolution_time ",E10.3,A)'
+  format_string = '("binary_c reset_stars M_1 ",E10.3," M_2 ",E10.3," metallicity ",E10.3," separation 0 orbital_period ",E10.3," eccentricity ",E10.3," max_evolution_time ",E10.3,A)'
   write(argstring,format_string) m1,m2,metallicity,per,ecc,maxt,char(0)
   
   ! create new stardata structure
@@ -93,8 +94,9 @@ program apitestf
 
      ! up our age counter
      t = t + dt
-
-     write(*,*)'At time ',t,binary_c_time,' M1=',m1,' Mc1=',Mc1,&
+     
+     write(*,*)'At time ',t,binary_c_time,&
+          ' M1=',m1,' Mc1=',Mc1,&
           'R=',r1,'omega=',omega1,&
           'logg=',logg1,'lum=',lum1,'teff1=',teff1,&
           ' stellar_type=',stellar_type1,&
diff --git a/apitest/apitest.sh b/apitest/apitest.sh
index 7ed5076df315240f7e2803a549cdccf1390319e2..7d105b9e90b9792d61dd753a9eb5200b1bff24ac 100755
--- a/apitest/apitest.sh
+++ b/apitest/apitest.sh
@@ -1,8 +1,7 @@
 #!/bin/bash
 
-export LD_LIBRARY_PATH=`../binary_c-config libdirs_list | tr " " ":"`.
-LD_LIBRARY_PATH="$LD_LIBRARY_PATH:../src"
-echo "LD_LIBRARY_PATH=$LD_LIBRARY_PATH"
+LD_LIBRARY_PATH="$(../binary_c-config libdirs_list | tr ' ' ':')"
+export LD_LIBRARY_PATH="$LD_LIBRARY_PATH"":../src"
 
 #gdb
 ./binary_c-apitest 2>&1
diff --git a/meson.build b/meson.build
index eaa47e2cc5ae36e771a8c1cce1572c6c42f66ea3..3eb318605e2ac08e99701396842e517f227e23ae 100644
--- a/meson.build
+++ b/meson.build
@@ -427,7 +427,8 @@ foreach libname : _required_libraries
             gsl_libdirs = []
             foreach gsl_libdir :run_command('sh','-c','gsl-config  --libs').stdout().strip().split(' ')
                 if gsl_libdir.startswith('-L')
-                    gsl_libdirs += gsl_libdir
+                    # add libdir without -L
+                    gsl_libdirs += gsl_libdir.strip().split('-L')[1]
                 endif
             endforeach
             libdirs += gsl_libdirs
@@ -614,7 +615,9 @@ cc_quoted = '-DCC="' + compiler.get_id() + '"'
 ld_quoted = '-DLD="' + compiler.get_id() + '"'
 incdirs_quoted = '-DINCDIRS=' + ''.join(['"-I', ' -I'.join(absolute_incdirs),'"']) # make -I... -I...
 incdirs_quoted = '_slash_'.join(incdirs_quoted.split('/')) # deslash (convert / to _slash_)
-libdirs_quoted = '-DLIBDIRS=' + ''.join(['"-L', ' -L'.join(absolute_libdirs),'"']) # make -L... -L...
+libdirs_quoted = '-DLIBDIRS=' + ''.join(['"-L', ' -L'.join(libdirs),'"']) # make -L... -L...
+message('LIBDIRS=' + libdirs_quoted)
+message('abs ' + ' '.join(libdirs))
 libdirs_quoted = '_slash_'.join(libdirs_quoted.split('/')) # deslash
 libs_quoted = '-DLIBS=' + ''.join(['"', ' '.join(libs),'"']) 
 libs_quoted = '_slash_'.join(libs_quoted.split('/'))
@@ -646,6 +649,13 @@ foreach _arg : [
 endforeach
 cflags += warn_flags
 
+############################################################
+# append library locations to cflags
+#
+foreach libdir : libdirs
+    _l = '-L' + libdir
+    cflags += _l
+endforeach
 
 ############################################################
 # make precompiled headers (PCH)
diff --git a/src/API/binary_c_API_fortran.c b/src/API/binary_c_API_fortran.c
index f73f4b8fb0cf1e76be70f2f981112468457fc48c..0bda0e6372607ccc91f64cee6ac558855ef905d2 100644
--- a/src/API/binary_c_API_fortran.c
+++ b/src/API/binary_c_API_fortran.c
@@ -93,7 +93,6 @@ void binary_c_API_function binary_c_fortran_api_stardata_info_(FORTRAN_OUT_VARLI
     *logg1 = logg(&(stardata->star[0]));
     *lum1 = stardata->star[0].luminosity;
     *teff1 = Teff(0);
-    stardata->star[0].angular_momentum=0.0;
 
     /* star 2 structure */
     *m2 = stardata->star[1].mass;
@@ -104,7 +103,6 @@ void binary_c_API_function binary_c_fortran_api_stardata_info_(FORTRAN_OUT_VARLI
     *logg2 = logg(&(stardata->star[1]));
     *lum2 = stardata->star[1].luminosity;
     *teff2 = Teff(1);
-    stardata->star[1].angular_momentum=0.0;
 
 #ifdef NUCSYN
     /* star 1 abundances */
diff --git a/src/API/binary_c_API_logging.c b/src/API/binary_c_API_logging.c
index 4e300375d8ebe3624aca1bab01a17b4f6ecbad1f..ffbb4d58695c80a7f45dc3e6c80b41ec469efc47 100644
--- a/src/API/binary_c_API_logging.c
+++ b/src/API/binary_c_API_logging.c
@@ -34,8 +34,8 @@ void binary_c_API_function binary_c_API_log(struct stardata_t * const stardata,
                       stardata->star[1].hybrid_HeCOWD,
                       stardata->common.orbit.separation,
                       stardata->common.orbit.eccentricity,
-                      stardata->star[0].radius/stardata->star[0].roche_radius,
-                      stardata->star[1].radius/stardata->star[1].roche_radius,
+                      Is_not_zero(stardata->star[0].roche_radius) ? (stardata->star[0].radius/stardata->star[0].roche_radius) : 0.0,
+                      Is_not_zero(stardata->star[1].roche_radius) ? (stardata->star[1].radius/stardata->star[1].roche_radius) : 0.0,
                       locations[loc],
                       stardata);
     fflush(stardata->model.api_log_fp);
diff --git a/src/binary_c_debug.h b/src/binary_c_debug.h
index adee6597609a9e2b2f9a9b8da2bee0d6ca25358a..4a417bf35bc6db6dc54bb1b3a0deb1d2cf24c5f1 100644
--- a/src/binary_c_debug.h
+++ b/src/binary_c_debug.h
@@ -100,7 +100,12 @@
  * this is not shown.
  */
 
-#define Debug_show_expression "  "
+#define Debug_show_expression " #=%ld Jorb=%g a=%g P=%g omegaorb=%g ",\
+        (long int)stardata->model.model_number,                       \
+        stardata->common.orbit.angular_momentum,                      \
+        stardata->common.orbit.separation,\
+        stardata->common.orbit.period,\
+        stardata->common.orbit.angular_frequency
 
 //#undef Debug_show_expression
 
diff --git a/src/binary_c_parameters.h b/src/binary_c_parameters.h
index 2c67ee6d95d4343a25aa90a74b47ff204f802ffc..9de8b52d3af1ff17621c791b8ca5f80bbdd78039 100644
--- a/src/binary_c_parameters.h
+++ b/src/binary_c_parameters.h
@@ -1522,9 +1522,6 @@
  */
 //#define FABIAN_IMF_LOG
 
-// reproduce Gijs' terminal mass plot
-//#define NELEMANS_MASS_PLOT
-
 // Stuff for Rob Detmers
 //#define DETMERS_LOG
 
diff --git a/src/evolution/evolve_system.c b/src/evolution/evolve_system.c
index 69484471ce41754d7e3ff4980f3a3e3866b3a23f..e66720d975876d7d06a10d330b33136be8b5c8cc 100644
--- a/src/evolution/evolve_system.c
+++ b/src/evolution/evolve_system.c
@@ -21,7 +21,12 @@
 
 int evolve_system(struct stardata_t * RESTRICT const stardata)
 {
-    pre_evolution(stardata);
+    
+    /* first time setup */
+    if(stardata->model.model_number == 0)
+    {
+        pre_evolution(stardata);
+    }
     
     stardata->evolving = TRUE;
 
diff --git a/src/evolution/initialize_parameters.c b/src/evolution/initialize_parameters.c
index a131d9b97102ce076a0e5a03c76b0485975bfa38..a51fd0182bc8795ab85ba87c410cbab7ec14ced5 100644
--- a/src/evolution/initialize_parameters.c
+++ b/src/evolution/initialize_parameters.c
@@ -100,12 +100,17 @@ void initialize_parameters(struct stardata_t * RESTRICT const stardata)
 #endif
     {
         common->orbit.angular_frequency = TWOPI/common->orbit.period;
+        Dprint("calc Jorb from a = %g, omega = %g\n",
+               common->orbit.separation,
+               common->orbit.angular_frequency);
         common->orbit.angular_momentum = MASS_PRODUCT/Total_mass*sqrt(1.0-ECC_SQUARED)*Separation_squared*common->orbit.angular_frequency;
-
+        
         if(star1->omega<0.0) star1->omega = common->orbit.angular_frequency;
         if(star2->omega<0.0) star2->omega = common->orbit.angular_frequency;
     }
 
+    Dprint("hence Jorb = %g\n",common->orbit.angular_momentum);
+
     const Boolean init_sgl = model->sgl;
 
     /**** per star setup ****/
@@ -519,5 +524,6 @@ void initialize_parameters(struct stardata_t * RESTRICT const stardata)
 #endif //NUCSYN
 
     Dprint("Entering first timestep\n");
+    
 }
 
diff --git a/src/evolution/pre_evolution.c b/src/evolution/pre_evolution.c
index 1b1f0506902e218adf78295ea1e3e4f55d761818..640d00101f82e9b29435c31aeeaa9015e9255939 100644
--- a/src/evolution/pre_evolution.c
+++ b/src/evolution/pre_evolution.c
@@ -9,14 +9,10 @@ void pre_evolution(struct stardata_t * RESTRICT const stardata)
     Boolean force_log = FALSE;
     if(strlen(stardata->preferences->stardata_load_filename))
     {
-        //char ** argv = stardata->common.argv;
-        //int argc = stardata->common.argc;
         load_stardata(stardata,
                       stardata->preferences->stardata_load_filename,
                       FALSE);
         stardata->model.log_fp = NULL;
-        //stardata->model.dt = stardata->model.dtm = 0.0;
-        //parse_arguments(1,argc,argv,stardata);
         force_log = TRUE;
     }
 
@@ -33,7 +29,6 @@ void pre_evolution(struct stardata_t * RESTRICT const stardata)
            stardata->model.log_fp);
     open_log_files(&(stardata->model.log_fp),stardata);   
 
-
 /* Set up the detailed log */
 #ifdef DETAILED_LOG
     /* Blank line! */
diff --git a/src/evolution/start_of_evolution.c b/src/evolution/start_of_evolution.c
index a57434bc4578686a0319d770bba2c2446438b104..9c1df139acd54fcd70fe4d54937c2c65a54f0ea3 100644
--- a/src/evolution/start_of_evolution.c
+++ b/src/evolution/start_of_evolution.c
@@ -6,6 +6,9 @@
  *
  * Note that the system does not have to have time == zero, 
  * it could start evolving half way through its lifetime.
+ *
+ * We only call initialize_parameters if model_number == 0
+ * (i.e. on the first timestep).
  */
 
 void start_of_evolution(struct stardata_t *stardata)
@@ -27,13 +30,14 @@ void start_of_evolution(struct stardata_t *stardata)
     
     /* 
      * Initialize stardata parameters
+     * for first model.
      */
-    Dprint("calling initialize_parameters\n");
-    initialize_parameters(stardata);
-
-#ifdef NELEMANS_MASS_PLOT
-    printf("GIJS0 %g\n",stardata->star[0].mass);
-#endif
+    if(stardata->model.model_number == 0)
+    {
+        Dprint("calling initialize_parameters\n");
+        initialize_parameters(stardata);
+    }
+    
    
 #ifdef NO_IMMEDIATE_MERGERS
     stardata->model.evolution_number = 0;
diff --git a/src/setup/parse_arguments.c b/src/setup/parse_arguments.c
index 98debb16e27b12e66986d0521e647a1d431e5651..1ce6238a3e31fc29caa6dfd2548484d4c52af229 100644
--- a/src/setup/parse_arguments.c
+++ b/src/setup/parse_arguments.c
@@ -540,6 +540,10 @@ static void derived_arguments(struct stardata_t * const stardata)
     stardata->star[1].stellar_type = stardata->preferences->initial_stellar_type2;
     
     stardata->common.orbit.angular_momentum = 0.0;
+    Dprint("masses %g %g\n",
+           stardata->star[0].mass,
+           stardata->star[1].mass);
+
     if(Is_not_zero(stardata->star[0].mass) &&
        Is_not_zero(stardata->star[1].mass))
     {
diff --git a/src/setup/version.c b/src/setup/version.c
index 00a1fa44b0655e1cf0998a0305115eb52b6b40d1..fb9c5c8ac0d634a7a61e3c0714b28e6824aa0e3c 100644
--- a/src/setup/version.c
+++ b/src/setup/version.c
@@ -1653,8 +1653,6 @@ void version(struct stardata_t * RESTRICT const stardata)
                 "NS/BH AIC log enabled");
     Macrotestif(NS_NS_BIRTH_LOG
                 "NS/NS birth log enabled");
-    Macrotestif(NELEMANS_MASS_PLOT,
-                "Mass plot to compare to Gijs' models is enabled");
     Macrotest(RS_LOG_ACCREITON_LOG);
     Macrotest(DETMERS_LOG);
     Macrotest(LOG_HERBERT);
diff --git a/src/stellar_structure/stellar_structure_BSE.c b/src/stellar_structure/stellar_structure_BSE.c
index f9f98fa5170f0c4641aaa5b1d2f2ad8184c32cdc..0a6ae5527c71537523bd635a8fc94326f56c27ff 100644
--- a/src/stellar_structure/stellar_structure_BSE.c
+++ b/src/stellar_structure/stellar_structure_BSE.c
@@ -32,6 +32,7 @@
 #include "../binary_c.h"
 #ifdef BSE
 #include "stellar_structure_debug.h"
+#define NEGATIVE_AGE_CHECK_MAGNITUDE 1e-10
 
 void stellar_structure_BSE(struct star_t *newstar,
                            struct star_t *oldstar,
@@ -84,7 +85,7 @@ void stellar_structure_BSE(struct star_t *newstar,
            newstar->num_thermal_pulses,
            newstar->angular_momentum);
     
-    if(newstar->age < -TINY)
+    if(newstar->age < -NEGATIVE_AGE_CHECK_MAGNITUDE)
     {
         printf("warning : Negative age (=%g) at stellar_structure entrance for star %d (M=%g M0=%g Mc=%g stellar_type=%d)\n",
                newstar->age,
@@ -549,7 +550,7 @@ void stellar_structure_BSE(struct star_t *newstar,
                           newstar->core_mass,
                           newstar->luminosity);     
         }
-        else if(newstar->age < 0.0)
+        else if(newstar->age < -NEGATIVE_AGE_CHECK_MAGNITUDE)
         {
             Exit_binary_c(BINARY_C_OUT_OF_RANGE,
                           "Negative age = %g at stellar_structure exit!\n",
diff --git a/src/stellar_structure/stellar_structure_HG.c b/src/stellar_structure/stellar_structure_HG.c
index 9dac236c62eee64f19953c52b758073380022e3b..70cfa2adeb7a3eb4bd06dccaf7425588c830831e 100644
--- a/src/stellar_structure/stellar_structure_HG.c
+++ b/src/stellar_structure/stellar_structure_HG.c
@@ -19,10 +19,6 @@ Stellar_type stellar_structure_HG(struct star_t * RESTRICT const newstar,
   
     Dprint("Hertzsprung Gap Star\n");
 
-#ifdef NELEMANS_MASS_PLOT
-    Boolean gijs_out=FALSE;
-    if(Is_zero(newstar->core_mass) && newstar->starnum==0) gijs_out=TRUE;
-#endif
 
 #ifdef MAIN_SEQUENCE_STRIP
     /*
@@ -40,14 +36,6 @@ Stellar_type stellar_structure_HG(struct star_t * RESTRICT const newstar,
     newstar->core_mass = mcbgbf * ((1.0-tau)*rho+tau);
     newstar->core_mass = Max(newstar->core_mass,mcx);
     
-#ifdef NELEMANS_MASS_PLOT
-    if(gijs_out==TRUE)
-    {
-        printf("GIJS1 %12.12e\n",mt);
-        printf("GIJS2 %12.12e\n",newstar->core_mass);
-    }
-#endif
-
 #ifdef NUCSYN
     if(!RLOF_CALL)
     {