From 09734e4bade69cd938913018fd587c7092ee9106 Mon Sep 17 00:00:00 2001
From: Robert Izzard <r.izzard@surrey.ac.uk>
Date: Mon, 28 Sep 2015 18:18:13 +0000
Subject: [PATCH] update with changes to match amanda's models

---
 amanda_sn_plot.pl              | 93 ++++++++++++++++++++++++++--------
 src/binary_structures.h        |  4 ++
 src/hrdiag/hrdiag_WD.c         |  3 --
 src/logging/open_log_files.c   |  3 +-
 src/nucsyn/nucsyn_mc1tp.c      | 20 ++++++++
 src/nucsyn/nucsyn_parameters.h |  2 +-
 src/setup/parse_arguments.c    |  5 ++
 src/zfuncs/dmcHe.c             | 32 ++++++++++++
 src/zfuncs/mc1tp.c             | 18 -------
 tbse                           |  2 +-
 10 files changed, 135 insertions(+), 47 deletions(-)
 create mode 100644 src/zfuncs/dmcHe.c

diff --git a/amanda_sn_plot.pl b/amanda_sn_plot.pl
index 66f85469e..83e59098d 100755
--- a/amanda_sn_plot.pl
+++ b/amanda_sn_plot.pl
@@ -13,43 +13,92 @@ use rob_misc;
 # mcx_eagbf
 # mctmsf (in hrdiag_HG?)
 
+# default core limits:
+
+# minimum_mass_for_carbon_ignition (ONe core) => 1.6
+# minimum_mass_for_neon_ignition (hence SNII) => 2.85
+
+my $opts = {
+#    0.28 => {
+ #       minimum_mass_for_carbon_ignition => 1.6,
+  #      minimum_mass_for_neon_ignition => 2.85
+   # },
+};
+
+my @masslist = (4,5,5.5,6.0,6.25,6.5,6.75,7.0,7.5,7.75,8.0,9.7,10,11,12);
+    
+
 foreach my $Y (0.28, 0.35, 0.40)
 #my $Y = 0.28;
 {
-    for(my $M=5.0; $M<=10; $M+=0.25)
+    if($Y==0.28)
     {
-        print "$Y $M ";
-        my $log = `tbse --init_abund 5 $Y --M_1 $M --log_filename /dev/stdout`;
-        my @x = grep {/^SUPERNOVA/} split(/\n/,$log);
-        
-        if($x[0] && $x[0]=~/in a type (\d+) \((\S+)\)/)
+        foreach my $M (@masslist)
         {
-            # supernova
-            print "$1 $2\n";
+            printf "% 10.2f",$M;
         }
-        else
+        print "\n";
+    }
+    printf "  % 5s\n ",'Y='.sprintf('%.2f',$Y);
+
+
+ #   foreach my $dM (-1.0,-0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0.0,0.1,0.2)
+    {
+#        printf "dM=% 5.2f ",$dM;
+
+        #for(my $M=5.0; $M<=12; $M+=$dM)
+        foreach my $M (@masslist)
         {
-            # WD : find out which type
-            if($log =~/ 13700.0000\s+\S+\s+\S+\s+(\d+)\s+0\s+\S+/)
+            my $args;
+            
+            #$args = "--amanda_dMHe $dM ";
+            $args .= "--init_abund 5 $Y --M_1 $M --log_filename /dev/stdout ";
+            if(defined $opts->{$Y})
             {
-                if($1 == 11)
+                map
                 {
-                    print "-1 COWD\n";
-                }
-                elsif($1 == 12)
+                    $args .= '--'.$_.' '.$opts->{$Y}->{$_}.' ';
+                }sort keys %{$opts->{$Y}};
+            }
+
+            #print "CALL $args\n";
+            my $log = `tbse $args`;
+            my @x = grep {/^SUPERNOVA/} split(/\n/,$log);
+            
+            my $type;
+            my $string;
+            if($x[0] && $x[0]=~/in a type (\d+) \((\S+)\)/)
+            {
+                # supernova
+                ($type,$string) = ($1,$2);
+            }
+            else
+            {
+                # WD : find out which type
+                if($log =~/ 13700.0000\s+\S+\s+\S+\s+(\d+)\s+0\s+\S+/)
                 {
-                    print "-2 ONeWD\n";
+                    if($1 == 11)
+                    {
+                        ($type,$string) = (-1,'COWD');
+                    }
+                    elsif($1 == 12)
+                    {
+                        ($type,$string) = (-2,'ONeWD');
+                    }
+                    else
+                    {
+                        ($type,$string) = (0,'FAIL');
+                    }
                 }
                 else
                 {
-                    print "WD TYPE FAILED\n";
+                    ($type,$string) = (0,'SCAN');
+                    #print "LOG SCAN FAILED : \n$log\n";
                 }
             }
-            else
-            {
-                print "LOG SCAN FAILED : \n$log\n";
-            }
+            
+            printf "% 9s ",$string;
         }
+        print "\n";
     }
-    print "\n";
 }
diff --git a/src/binary_structures.h b/src/binary_structures.h
index eb549aad3..9321c1998 100644
--- a/src/binary_structures.h
+++ b/src/binary_structures.h
@@ -447,6 +447,10 @@ struct preferences_t {
     Boolean show_minimum_separation_for_instant_RLOF;
     Boolean show_minimum_orbital_period_for_instant_RLOF;
 
+#ifdef AMANDA_HELIUM_CORE_SHIFT
+    double amanda_DMHe;
+#endif
+
 };
 
 /************************************************************/
diff --git a/src/hrdiag/hrdiag_WD.c b/src/hrdiag/hrdiag_WD.c
index 5be5c2325..783b3270e 100644
--- a/src/hrdiag/hrdiag_WD.c
+++ b/src/hrdiag/hrdiag_WD.c
@@ -42,8 +42,6 @@ void hrdiag_WD(double *aj,
              * Note that the RLOF cases are dealt with 
              * in RLOF/SNIa_check
              */
-        
-
 #ifdef LOG_SUPERNOVAE
             log_sn(stardata,star->starnum,SN_IA_CHAND,PRE_SN);
 #endif  
@@ -74,7 +72,6 @@ void hrdiag_WD(double *aj,
             *r = (*mt < 5e-6) ? 0.009 : 0.09;
         }
         Dprint("r=%.12g\n",*r);
-
     }
 
 #ifdef LOG_ACCRETED_MASS
diff --git a/src/logging/open_log_files.c b/src/logging/open_log_files.c
index 71036606c..a259c29ee 100644
--- a/src/logging/open_log_files.c
+++ b/src/logging/open_log_files.c
@@ -1,6 +1,5 @@
 #include "../binary_c.h"
 
-
 /*
  * Function to open binary_c's log files
  */
@@ -10,7 +9,7 @@ void open_log_files(FILE *(*log_fp),
 #ifdef FILE_LOG
   
 #ifdef FILE_LOG_REWIND
-    if(*log_fp!=NULL)
+    if(*log_fp!=NULL && *log_fp!=stdout && *log_fp!=stderr)
     {
         /* log is already open, just rewind and truncate it */
         Dprint("rewind log_fp\n");
diff --git a/src/nucsyn/nucsyn_mc1tp.c b/src/nucsyn/nucsyn_mc1tp.c
index d0a0f9cdb..eb72b6623 100644
--- a/src/nucsyn/nucsyn_mc1tp.c
+++ b/src/nucsyn/nucsyn_mc1tp.c
@@ -228,6 +228,26 @@ double nucsyn_mc1tp (double m,
     y=(1.0-fl)*y + fl * mc1tpj;
     Dprint("%g\n",y);
 
+
+    
+#ifdef AMANDA_HELIUM_CORE_SHIFT
+        {
+            /* 
+             * Approximate shift in total stellar mass that makes an ONeWD
+             * according to Amanda's helium-rich, Z=0.014 models:
+             * Y=0.28 dm = -0.6
+             * Y=0.35 dm = -0.2
+             * Y=0.38 dm = 0.0
+             * Fit : dm = (-2.26960e+00)+(5.94940e+00)*Y
+             *
+             * see http://bulgepositrons.pbworks.com/w/file/fetch/100788601/grid.pdf
+             */
+           
+           y += dmcHe(m,star,stardata);
+        }
+#endif
+
+
 #ifdef NANCHECKS
     if(isnan(y))
     {
diff --git a/src/nucsyn/nucsyn_parameters.h b/src/nucsyn/nucsyn_parameters.h
index d92232582..b42cfdf47 100644
--- a/src/nucsyn/nucsyn_parameters.h
+++ b/src/nucsyn/nucsyn_parameters.h
@@ -128,7 +128,7 @@
  * mc1tp, mcmin, lambda etc.
  */
 
-//#define AMANDA_HELIUM_CORE_SHIFT
+#define AMANDA_HELIUM_CORE_SHIFT
 
 /* define this to use fits to the Padova group Mc1tp over the mass
  * range M<3 Msun rather than Amanda's fits */
diff --git a/src/setup/parse_arguments.c b/src/setup/parse_arguments.c
index 24e6f7b9f..0f5bd0dea 100644
--- a/src/setup/parse_arguments.c
+++ b/src/setup/parse_arguments.c
@@ -1268,6 +1268,11 @@ void parse_arguments(int start,
 #ifdef BATCHMODE
                 stardata->preferences->no_arg_checking=FALSE;
 #endif
+            } 
+            EARGIS("--amanda_dMHe")
+            {
+                stardata->preferences->amanda_DMHe = atof(next_argument);
+                ++i;
             }
             BREAK_OR_ERROR;
          
diff --git a/src/zfuncs/dmcHe.c b/src/zfuncs/dmcHe.c
new file mode 100644
index 000000000..d490f7c77
--- /dev/null
+++ b/src/zfuncs/dmcHe.c
@@ -0,0 +1,32 @@
+#include "../binary_c.h"
+
+#ifdef AMANDA_HELIUM_CORE_SHIFT
+
+double dmcHe(const double m,
+             const struct star_t * star,
+             const struct stardata_t * RESTRICT stardata)
+{
+    double Y = stardata->preferences->the_initial_abundances[XHe4];
+    
+    // first, from cmd line
+    double dm = stardata->preferences->amanda_DMHe;
+
+    /*
+     * Y=0.28 -> -0.7
+     * Y=0.35 -> -0.2
+     * Y=0.40 -> -0.1
+     */
+
+    /* exact fit to the above, modulated by turn-on function */
+    dm  = m>5.0 ? 
+        ( 1.0/(1.0+pow(0.01,m-5.0)) * (-6.90 + 3.41430e+01*Y + -4.28570e+01*Y*Y) ) 
+        : 0.0;
+
+    /* limit to prevent problems */
+    dm = MAX(dm,-0.5*star->core_mass);
+    
+    //fprintf(stderr,"Y=%g M=%g Mc=%g -> DM=%g\n",Y,m,star->core_mass,dm);
+    return (dm);
+}
+
+#endif// AMANDA_HELIUM_CORE_SHIFT
diff --git a/src/zfuncs/mc1tp.c b/src/zfuncs/mc1tp.c
index d8987e17f..e6ebee3c0 100644
--- a/src/zfuncs/mc1tp.c
+++ b/src/zfuncs/mc1tp.c
@@ -92,24 +92,6 @@ double guess_mc1tp(double mass,
     Dprint("no nucsyn use mcagbf %g\n",mc1);
 #endif
 
-#ifdef AMANDA_HELIUM_CORE_SHIFT
-        {
-            /* 
-             * Approximate shift in total stellar mass that makes an ONeWD
-             * according to Amanda's helium-rich, Z=0.014 models:
-             * Y=0.28 dm = -0.6
-             * Y=0.35 dm = -0.2
-             * Y=0.38 dm = 0.0
-             * Fit : dm = (-2.26960e+00)+(5.94940e+00)*Y
-             *
-             * see http://bulgepositrons.pbworks.com/w/file/fetch/100788601/grid.pdf
-             */
-           
-            mc1 += dmcHe(mass,star,stardata);
-        }
-#endif
-
-
     /* the initial CO core mass must be > the "compact core mass" */
     mc1=MAX(star->compact_core_mass,mc1);
 
diff --git a/tbse b/tbse
index 33656ae25..80c8e7a9a 100755
--- a/tbse
+++ b/tbse
@@ -409,7 +409,7 @@ MAX_TPAGB_CORE_MASS=1.38
 CHANDRASEKHAR_MASS=1.44
 MAX_NEUTRON_STAR_MASS=1.8
 MINIMUM_MASS_FOR_CARBON_IGNITION=1.6
-MINIMUM_MASS_FOR_NEON_IGNITION=1.6
+MINIMUM_MASS_FOR_NEON_IGNITION=2.85
 
 # include nuclear reaction rate multipliers in the form
 # --nucreacmult a f
-- 
GitLab