diff --git a/examples/notebook_luminosity_function.ipynb b/examples/notebook_luminosity_function_single.ipynb
similarity index 96%
rename from examples/notebook_luminosity_function.ipynb
rename to examples/notebook_luminosity_function_single.ipynb
index 824f6301e362622d16d650db38e6cea9edfa6048..d2984104dc58d191dbef0321cb92d3179dd9adcb 100644
--- a/examples/notebook_luminosity_function.ipynb
+++ b/examples/notebook_luminosity_function_single.ipynb
@@ -7,7 +7,7 @@
"tags": []
},
"source": [
- "# Stellar luminosity function\n",
+ "# Zero-age stellar luminosity function\n",
"\n",
"In this notebook we compute the luminosity function of the zero-age main-sequence by running a population of single stars using binary_c. \n",
"\n",
@@ -30,7 +30,7 @@
"\n",
"TMP_DIR = temp_dir(\"notebooks\", \"notebook_luminosity\")\n",
"\n",
- "# help(Population) # Uncomment to see the public functions of this object"
+ "# help(Population) # Uncomment this line to see the public functions of this object"
]
},
{
@@ -39,14 +39,9 @@
"metadata": {},
"source": [
"## Setting up the Population object\n",
- "To set up and configure the population object we need to make an object instance of the `Population` object, and add configuration via the `.set()` function.\n",
+ "To set up and configure the population object we need to make a new instance of the `Population` object and configure it with the `.set()` function.\n",
"\n",
- "There are three categories of options that the Population object can set:\n",
- "- BSE (stellar evolution) options: these options will be used for the binary_c calls, and are recognized by comparing the arguments to a known list of available arguments of binary_c. To see which options are available, see section [`binary_c parameters` in the documentation](https://ri0005.pages.surrey.ac.uk/binary_c-python/binary_c_parameters.html). You can access these through `population.bse_options['<bse option name>']` after you have set them. \n",
- "\n",
- "- Grid options: these options will be used to configure the behaviour of the Population object. To see which options are available, see section [`Population grid code options` in the documentation](https://ri0005.pages.surrey.ac.uk/binary_c-python/grid_options_descriptions.html). They can be accessed via `population.grid_options['<grid option name>']` after you have set them. \n",
- "\n",
- "- Custom options: these options are not recognized as either of the above, so they will be stored in the custom_options, and can be accessed via `population.custom_options['<custom option name>']`"
+ "In our case, we only need to set the maximum evolution time to something short, because we care only about zero-age main sequence stars which have, by definition, age zero."
]
},
{
@@ -59,8 +54,8 @@
"name": "stdout",
"output_type": "stream",
"text": [
- "adding: max_evolution_time=15000 to BSE_options\n",
"adding: tmp_dir=/tmp/binary_c_python/notebooks/notebook_luminosity to grid_options\n",
+ "adding: max_evolution_time=0.1 to BSE_options\n",
"verbosity is 1\n"
]
}
@@ -78,7 +73,7 @@
"# All other values that you set are put in a custom_options dict\n",
"population.set(\n",
" # binary_c physics options\n",
- " max_evolution_time=15000, # maximum stellar evolution time in Myr\n",
+ " max_evolution_time=0.1, # maximum stellar evolution time in Myr\n",
" tmp_dir=TMP_DIR,\n",
")\n",
"\n",
@@ -94,14 +89,9 @@
"## Adding grid variables\n",
"The main purpose of the Population object is to handle the population synthesis side of running a set of stars. The main method to do this with binarycpython, as is the case with Perl binarygrid, is to use grid variables. These are loops over a predefined range of values, where a probability will be assigned to the systems based on the chosen probability distributions.\n",
"\n",
- "Usually we use either 1 mass grid variable, or a trio of mass, mass ratio and period (See below for full examples of all of these). We can, however, also add grid sampling for e.g. eccentricity, metallicity or other parameters. \n",
- "\n",
- "In some cases it could be easier to set up a for loop that sets that parameter and calls the evolve function several times, e.g. when you want to vary a prescription (usually a discrete, unweighted parameter) \n",
- "\n",
+ "Usually we use either 1 mass grid variable, or a trio of mass, mass ratio and period (other notebooks cover these examples). We can, however, also add grid sampling for e.g. eccentricity, metallicity or other parameters. \n",
"\n",
- "A notable special type of grid variable is that of the Moe & di Stefano 2017 dataset (see further down in the notebook).\n",
- "\n",
- "To add a grid variable to the population object we use `population.add_grid_variable` (see next cell)"
+ "To add a grid variable to the population object we use `population.add_grid_variable`"
]
},
{
@@ -136,12 +126,40 @@
"# help(binarycpython.utils.distribution_functions)"
]
},
+ {
+ "cell_type": "markdown",
+ "id": "2a9104fc-4136-4e53-8604-f24ad52fbe56",
+ "metadata": {},
+ "source": [
+ "First let us set up some global variables that will be useful throughout. \n",
+ "* The resolution is the number of stars we simulate in our model population.\n",
+ "* The massrange is a list of the min and max masses\n",
+ "* The total_probability is the theoretical integral of a probability density function, i.e. 1.0.\n",
+ "* The binwidth sets the resolution of the final distribution. If set to 0.5, the bins in log*L* are 0.5dex wide."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "aba3fe4e-18f2-4bb9-8e5c-4c6007ab038b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Set resolution and mass range that we simulate\n",
+ "resolution = {\"M_1\": 40} # start with resolution = 10, and increase later if you want \"more accurate\" data\n",
+ "massrange = (0.07, 100.0) # we work with stars of mass 0.07 to 100 Msun\n",
+ "total_probability = 1.0 # theoretical integral of the mass probability density function over all masses \n",
+ "# distribution binwidths : \n",
+ "# (log10) luminosity distribution\n",
+ "binwidth = { 'luminosity' : 0.5 }"
+ ]
+ },
{
"cell_type": "markdown",
"id": "1b3a007b-5c17-42a7-a981-7e268e6f545c",
"metadata": {},
"source": [
- "The next cell contains an example of adding the mass grid variable, but sampling in log mass. The commented grid variables are examples of the mass ratio sampling and the period sampling."
+ "The next cell contains an example of adding the mass grid variable, sampling the phase space in linear mass *M*_1."
]
},
{
@@ -151,11 +169,6 @@
"metadata": {},
"outputs": [],
"source": [
- "# Add grid variables\n",
- "resolution = {\"M_1\": 40} # start with resolution = 10, and increase later if you want \"more accurate\" data\n",
- "massrange = [0.07, 100.0] # we work with stars of mass 0.07 to 100 Msun\n",
- "total_probability = 1.0\n",
- "\n",
"# Mass\n",
"population = Population()\n",
"population.set(\n",
@@ -171,7 +184,7 @@
" dphasevol=\"dM_1\",\n",
" parameter_name=\"M_1\",\n",
" condition=\"\", # Impose a condition on this grid variable. Mostly for a check for yourself\n",
- ")\n"
+ ")"
]
},
{
@@ -180,21 +193,28 @@
"metadata": {},
"source": [
"## Setting logging and handling the output\n",
- "On default, binary_c will not output anything (except for 'SINGLE STAR LIFETIME'). It is up to us to determine what will be printed. We can either do that by hardcoding the print statements into `binary_c` (see documentation binary_c). Or, we can use the custom logging functionality of binarycpython (see notebook `notebook_custom_logging.ipynb`), which is faster to set up and requires no recompilation of binary_c, but is somewhat more limited in its functionality. \n",
+ "By default, binary_c will not output anything (except for 'SINGLE STAR LIFETIME'). It is up to us to determine what will be printed. We can either do that by hardcoding the print statements into `binary_c` (see documentation binary_c) or we can use the custom logging functionality of binarycpython (see notebook `notebook_custom_logging.ipynb`), which is faster to set up and requires no recompilation of binary_c, but is somewhat more limited in its functionality. For our current purposes, it works perfectly well.\n",
"\n",
"After configuring what will be printed, we need to make a function to parse the output. This can be done by setting the parse_function parameter in the population object (see also notebook `notebook_individual_systems.ipynb`). \n",
"\n",
- "In the code below we will set up both the custom logging, and a parse function to handle that output"
+ "In the code below we will set up both the custom logging and a parse function to handle that output."
]
},
{
"cell_type": "code",
- "execution_count": 6,
+ "execution_count": 18,
"id": "0c986215-93b1-4e30-ad79-f7c397e9ff7d",
"metadata": {},
"outputs": [],
"source": [
- "# Create custom logging statement: in this case we will log when the star turns into a compact object, and then terminate the evolution.\n",
+ "# Create custom logging statement\n",
+ "#\n",
+ "# we check that the model number is zero, i.e. we're on the first timestep (stars are born on the ZAMS)\n",
+ "# we make sure that the stellar type is <= MAIN_SEQUENCE, i.e. the star is a main-sequence star\n",
+ "# we also check that the time is 0.0 (this is not strictly required, but good to show how it is done)\n",
+ "#\n",
+ "# The Printf statement does the outputting: note that the header string is ZERO_AGE_MAIN_SEQUENCE_STAR\n",
+ "\n",
"custom_logging_statement = \"\"\"\n",
"if(stardata->model.model_number == 0 &&\n",
" stardata->star[0].stellar_type <= MAIN_SEQUENCE &&\n",
@@ -232,12 +252,8 @@
"outputs": [],
"source": [
"# import the bin_data function so we can construct finite-resolution probability distributions\n",
+ "# import the datalinedict to make a dictionary from each line of data from binary_c\n",
"from binarycpython.utils.functions import bin_data,datalinedict\n",
- " \n",
- "# distribution binwidths : \n",
- "# (log10) luminosity distribution\n",
- "binwidth = { 'luminosity' : 0.5 }\n",
- "\n",
"\n",
"def parse_function(self, output):\n",
" \"\"\"\n",
@@ -323,8 +339,8 @@
"# Evolve the population - this is the slow, number-crunching step\n",
"analytics = population.evolve() \n",
"\n",
- "# Show the results\n",
- "print (population.grid_results)"
+ "# Show the results (debugging)\n",
+ "# print (population.grid_results)"
]
},
{
@@ -403,8 +419,15 @@
"p = sns.lineplot(data=plot_data)\n",
"p.set_xlabel(\"$\\log_{10}$ ($L_\\mathrm{ZAMS}$ / L$_{☉}$)\")\n",
"p.set_ylabel(\"Number of stars\")\n",
- "p.set(yscale=\"log\")\n",
- "\n"
+ "p.set(yscale=\"log\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "7d7b275e-be92-4d59-b44d-ef6f24023cc3",
+ "metadata": {},
+ "source": [
+ "Does this look like a reasonable stellar luminosity function to you? The implication is that the most likely stellar luminosity is 10<sup>5.8</sup> L<sub>☉</sub>! Clearly, this is not very realistic... let's see what went wrong."
]
},
{
@@ -422,7 +445,7 @@
"source": [
"## ZAMS Luminosity distribution with the initial mass function\n",
"\n",
- "We now include an initial mass function as a three-part power law based on Kroupa (2001).\n"
+ "In the previous example, all the stars in our grid had an equal weighting. This is very unlikely to be true in reality: indeed, we know that low mass stars are far more likely than high mass stars. So we now include an initial mass function as a three-part power law based on Kroupa (2001). Kroupa's distribution is a three-part power law: we have a function that does this for us (it's very common to use power laws in astrophysics).\n"
]
},
{
@@ -467,8 +490,8 @@
"population.clean()\n",
"analytics = population.evolve() \n",
"\n",
- "# Show the results\n",
- "print (population.grid_results)"
+ "# Show the results (debugging)\n",
+ "# print (population.grid_results)"
]
},
{
@@ -521,7 +544,11 @@
"id": "0546f2f3-4732-4841-8ef3-565fbf6b9961",
"metadata": {},
"source": [
- "This distribution is peaked at low luminosity, as one expects from observations, but the resolution is clearly not great. We could simply increase the resolution to compensate, but this is very CPU intensive and a waste of time. Instead, let's try sampling the masses of the stars in a more smart way."
+ "This distribution is peaked at low luminosity, as one expects from observations, but the resolution is clearly not great because it's not smooth - it's spiky! \n",
+ "\n",
+ "If you noticed above, the total probability of the grid was about 0.2. Given that the total probability of a probability distribution function should be 1.0, this shows that our sampling is (very) poor. \n",
+ "\n",
+ "We could simply increase the resolution to compensate, but this is very CPU intensive and a complete waste of time and resources. Instead, let's try sampling the masses of the stars in a smarter way."
]
},
{
@@ -543,7 +570,7 @@
"metadata": {},
"outputs": [],
"source": [
- "# Rename the old variable because we want it to be called lnM_1 now\n",
+ "# Rename the old variable (M_1) because we want it to be called lnM_1 now\n",
"population.rename_grid_variable(\"M_1\",\"lnM_1\")"
]
},
@@ -560,15 +587,7 @@
"execution_count": 15,
"id": "108d470a-bb21-40b0-8387-2caa7ab0f923",
"metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "{'lnM_1': {'name': 'lnM_1', 'longname': 'Primary mass', 'valuerange': [0.07, 100.0], 'resolution': '40', 'spacingfunc': 'const(math.log(0.07), math.log(100.0), 40)', 'precode': 'M_1=math.exp(lnM_1)', 'probdist': 'three_part_powerlaw(M_1, 0.1, 0.5, 1.0, 150, -1.3, -2.3, -2.3)*M_1', 'dphasevol': 'dlnM_1', 'parameter_name': 'M_1', 'condition': '', 'gridtype': 'centred', 'branchpoint': 0, 'grid_variable_number': 0}}\n"
- ]
- }
- ],
+ "outputs": [],
"source": [
"# update the sampling, note that the IMF is dprob/dM1, and the phase \n",
"# space is now sampled in lnM1, so we multiply by M_1 to \n",
@@ -581,7 +600,7 @@
" parameter_name=\"M_1\",\n",
" precode=\"M_1=math.exp(lnM_1)\",\n",
")\n",
- "print(population.grid_options[\"_grid_variables\"])"
+ "# print(population.grid_options[\"_grid_variables\"]) # debugging"
]
},
{
@@ -612,8 +631,16 @@
"population.clean()\n",
"analytics = population.evolve() \n",
"\n",
- "# Show the results\n",
- "print (population.grid_results)"
+ "# Show the results (debugging)\n",
+ "# print (population.grid_results)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "182b1094-5057-4ccf-bac6-9b0e560ad4f6",
+ "metadata": {},
+ "source": [
+ "You should see that the total probability is very close to 1.0, as you would expect for a well-sampled grid. The total will never be exactly 1.0, but that is because we are running a simulation, not a perfect copy of reality."
]
},
{
@@ -659,9 +686,9 @@
"source": [
"Most stars are low mass red dwarfs, with small luminosities. Without the IMF weighting, our model population would have got this completely wrong! \n",
"\n",
- "As you increase the resolution, you will see this curve becomes even smoother. The wiggles in the curve are (usually) sampling artefacts because the curve should monotonically brighter than about log(*L*/L<sub>☉</sub>)=-2. \n",
+ "As you increase the resolution, you will see this curve becomes even smoother. The wiggles in the curve are (usually) sampling artefacts because the curve should monotonically brighten above about log(*L*/L<sub>☉</sub>)=-2. \n",
" \n",
- "Remember you can play with the binwidth too. If you want a very accurate distribution you need a narrow binwidth, but then you'll also need high resolution so lots of CPU time, hence cost, CO<sub>2</sub>, etc."
+ "Remember you can play with the binwidth too. If you want a very accurate distribution you need a narrow binwidth, but then you'll also need high resolution (lots of stars) so lots of CPU time, hence cost, CO<sub>2</sub>, etc."
]
}
],