From 3a7650bee5184c06f99bb973c7ca011aa5fdcf7a Mon Sep 17 00:00:00 2001
From: Robert Izzard <r.izzard@surrey.ac.uk>
Date: Fri, 29 Sep 2017 09:58:50 +0000
Subject: [PATCH] code is behaving oddly when counting CEMP stars and changing
 the lower mass limit from 0.1 to 0.08

commit so I can work on it later today
---
 cemps                          | 169 +++++++++++++++++++++++++++++++++
 src/nucsyn/nucsyn_cemp_log.c   |  70 +++-----------
 src/perl/scripts2/cemps2017.pl |  81 ++++++++--------
 tbse                           |   2 -
 4 files changed, 223 insertions(+), 99 deletions(-)
 create mode 100644 cemps

diff --git a/cemps b/cemps
new file mode 100644
index 000000000..35c9f98e5
--- /dev/null
+++ b/cemps
@@ -0,0 +1,169 @@
+;; This buffer is for notes you don't want to save, and for Lisp evaluation.
+;; If you want to create a file, visit that file with C-x C-f,
+;; then enter the text in that file's own buffer.
+
+40-cubed grid with mmin = 0.1?
+
+standard 
+      CEMP  1.15 % 0.223161
+     CNEMP  0.02 % 0.00389339
+      XEMP 98.83 % 19.2273
+
+flat-lna
+      CEMP  1.53 % 0.25037
+     CNEMP  0.02 % 0.00261962
+      XEMP 98.45 % 16.0869
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+40-cubed grid with mmin = 0.1
+----------
+Redshift 0, age 0
+
+      CEMP  1.15 % 0.226738
+     CNEMP  0.02 % 0.0039558
+      XEMP 98.83 % 19.5355
+----------
+Redshift 0.073271, age 1
+
+      CEMP  1.15 % 0.226738
+     CNEMP  0.02 % 0.0039558
+      XEMP 98.83 % 19.5355
+----------
+Redshift 1.70363, age 10
+
+      CEMP  1.15 % 0.226738
+     CNEMP  0.02 % 0.0039558
+      XEMP 98.83 % 19.5355
+----------
+Redshift 2.31993, age 11
+
+      CEMP  1.15 % 0.226738
+     CNEMP  0.02 % 0.0039558
+      XEMP 98.83 % 19.5355
+----------
+Redshift 3.43854, age 12
+
+      CEMP  1.15 % 0.226738
+     CNEMP  0.02 % 0.0039558
+      XEMP 98.83 % 19.5355
+----------
+Redshift 6.42945, age 13
+
+      CEMP  1.15 % 0.232298
+     CNEMP  0.02 % 0.00405279
+      XEMP 98.83 % 20.0118
+----------
+Redshift 8.87862, age 13.3
+
+      CEMP  1.30 % 0.396665
+     CNEMP  0.02 % 0.00692073
+      XEMP 98.67 % 30.0433
+----------
+Redshift 10.2538, age 13.4
+
+      CEMP  1.66 % 0.531308
+     CNEMP  0.03 % 0.00928702
+      XEMP 98.31 % 31.5082
+----------
+Redshift 12.2342, age 13.5
+
+      CEMP  2.08 % 0.712213
+     CNEMP  0.04 % 0.0136675
+      XEMP 97.88 % 33.5596
+----------
+Redshift 15.4011, age 13.6
+
+      CEMP  2.24 % 0.843418
+     CNEMP  0.06 % 0.0232776
+      XEMP 97.69 % 36.7075
+----------
+Redshift 21.5464, age 13.7
+
+      CEMP  2.38 % 0.972373
+     CNEMP  0.07 % 0.0298953
+      XEMP 97.55 % 39.8494
+----------
+Redshift 41.3612, age 13.8
+
+      CEMP  3.95 % 0.796859
+     CNEMP  0.20 % 0.0408007
+      XEMP 95.85 % 19.3552
+----------
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+40-cubed grid with mmin=0.08
+
+Redshift 0, age 0
+
+      CEMP  0.71 % 0.137588
+     CNEMP  0.02 % 0.00439895
+      XEMP 99.27 % 19.2267
+----------
+Redshift 0.073271, age 1
+
+      CEMP  0.71 % 0.137588
+     CNEMP  0.02 % 0.00439895
+      XEMP 99.27 % 19.2267
+----------
+Redshift 1.70363, age 10
+
+      CEMP  0.71 % 0.137588
+     CNEMP  0.02 % 0.00439895
+      XEMP 99.27 % 19.2267
+----------
+Redshift 2.31993, age 11
+
+      CEMP  0.71 % 0.137588
+     CNEMP  0.02 % 0.00439895
+      XEMP 99.27 % 19.2267
+----------
+Redshift 3.43854, age 12
+
+      CEMP  0.71 % 0.137588
+     CNEMP  0.02 % 0.00439895
+      XEMP 99.27 % 19.2267
+----------
+Redshift 6.42945, age 13
+
+      CEMP  0.71 % 0.141989
+     CNEMP  0.02 % 0.00453966
+      XEMP 99.27 % 19.8398
+----------
+Redshift 8.87862, age 13.3
+
+      CEMP  0.81 % 0.259519
+     CNEMP  0.03 % 0.00829788
+      XEMP 99.16 % 31.7348
+----------
+Redshift 10.2538, age 13.4
+
+      CEMP  1.04 % 0.353635
+     CNEMP  0.03 % 0.0113082
+      XEMP 98.93 % 33.5882
+----------
+Redshift 12.2342, age 13.5
+
+      CEMP  1.35 % 0.491915
+     CNEMP  0.05 % 0.0168624
+      XEMP 98.61 % 36.0149
+----------
+Redshift 15.4011, age 13.6
+
+      CEMP  1.48 % 0.598692
+     CNEMP  0.07 % 0.0290474
+      XEMP 98.45 % 39.8051
+----------
+Redshift 21.5464, age 13.7
+
+      CEMP  1.58 % 0.697132
+     CNEMP  0.08 % 0.0345455
+      XEMP 98.34 % 43.4588
+----------
+Redshift 41.3612, age 13.8
+
+      CEMP  2.78 % 0.596026
+     CNEMP  0.21 % 0.0440141
+      XEMP 97.01 % 20.7733
+----------
diff --git a/src/nucsyn/nucsyn_cemp_log.c b/src/nucsyn/nucsyn_cemp_log.c
index c60068469..bf77b45e4 100644
--- a/src/nucsyn/nucsyn_cemp_log.c
+++ b/src/nucsyn/nucsyn_cemp_log.c
@@ -38,52 +38,16 @@ void nucsyn_cemp_log(struct stardata_t * RESTRICT stardata)
     unsigned int k;
     char s[100];
   
-/*  { */
-/*     star=&(stardata->star[0]); */
-/*     if(star->stellar_type == 6) */
-/*       { */
-	
-/* 	Abundance * Xenv=nucsyn_observed_surface_abundances(star); */
-/* 	printf("AGB %g %g\n", */
-/* 	       stardata->model.time, */
-/* 	       nucsyn_square_multibracket(Xenv,Xsolar, */
-/* 					  (Isotope[]){XC12,XC13},2, */
-/* 					  (Isotope[]){XFe52,XFe54,XFe56,XFe57,XFe58},5)); */
-/*       } */
-/*   } */
-
-    //stardata->preferences->CEMP_minimum_age=0.0;
-
-    if(stardata->model.time<TINY)stardata->model.cemp_log_prevt=0.0; // first timestep
-    dt=stardata->model.time-MAX(stardata->model.cemp_log_prevt,
-                                stardata->preferences->CEMP_minimum_age)
-
+    if(stardata->model.time<TINY) 
+    {
+        stardata->model.cemp_log_prevt=0.0; // first timestep
+    }
+    dt = stardata->model.time - MAX(stardata->model.cemp_log_prevt,
+                                    stardata->preferences->CEMP_minimum_age)
         ; // calculate time since previous output
-    //printf("dt=%g from t=%g prevt=%g\n",dt,stardata->model.time,stardata->model.cemp_log_prevt);
-    //if(dt<CEMP_DT) return;
-    stardata->model.cemp_log_prevt=stardata->model.time;
-  
-    /*
-      printf("PAPERLOG %g %g %g %g %g %g %g %g\n",
-      stardata->model.time,
-      stardata->common.orbital_period,
-      stardata->star[0].mass,
-      stardata->star[1].mass,
-      nucsyn_square_multibracket(stardata->star[0].Xenv,Xsolar,
-      (Isotope[]){XC12,XC13},2,
-      (Isotope[]){XFe52,XFe54,XFe56,XFe57,XFe58},5),
-      nucsyn_square_multibracket(stardata->star[1].Xenv,Xsolar,
-      (Isotope[]){XC12,XC13},2,
-      (Isotope[]){XFe52,XFe54,XFe56,XFe57,XFe58},5),	 
-      nucsyn_square_multibracket(stardata->star[0].Xenv,Xsolar,
-      (Isotope[]){XN14},1,
-      (Isotope[]){XFe52,XFe54,XFe56,XFe57,XFe58},5),
-      nucsyn_square_multibracket(stardata->star[1].Xenv,Xsolar,
-      (Isotope[]){XN14},1,
-      (Isotope[]){XFe52,XFe54,XFe56,XFe57,XFe58},5));
-
-    */
 
+    stardata->model.cemp_log_prevt = stardata->model.time;
+  
     /*
      * Look for CEMP stars
      */
@@ -183,9 +147,9 @@ void nucsyn_cemp_log(struct stardata_t * RESTRICT stardata)
 	   * log_10 (g) < 4.0 - same as the observations
 	   * and pre-WD
 	   */
-	  if((star->mass > 0.1)&&
-	     (log10g<stardata->preferences->CEMP_logg_maximum)&&
-	     (star->stellar_type<HeWD))
+	  if(star->stellar_type < HeWD &&
+             star->mass > 0.1 &&
+	     log10g < stardata->preferences->CEMP_logg_maximum)
           {
               // EMP giant
 
@@ -198,18 +162,6 @@ void nucsyn_cemp_log(struct stardata_t * RESTRICT stardata)
 	       * (usually 1.0)
 	       */
 	      
-	      
-	   
-	      
-	      /*
-                for(cfe_bracket=0.01;cfe_bracket<=3.0;cfe_bracket+=0.1)
-		{
-                for(nfe_bracket=0.01;nfe_bracket<=3.0;nfe_bracket+=0.1)
-                {
-                cn_bracket=cfe_bracket - nfe_bracket;
-
-                PRINTF("%g %g ",CEMP_CRITERIA,cfe_bracket,nfe_bracket);
-	      */
 
 	      if(CEMP_CRITERIA)
               {
diff --git a/src/perl/scripts2/cemps2017.pl b/src/perl/scripts2/cemps2017.pl
index 86663d8b5..ebe07f5a4 100755
--- a/src/perl/scripts2/cemps2017.pl
+++ b/src/perl/scripts2/cemps2017.pl
@@ -93,6 +93,11 @@ $population->set(
     AGB_luminosity_algorithm=>2,
     AGB_3dup_algorithm=>2,
 
+    # CEMP options
+    CEMP_minimum_age => 10e3, # 10 Gyr 
+    CEMP_logg_maximum => 4.0, # logg < 4.0 for CEMP stars 
+    CEMP_cfe_minimum => 1.0, # [C/Fe] > 1 for CEMP stars
+
     );
 # scan command line arguments for extra options
 $population->parse_args();  
@@ -100,14 +105,14 @@ $population->parse_args();
 # duplicity is 0 for single stars, 1 for binary stars
 # and 2 for a mixed population sampled at equal times
 my $duplicity = 1;
-my $mmin = 0.1;
+my $mmin = 0.08;
 
 if($duplicity == 0)
 {
     # make a grid of $nstars single binary stars, log-spaced,
     # with masses between $mmin and $mmax
     my $nstars = 40;
-   
+    
     my $mmax = 80.0;
     $population->add_grid_variable(
         'name'       => 'lnm1', 
@@ -116,7 +121,7 @@ if($duplicity == 0)
         'resolution' => $nstars, # just a counter for the grid
         'spacingfunc'=>"const(log($mmin),log($mmax),$nstars)",
         'precode'    =>'$m1=exp($lnm1);',
-        'probdist'   =>'Kroupa2001($m1)*$m1',
+        'probdist'   =>'Kroupa2001($m1,{m0=>'.$mmin.'})*$m1',
         'dphasevol'  =>'$dlnm1',
         );
 }
@@ -135,30 +140,30 @@ elsif($duplicity == 1)
     
     $population->add_grid_variable
         (
-        'name'       => 'lnm1', 
-        'longname'   =>'Primary mass', 
-        'range'      =>[log($mmin),log($mmax)],
-        'resolution' => $resolution->{m1},
-        'spacingfunc'=>"const(log($mmin),log($mmax),$resolution->{m1})",
-        'precode'    =>'$m1=exp($lnm1);',
-        'probdist'   =>'Kroupa2001($m1)*$m1',
-        'dphasevol'  =>'$dlnm1',
+         'name'       => 'lnm1', 
+         'longname'   =>'Primary mass', 
+         'range'      =>[log($mmin),log($mmax)],
+         'resolution' => $resolution->{m1},
+         'spacingfunc'=>"const(log($mmin),log($mmax),$resolution->{m1})",
+         'precode'    =>'$m1=exp($lnm1);',
+         'probdist'   =>'Kroupa2001($m1,{m0=>'.$mmin.'})*$m1',
+         'dphasevol'  =>'$dlnm1',
         );
     # q=M1/M2 distribution flat in q between $mmin/M1 and 1.0
     $population->add_grid_variable
         (
-        'condition'  =>'$self->{_grid_options}{binary}==1',
-        'name'       =>'q',
-        'longname'   =>'Mass ratio',
-        'range'      =>[$mmin.'/$m1',1.0],
+         'condition'  =>'$self->{_grid_options}{binary}==1',
+         'name'       =>'q',
+         'longname'   =>'Mass ratio',
+         'range'      =>[$mmin.'/$m1',1.0],
          'resolution'=>$resolution->{q},
-        'spacingfunc'=>"const($mmin/\$m1,1.0,$resolution->{q})",
-        'probdist'   =>"flatsections\(\$q,\[
+         'spacingfunc'=>"const($mmin/\$m1,1.0,$resolution->{q})",
+         'probdist'   =>"flatsections\(\$q,\[
 \{min=>$mmin/\$m1,max=>0.8,height=>1.0\},
 \{min=>0.8,max=>1.0,height=>1.0\},
 \]\)",
-        precode     =>'$m2=$q*$m1;',
-        dphasevol   =>'$dq',
+         precode     =>'$m2=$q*$m1;',
+         dphasevol   =>'$dq',
         );
 
 
@@ -219,24 +224,24 @@ else
     my $time_adaptive_options = {
         max_evolution_time =>
             $population->{_bse_options}{max_evolution_time},
-        stellar_lifetime_table_nm=>1000,
-        time_adaptive_mass_grid_log10_time=>0,
-        time_adaptive_mass_grid_step=>$dt,
-        time_adaptive_mass_grid_nlow_mass_stars=>10,
-        nthreads           =>$nthreads,
-        thread_sleep       =>1,
-        mmin               =>$mmin,
-        mmax               =>80.0,
-        mass_grid_log10_time=>0,
-        mass_grid_step      =>$dt*$sampling_factor,
-        extra_flash_resolution=>0, # 1 = broken?
-        mass_grid_nlow_mass_stars=>10,
-        debugging_output_directory=>'/tmp',
-        max_delta_m         =>1.0,
-        savegrid            =>1,
-        vb                  =>0,
-        metallicity=>$population->{_bse_options}{metallicity},
-        agbzoom             =>0,
+            stellar_lifetime_table_nm=>1000,
+            time_adaptive_mass_grid_log10_time=>0,
+            time_adaptive_mass_grid_step=>$dt,
+            time_adaptive_mass_grid_nlow_mass_stars=>10,
+            nthreads           =>$nthreads,
+            thread_sleep       =>1,
+            mmin               =>$mmin,
+            mmax               =>80.0,
+            mass_grid_log10_time=>0,
+            mass_grid_step      =>$dt*$sampling_factor,
+            extra_flash_resolution=>0, # 1 = broken?
+            mass_grid_nlow_mass_stars=>10,
+            debugging_output_directory=>'/tmp',
+            max_delta_m         =>1.0,
+            savegrid            =>1,
+            vb                  =>0,
+            metallicity=>$population->{_bse_options}{metallicity},
+            agbzoom             =>0,
     };
     distribution_functions::bastard_distribution(
         $population, {
@@ -323,7 +328,7 @@ sub parse_data
     my ($population, $results) = @_;
 
     my $m1init = $population->{_grid_options}->{progenitor_hash}->{m1}; 
-    my $p_Kroupa = distribution_functions::Kroupa2001($m1init);
+    my $p_Kroupa = distribution_functions::Kroupa2001($m1init,{m0=>$mmin});
 
     while(1)
     {
diff --git a/tbse b/tbse
index 38e64d984..12e624ab3 100755
--- a/tbse
+++ b/tbse
@@ -950,8 +950,6 @@ CEMP_CFE_MINIMUM=1.0
 # maximum logg for CEMP stars
 CEMP_LOGG_MAXIMUM=4.0
 
-# minimum stellar type for cemps (2 for giants, 0 for all)
-#CEMP_MINIMUM_STELLAR_TYPE=2
 # minimum age for cemps (12000 Myr)
 CEMP_MINIMUM_AGE=10000
 
-- 
GitLab