fix issues with start of evolution in the FORTRAN API: we now only

call the initialization function if stardata->model.model_number == 0

updates for apitest to use the latest code

fixed Meson's handling of libdirs so binary_c-config works

cleaned up some unused logging options

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"

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
 }
 

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,&

apitest/apitest.sh

 #!/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

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)

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 */

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);

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
 

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
 

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;
 

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");
+    
 }
 

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! */

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;

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))
     {

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);

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",

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)
     {