From 2a8356f1a25a50b82937ef8d019a170002ff5588 Mon Sep 17 00:00:00 2001
From: David Hendriks <davidhendriks93@gmail.com>
Date: Fri, 26 Mar 2021 11:32:23 +0000
Subject: [PATCH] updated example population

---
 examples/example_population.py | 90 ++++++++++++++++++++++++----------
 1 file changed, 64 insertions(+), 26 deletions(-)

diff --git a/examples/example_population.py b/examples/example_population.py
index 41567d6ba..7f5f776ee 100644
--- a/examples/example_population.py
+++ b/examples/example_population.py
@@ -2,10 +2,8 @@ import os
 
 # import json
 # import time
-# import sys
-
 from binarycpython.utils.grid import Population
-from binarycpython.utils.functions import get_help_all, get_help, create_hdf5
+from binarycpython.utils.functions import get_help_all, get_help, create_hdf5, output_lines
 from binarycpython.utils.custom_logging_functions import temp_dir
 
 #########################################################
@@ -13,15 +11,9 @@ from binarycpython.utils.custom_logging_functions import temp_dir
 # The use of help(<function>) is a good way to inspect what parameters are there to use
 #########################################################
 
-## Quick script to get some output
-def output_lines(output):
-    """
-    Function that outputs the lines that were recieved from the binary_c run.
-    """
-    return output.splitlines()
-
-
 def parse_function(self, output):
+    # EXAMPLE PARSE_FUNCTION
+
     # extract info from the population instance
 
     # Get some information from the
@@ -36,14 +28,20 @@ def parse_function(self, output):
     # Create filename
     outfilename = os.path.join(data_dir, base_filename)
 
+    parameters = ["time", "mass", "zams_mass", "probability", "radius", "stellar_type"]
+
     # Go over the output.
     for el in output_lines(output):
         headerline = el.split()[0]
 
         # CHeck the header and act accordingly
         if headerline == "MY_STELLAR_DATA":
-            parameters = ["time", "mass", "zams_mass", "probability", "radius"]
             values = el.split()[1:]
+            print(values)
+
+            if not len(parameters)==len(values):
+                print("Amount of column names isnt equal to amount of columns")
+                raise ValueError
 
             if not os.path.exists(outfilename):
                 with open(outfilename, "w") as f:
@@ -52,7 +50,6 @@ def parse_function(self, output):
             with open(outfilename, "a") as f:
                 f.write(seperator.join(values) + "\n")
 
-
 # Create population object
 example_pop = Population()
 
@@ -70,9 +67,14 @@ example_pop.set(
     orbital_period=45000000080,  # bse_options
     max_evolution_time=15000,  # bse_options
     eccentricity=0.02,  # bse_options
+
+    # Set companion to low mass
+    M_2=0.08, # Since in the example we run a single system, we should set the companion mass here. If we donm't do this, the code will complain. 
+
     # grid_options
-    amt_cores=1,  # grid_options
+    amt_cores=2,  # grid_options
     verbose=1,  # verbosity. Not fully configured correctly yet but having it value of 1 prints alot of stuff
+
     # Custom options # TODO: need to be set in grid_options probably
     data_dir=os.path.join(
         temp_dir(), "example_python_population_result"
@@ -80,7 +82,6 @@ example_pop.set(
     base_filename="example_pop.dat",  # custom_options
 )
 
-
 # Creating a parsing function
 example_pop.set(
     parse_function=parse_function,  # Setting the parse function thats used in the evolve_population
@@ -97,20 +98,21 @@ example_pop.set(
 # )
 
 
-# Log the moment when the star turns into a hertzsprung-gap
+# Log the moment when the star turns into neutron
 example_pop.set(
     C_logging_code="""
-if(stardata->star[0].stellar_type >= 2)    
+if(stardata->star[0].stellar_type >= 13)    
 {
     if (stardata->model.time < stardata->model.max_evolution_time)
     {
-        Printf("MY_STELLAR_DATA %30.12e %g %g %g %g\\n",
+        Printf("MY_STELLAR_DATA %30.12e %g %g %g %g %d\\n",
             // 
             stardata->model.time, // 1
-            stardata->star[0].mass, //2
-            stardata->star[0].pms_mass, //4
-            stardata->model.probability, //5
-            stardata->star[0].radius // 6
+            stardata->star[0].mass, // 2
+            stardata->common.zero_age.mass[0], // 4
+            stardata->model.probability, // 5
+            stardata->star[0].radius, // 6
+            stardata->star[0].stellar_type // 7
       );
     };
     /* Kill the simulation to save time */
@@ -120,7 +122,11 @@ if(stardata->star[0].stellar_type >= 2)
 )
 
 # Add grid variables
-resolution = {"M_1": 100}
+resolution = {
+    "M_1": 20, 
+    'q': 20, 
+    'per': 40
+}
 
 # Mass
 example_pop.add_grid_variable(
@@ -136,6 +142,38 @@ example_pop.add_grid_variable(
     condition="",  # Impose a condition on this grid variable. Mostly for a check for yourself
 )
 
+# # Mass ratio
+# test_pop.add_grid_variable(
+#     name="q",
+#     longname="Mass ratio",
+#     valuerange=["0.1/M_1", 1],
+#     resolution="{}".format(resolution['q']),
+#     spacingfunc="const(0.1/M_1, 1, {})".format(resolution['q']),
+#     probdist="flatsections(q, [{'min': 0.1/M_1, 'max': 1.0, 'height': 1}])",
+#     dphasevol="dq",
+#     precode="M_2 = q * M_1",
+#     parameter_name="M_2",
+#     condition="",  # Impose a condition on this grid variable. Mostly for a check for yourself    
+# )
+
+# # 
+# test_pop.add_grid_variable(
+#    name="log10per", # in days
+#    longname="log10(Orbital_Period)", 
+#    valuerange=[0.15, 5.5],
+#    resolution="{}".format(resolution["per"]),
+#    spacingfunc="const(0.15, 5.5, {})".format(resolution["per"]),
+#    precode="""orbital_period = 10** log10per
+# sep = calc_sep_from_period(M_1, M_2, orbital_period)
+# sep_min = calc_sep_from_period(M_1, M_2, 10**0.15)
+# sep_max = calc_sep_from_period(M_1, M_2, 10**5.5)""", 
+#    probdist="sana12(M_1, M_2, sep, orbital_period, sep_min, sep_max, math.log10(10**0.15), math.log10(10**5.5), -0.55)",
+#    parameter_name="orbital_period",
+#    dphasevol="dlog10per",
+# )
+
+
+
 # Exporting of all the settings can be done with .export_all_info()
 # on default it exports everything, but can be supressed by turning it off:
 #   population settings (bse_options, grid_options, custom_options), turn off with include_population
@@ -152,12 +190,12 @@ example_pop.export_all_info()
 
 ## Executing a single system
 ## This uses the M_1 orbital period etc set with the set function
-output = example_pop.evolve_single()
-print(output)
+# output = example_pop.evolve_single()
+# print(output)
 
 ## Executing a population
 ## This uses the values generated by the grid_variables
-# example_pop.evolve_population_mp_chunks() # TODO: update this function call
+example_pop.evolve() # TODO: update this function call
 
 # Wrapping up the results to an hdf5 file can be done by using the create_hdf5
 # (<directory containing data and settings>) This function takes the settings file
-- 
GitLab