From 5a8cf77d0e40f61f356c3a82494a59ae8e24907a Mon Sep 17 00:00:00 2001
From: Robert Izzard <r.izzard@surrey.ac.uk>
Date: Tue, 4 Jun 2019 13:36:10 +0100
Subject: [PATCH] fix mass ratio used in Claeys factor (thanks to Karel for
 spotting this!)

---
 src/RLOF/RLOF_mass_transfer_rate.c             |  8 ++++----
 src/binary_c_debug.h                           | 11 +++++------
 src/binary_star_functions/tides.c              |  2 +-
 .../common_envelope_evolution_BSE.c            |  4 ++--
 src/evolution/evolve_system_binary_c.c         |  3 +--
 src/evolution/stellar_evolution.c              | 10 ++--------
 src/logging/log_every_timestep.c               |  4 ++++
 tbse                                           |  2 +-
 test_new.sh                                    | 18 ++++++++++++------
 valgrind.pl                                    |  2 +-
 10 files changed, 33 insertions(+), 31 deletions(-)

diff --git a/src/RLOF/RLOF_mass_transfer_rate.c b/src/RLOF/RLOF_mass_transfer_rate.c
index eee78d1eb..a075a2ec9 100644
--- a/src/RLOF/RLOF_mass_transfer_rate.c
+++ b/src/RLOF/RLOF_mass_transfer_rate.c
@@ -549,8 +549,8 @@ double hurley_rate(double r,
 static double Claeys_factor(struct star_t * loser)
 {
     double f;
- 
-    if(loser->q < 1.0+TINY)
+    const double Q = 1.0 / loser->q;
+    if(Q < 1.0+TINY)
     {
         f = 1000.0;
     }
@@ -562,8 +562,8 @@ static double Claeys_factor(struct star_t * loser)
          * very rare.
          */
         f = MAX(1.0,
-                1000.0/loser->q*
-                exp(-0.5*POW2(log(MIN(100.0,loser->q))/0.15))
+                1000.0/Q *
+                exp(-0.5*POW2(log(MIN(100.0,Q))/0.15))
             );
     }
     return f;
diff --git a/src/binary_c_debug.h b/src/binary_c_debug.h
index b48c75baa..48e529d3f 100644
--- a/src/binary_c_debug.h
+++ b/src/binary_c_debug.h
@@ -20,7 +20,7 @@
  *
  * See below for further options.
  */
-#define DEBUG 0
+#define DEBUG 1
 
 /*
  * To remove debugging possibility from the code, globally, hence
@@ -57,7 +57,7 @@
  * variable.
  */
 //#define Debug_expression (1)
-#define Debug_expression (stardata!=NULL && stardata->model.time>1419)
+#define Debug_expression (stardata!=NULL && stardata->model.time>9.0)
 
 /*
  * If you define Debug_stop_expression, and it is at any time TRUE,
@@ -74,7 +74,7 @@
 
 #define Debug_stop_expression \
     (stardata->star[0].stellar_type==7)
-//#undef Debug_stop_expression 
+#undef Debug_stop_expression 
 
 /*
  * Debug_show_expression, if defined, is shown at the beginning
@@ -86,9 +86,8 @@
  * this is not shown.
  */
 
-#define Debug_show_expression " st 0=%d 1=%d ",            \
-        stardata->star[0].stellar_type,                    \
-        stardata->star[1].stellar_type
+#define Debug_show_expression " sep=%g ",       \
+        stardata->common.orbit.separation
 
 //#undef Debug_show_expression
 
diff --git a/src/binary_star_functions/tides.c b/src/binary_star_functions/tides.c
index e17b457ab..c2c5ede53 100644
--- a/src/binary_star_functions/tides.c
+++ b/src/binary_star_functions/tides.c
@@ -264,7 +264,7 @@ void tides(struct stardata_t *stardata,
              */
             const double omega_estimate = star->omega + 
                 star->derivative[DERIVATIVE_STELLAR_ANGULAR_VELOCITY_TIDES] * stardata->model.dt;
-            if(star->omega * omega_estimate < 0.0)
+            if(0 && star->omega * omega_estimate < 0.0)
             {
                 /*
                  * Instead, assume omega = 0 on the next step
diff --git a/src/common_envelope/common_envelope_evolution_BSE.c b/src/common_envelope/common_envelope_evolution_BSE.c
index 44275a578..a291c0b65 100644
--- a/src/common_envelope/common_envelope_evolution_BSE.c
+++ b/src/common_envelope/common_envelope_evolution_BSE.c
@@ -597,8 +597,8 @@ int common_envelope_evolution_BSE (
             final_separation = donor->core_mass * accretor->mass/(2.0*EorbF);
         }
 
-#warning REMOVE ME
-        final_separation =4.4468;
+//#warning REMOVE ME
+        //final_separation =4.4468;
 
         Dprint("Final separation = %g (mc1=%g m2=%g EorbF=%g)\n",
                final_separation,
diff --git a/src/evolution/evolve_system_binary_c.c b/src/evolution/evolve_system_binary_c.c
index 5442be1a3..896e4bba8 100644
--- a/src/evolution/evolve_system_binary_c.c
+++ b/src/evolution/evolve_system_binary_c.c
@@ -115,8 +115,7 @@ int evolve_system_binary_c(struct stardata_t * RESTRICT stardata)
                            stardata->star[1].radius / stardata->star[1].roche_radius
                         );
 
-                    if(1 &&
-//stardata->model.in_RLOF == TRUE &&
+                    if(//stardata->model.in_RLOF == TRUE &&
                        test_if_primary_still_fills_roche_lobe(stardata)==CONTACT)
                     {
                         printf("EVRLOF -> contact\n");
diff --git a/src/evolution/stellar_evolution.c b/src/evolution/stellar_evolution.c
index f31a97ccc..9a76db79e 100644
--- a/src/evolution/stellar_evolution.c
+++ b/src/evolution/stellar_evolution.c
@@ -23,8 +23,8 @@ static void set_derivatives(struct star_t * old,
 int stellar_evolution(struct stardata_t * RESTRICT stardata,
                       const int codetype)
 {
-    struct common_t *common=&(stardata->common);
-    struct model_t *model=&(stardata->model);
+    struct common_t *common MAYBE_UNUSED = &stardata->common;
+    struct model_t *model = &stardata->model;
     double age,m0,mc_1tp,mt,mc,lum,rm,rc,mcCO,mcGB,mcmaxMS;
     age=m0=mc_1tp=mt=mc=lum=rm=rc=mcCO=mcGB=mcmaxMS=0.0;
     Star_number k;
@@ -362,12 +362,6 @@ int stellar_evolution(struct stardata_t * RESTRICT stardata,
                dt,
                model->time_remaining,
                MINIMUM_STELLAR_TIMESTEP);
-
-        if(0&&stardata->model.time > 3812.2)
-        {
-            fflush(NULL);
-            _exit(0);
-        }
         
         /*
          * Update stellar structure
diff --git a/src/logging/log_every_timestep.c b/src/logging/log_every_timestep.c
index 1599208c9..abb1e224d 100644
--- a/src/logging/log_every_timestep.c
+++ b/src/logging/log_every_timestep.c
@@ -2471,5 +2471,9 @@ void log_every_timestep(struct stardata_t * RESTRICT stardata)
                stardata->star[0].time_next_pulse
             );
     }
+
+    printf("At %g sep = %g\n",
+           stardata->model.time,
+           stardata->common.orbit.separation);
     
 }
diff --git a/tbse b/tbse
index e83dc7077..51a70ee8e 100755
--- a/tbse
+++ b/tbse
@@ -65,7 +65,7 @@ RANDOM_SYSTEMS=0
 # 1 = RK2
 # 2 = RK4
 # 3 = Predictor-corrector
-SOLVER=2
+SOLVER=0
 
 # use new events framework?
 USE_EVENTS=1
diff --git a/test_new.sh b/test_new.sh
index 388c44de6..e07284c62 100755
--- a/test_new.sh
+++ b/test_new.sh
@@ -1,16 +1,22 @@
 #!/bin/bash
 
-INARGS="$@_"
+INARGS="$@"
+RANDOM_ARGS=`tbse random_systems 1  max_evolution_time 0|grep At|head -1|perl -lane 'print "@F[4..$#F]"'`
+ARGS="$RANDOM_ARGS maximum_timestep 0.1 solver 0 $INARGS"
 
-ARGS="maximum_timestep 0.1 use_events 1"
+echo
+echo "Args random : $RANDOM_ARGS"
+echo
+echo "Args first : $ARGS"
+echo
+echo "Args in : $INARGS"
+echo
 
 echo "Run 0"
-tbse $INARGS $ARGS solver 0 evolution_algorithm 0 log_filename /tmp/c_log2.0 use_events 0 > 0.out && echo "0 done" &
+tbse $INARGS $ARGS evolution_algorithm 0 log_filename /tmp/c_log2.0 use_events 0 > 0.out && echo "0 done" &
 
 echo "Run 1"
-tbse $INARGS $ARGS solver 0 evolution_algorithm 1 log_filename /tmp/c_log2.1 use_events 1 | filter MMM > mmm.dat
-
-#> 1.out && echo "1 done" &
+tbse $INARGS $ARGS evolution_algorithm 1 log_filename /tmp/c_log2.1 use_events 1 > 1.out && echo "1 done" &
 
 wait
 
diff --git a/valgrind.pl b/valgrind.pl
index b7b84b663..c6ec87e1d 100755
--- a/valgrind.pl
+++ b/valgrind.pl
@@ -26,7 +26,7 @@ my $settle = 10; # wait this many seconds before output
 my $sleeptime = 0.25; # number of s to wait between checking for binary_c output 
 my $twarn = 180.0; # warn if no output for this time
 my $warn_every = 10.0; # warn every this many seconds
-my $n_binary_c_threads = 8;#rob_misc::ncpus()-1;    
+my $n_binary_c_threads = 12;#rob_misc::ncpus()-1;    
 my $nstore = 1000;# store this many lines of output
 my $q = robqueue->new(
     nthreads=>$n_binary_c_threads+1,
-- 
GitLab