# Example use case: Zero-age stellar luminosity function in binaries

In this notebook we compute the luminosity function of the zero-age main-sequence by running a population of binary stars using binary_c.

Before you go through this notebook, you should look at notebook_luminosity_function.ipynb which is for the - conceptually more simple - single stars.

We start by loading in some standard Python modules and the binary_c module.

    [1]:

    import os
    import math
    from binarycpython.utils.grid import Population

    # help(Population) # Uncomment this line to see the public functions of this object

## Setting up the Population object

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.

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.

    [2]:

    # Create population object
    population = Population()

    # If you want verbosity, set this before other things
    population.set(verbosity=1)

    # Setting values can be done via .set(<parameter_name>=<value>)
    # Values that are known to be binary_c_parameters are loaded into bse_options.
    # Those that are present in the default grid_options are set in grid_options
    # All other values that you set are put in a custom_options dict
    population.set(
        # binary_c physics options
        max_evolution_time=0.1,  # maximum stellar evolution time in Myr
     )

    # We can access the options through
    print("verbosity is", population.grid_options['verbosity'])

    adding: max_evolution_time=0.1 to BSE_options
    verbosity is 1

## Adding grid variables

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.

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.

To add a grid variable to the population object we use `population.add_grid_variable`

    [3]:

    # help(population.add_grid_variable)

All the distribution functions that we can use are stored in the `binarycpython.utils.distribution_functions` or `binarycpython/utils/distribution_functions.py` on git. If you uncomment the help statement below you can see which functions are available now:

    [4]:

    import binarycpython.utils.distribution_functions
    # help(binarycpython.utils.distribution_functions)

First let us set up some global variables that will be useful throughout.

- The resolution is the number of stars we simulate in our model population.

- The massrange is a list of the min and max masses

- The total_probability is the theoretical integral of a probability density function, i.e. 1.0.

- The binwidth sets the resolution of the final distribution. If set to 0.5, the bins in logL are 0.5dex wide.

    [5]:

    # Set resolution and mass range that we simulate
    resolution = {"M_1": 40} # start with resolution = 10, and increase later if you want "more accurate" data
    massrange = (0.07, 100.0) # we work with stars of mass 0.07 to 100 Msun
    total_probability = 1.0 # theoretical integral of the mass probability density function over all masses
    # distribution binwidths :
    # (log10) luminosity distribution
    binwidth = { 'luminosity' : 1.0 }

The next cell contains an example of adding the mass grid variable, sampling the phase space in linear mass M_1.

    [6]:

    # Set up the binary grid in "cubic" M1 - M2=q*M1 - log10 period space

    population = Population()

    # resolution on each side of the cube, with more stars for the primary mass
    nres = 10
    resolution = {"M_1": 4*nres,
                  "q": nres,
                  "per": nres}

    massrange = [0.07,100]
    logperrange = [0.15, 5.5]

    population.add_grid_variable(
        name="lnm1",
        longname="Primary mass",
        valuerange=massrange,
        samplerfunc="const(math.log({min}), math.log({max}), {res})".format(min=massrange[0],max=massrange[1],res=resolution["M_1"]),
        precode="M_1=math.exp(lnm1)",
        probdist="three_part_powerlaw(M_1, 0.1, 0.5, 1.0, 150, -1.3, -2.3, -2.3)*M_1",
        dphasevol="dlnm1",
        parameter_name="M_1",
        condition="",  # Impose a condition on this grid variable. Mostly for a check for yourself
    )

    # Mass ratio
    population.add_grid_variable(
         name="q",
         longname="Mass ratio",
         valuerange=["0.1/M_1", 1],
         samplerfunc="const({}/M_1, 1, {})".format(massrange[0],resolution['q']),
         probdist="flatsections(q, [{{'min': {}/M_1, 'max': 1.0, 'height': 1}}])".format(massrange[0]),
         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
     )

    # Orbital period
    population.add_grid_variable(
        name="log10per", # in days
        longname="log10(Orbital_Period)",
        valuerange=[0.15, 5.5],
        samplerfunc="const({}, {}, {})".format(logperrange[0],logperrange[1],resolution["per"]),
        precode="""orbital_period = 10.0 ** log10per
    sep = calc_sep_from_period(M_1, M_2, orbital_period)
    sep_min = calc_sep_from_period(M_1, M_2, 10**{})
    sep_max = calc_sep_from_period(M_1, M_2, 10**{})""".format(logperrange[0],logperrange[1]),
        probdist="sana12(M_1, M_2, sep, orbital_period, sep_min, sep_max, math.log10(10**{}), math.log10(10**{}), {})".format(logperrange[0],logperrange[1],-0.55),
        parameter_name="orbital_period",
        dphasevol="dlog10per",
     )

## Setting logging and handling the output

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.

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`).

In the code below we will set up both the custom logging and a parse function to handle that output.

    [7]:

    # Create custom logging statement
    #
    # we check that the model number is zero, i.e. we're on the first timestep (stars are born on the ZAMS)
    # we make sure that the stellar type is <= MAIN_SEQUENCE, i.e. the star is a main-sequence star
    # we also check that the time is 0.0 (this is not strictly required, but good to show how it is done)
    #
    # The
    #
    # The Printf statement does the outputting: note that the header string is ZERO_AGE_MAIN_SEQUENCE_STARn
    #
    # where:
    #
    # n = PRIMARY    = 0 is star 0 (primary star)
    # n = SECONDARY  = 1 is star 1 (secondary star)
    # n = UNRESOLVED = 2 is the unresolved system (both stars added)

    PRIMARY = 0
    SECONDARY = 1
    UNRESOLVED = 2

    custom_logging_statement = """
    // select ZAMS
    if(stardata->model.model_number == 0 &&
       stardata->model.time == 0)
    {
        // loop over the stars individually (equivalent to a resolved binary)
        Foreach_star(star)
        {
            // select main-sequence stars
            if(star->stellar_type <= MAIN_SEQUENCE)
            {
                /* Note that we use Printf - with a capital P! */
               Printf("ZERO_AGE_MAIN_SEQUENCE_STAR%d %30.12e %g %g %g %g\\n",
                      star->starnum,
                      stardata->model.time, // 1
                      stardata->common.zero_age.mass[0], // 2
                      star->mass, // 3
                      star->luminosity, // 4
                      stardata->model.probability // 5
               );
            }
        }

        // unresolved MS-MS binary
        if(stardata->star[0].stellar_type <= MAIN_SEQUENCE &&
           stardata->star[1].stellar_type <= MAIN_SEQUENCE)
        {
            Printf("ZERO_AGE_MAIN_SEQUENCE_STAR%d %30.12e %g %g %g %g\\n",
                      2,
                      stardata->model.time, // 1
                      stardata->common.zero_age.mass[0] + stardata->common.zero_age.mass[1], // 2
                      stardata->star[0].mass + stardata->star[1].mass, // 3
                      stardata->star[0].luminosity + stardata->star[1].luminosity, // 4
                      stardata->model.probability // 5
               );
        }
    }
    """

    population.set(
        C_logging_code=custom_logging_statement
    )

The parse function must now catch lines that start with “ZERO_AGE_MAIN_SEQUENCE_STAR” and process the associated data.

    [8]:

    # import the bin_data function so we can construct finite-resolution probability distributions
    # import the datalinedict to make a dictionary from each line of data from binary_c
    from binarycpython.utils.functions import bin_data,datalinedict
    import re

    def parse_function(self, output):
        """
        Example parse function
        """

        # list of the data items
        parameters = ["header", "time", "zams_mass", "mass", "luminosity", "probability"]

        # Loop over the output.
        for line in output.splitlines():

            # check if we match a ZERO_AGE_MAIN_SEQUENCE_STAR
            match = re.search('ZERO_AGE_MAIN_SEQUENCE_STAR(\d)',line)
            if match:
                nstar = match.group(1)
                #print("matched star",nstar)

                # obtain the line of data in dictionary form
                linedata = datalinedict(line,parameters)

                # bin the log10(luminosity) to the nearest 0.1dex
                binned_log_luminosity = bin_data(math.log10(linedata['luminosity']),
                                                 binwidth['luminosity'])

                # append the data to the results_dictionary
                self.grid_results['luminosity distribution'][int(nstar)][binned_log_luminosity] += linedata['probability']

                #print (self.grid_results)

        # verbose reporting
        #print("parse out results_dictionary=",self.grid_results)

    # Add the parsing function
    population.set(
        parse_function=parse_function,
    )

## Evolving the grid

Now that we configured all the main parts of the population object, we can actually run the population! Doing this is straightforward: `population.evolve()`

This will start up the processing of all the systems. We can control how many cores are used by settings `num_cores`. By setting the `verbosity` of the population object to a higher value we can get a lot of verbose information about the run, but for now we will set it to 0.

There are many grid_options that can lead to different behaviour of the evolution of the grid. Please do have a look at those: grid options docs, and try

    [9]:

    # set number of threads
    population.set(
        # verbose output is not required
        verbosity=1,
        # set number of threads (i.e. number of CPU cores we use)
        num_cores=4,
        )

    # Evolve the population - this is the slow, number-crunching step
    print("Running the population now, this may take a little while...")
    analytics = population.evolve()
    print("Done population run!")

    # Show the results (debugging)
    # print (population.grid_results)

    adding: num_cores=4 to grid_options
    Running the population now, this may take a little while...
    Creating and loading custom logging functionality
    Doing dry run to calculate total starcount and probability
    Generating grid code
    Generating grid code
    Saving grid code to grid_options
    Writing grid code to /tmp/binary_c_python-izzard/binary_c_grid_25014bc73b334765a1c09a4e4a97ed66.py [dry_run = True]
    Symlinked grid code to /tmp/binary_c_python-izzard/binary_c_grid-latest0
    Loading grid code function from /tmp/binary_c_python-izzard/binary_c_grid_25014bc73b334765a1c09a4e4a97ed66.py
    Grid code loaded
    Dry run of the grid
    Grid has handled 4000 stars with a total probability of 0.648566
    ****************************************
    * Total starcount for this run is 4000 *
    *     Total probability is 0.648566    *
    ****************************************

    [2021-11-01 09:56:53,685 DEBUG    Process-2] --- Setting up processor: process-0
    [2021-11-01 09:56:53,690 DEBUG    Process-3] --- Setting up processor: process-1
    [2021-11-01 09:56:53,694 DEBUG    Process-4] --- Setting up processor: process-2
    [2021-11-01 09:56:53,698 DEBUG    MainProcess] --- setting up the system_queue_filler now

    Generating grid code
    Generating grid code
    Saving grid code to grid_options
    Writing grid code to /tmp/binary_c_python-izzard/binary_c_grid_25014bc73b334765a1c09a4e4a97ed66.py [dry_run = False]
    Symlinked grid code to /tmp/binary_c_python-izzard/binary_c_grid-latest1

    [2021-11-01 09:56:53,707 DEBUG    Process-5] --- Setting up processor: process-3

    Loading grid code function from /tmp/binary_c_python-izzard/binary_c_grid_25014bc73b334765a1c09a4e4a97ed66.py
    Grid code loaded
    1354/4000  33.9% complete 09:56:57 ETA=    9.8s tpr=1.48e-02 ETF=09:57:07 mem:856.4MB M1=0.81 M2=0.48 P=3.6e+2
    1355/4000  33.9% complete 09:56:57 ETA=    7.3s tpr=1.11e-02 ETF=09:57:05 mem:856.4MB M1=0.81 M2=0.48 P=1.2e+3
    1449/4000  36.2% complete 09:57:02 ETA=    7.1s tpr=1.11e-02 ETF=09:57:10 mem:596.3MB M1=0.97 M2=0.48 P=1.7e+5
    1494/4000  37.4% complete 09:57:08 ETA=    2.3m tpr=2.20e-01 ETF=09:59:26 mem:603.7MB M1=0.97 M2=0.93 P=3.6e+2
    1539/4000  38.5% complete 09:57:13 ETA=    3.0m tpr=2.97e-01 ETF=10:00:16 mem:606.6MB M1=1.2 M2=0.45 P=1.7e+5
    1585/4000  39.6% complete 09:57:18 ETA=    3.4m tpr=3.36e-01 ETF=10:00:40 mem:607.6MB M1=1.2 M2=1 P=1.2e+3
    1624/4000  40.6% complete 09:57:23 ETA=    3.5m tpr=3.56e-01 ETF=10:00:54 mem:608.5MB M1=1.4 M2=0.4 P=3.6e+2
    1674/4000  41.9% complete 09:57:28 ETA=    3.7m tpr=3.79e-01 ETF=10:01:08 mem:611.5MB M1=1.4 M2=1.1 P=3.6e+2
    1703/4000  42.6% complete 09:57:33 ETA=    3.7m tpr=3.83e-01 ETF=10:01:13 mem:612.4MB M1=1.7 M2=0.15 P=1.1e+2
    1749/4000  43.7% complete 09:57:38 ETA=    3.9m tpr=4.11e-01 ETF=10:01:30 mem:613.0MB M1=1.7 M2=0.79 P=1.7e+5
    1780/4000  44.5% complete 09:57:43 ETA=    3.8m tpr=4.14e-01 ETF=10:01:33 mem:613.0MB M1=1.7 M2=1.4 P=2.6
    1814/4000  45.4% complete 09:57:48 ETA=    3.9m tpr=4.32e-01 ETF=10:01:44 mem:613.2MB M1=2 M2=0.36 P=3.6e+2
    1851/4000  46.3% complete 09:57:53 ETA=    4.0m tpr=4.44e-01 ETF=10:01:52 mem:613.6MB M1=2 M2=1.1 P=9
    1879/4000  47.0% complete 09:57:59 ETA=    4.0m tpr=4.52e-01 ETF=10:01:58 mem:614.8MB M1=2 M2=1.5 P=1.7e+5
    1910/4000  47.8% complete 09:58:04 ETA=    4.1m tpr=4.67e-01 ETF=10:02:08 mem:614.8MB M1=2.4 M2=0.42 P=2.6
    1947/4000  48.7% complete 09:58:09 ETA=    4.1m tpr=4.78e-01 ETF=10:02:14 mem:614.8MB M1=2.4 M2=1.1 P=1.5e+4
    1968/4000  49.2% complete 09:58:14 ETA=    4.1m tpr=4.82e-01 ETF=10:02:19 mem:616.0MB M1=2.4 M2=1.6 P=5e+4
    1989/4000  49.7% complete 09:58:19 ETA=    4.2m tpr=4.98e-01 ETF=10:02:30 mem:616.1MB M1=2.4 M2=2.1 P=1.7e+5
    2017/4000  50.4% complete 09:58:24 ETA=    4.3m tpr=5.15e-01 ETF=10:02:40 mem:617.1MB M1=2.9 M2=0.49 P=1.5e+4
    2048/4000  51.2% complete 09:58:29 ETA=    4.3m tpr=5.24e-01 ETF=10:02:45 mem:617.9MB M1=2.9 M2=1.3 P=5e+4
    2071/4000  51.8% complete 09:58:34 ETA=    4.3m tpr=5.31e-01 ETF=10:02:51 mem:619.7MB M1=2.9 M2=2.2 P=9
    2091/4000  52.3% complete 09:58:40 ETA=    4.3m tpr=5.42e-01 ETF=10:02:58 mem:620.5MB M1=2.9 M2=2.8 P=9
    2122/4000  53.0% complete 09:58:45 ETA=    4.3m tpr=5.55e-01 ETF=10:03:05 mem:620.5MB M1=3.5 M2=0.92 P=31
    2143/4000  53.6% complete 09:58:50 ETA=    4.3m tpr=5.59e-01 ETF=10:03:09 mem:623.3MB M1=3.5 M2=1.6 P=1.1e+2
    2164/4000  54.1% complete 09:58:55 ETA=    4.4m tpr=5.70e-01 ETF=10:03:16 mem:623.9MB M1=3.5 M2=2.3 P=3.6e+2
    2183/4000  54.6% complete 09:59:00 ETA=    4.4m tpr=5.80e-01 ETF=10:03:23 mem:624.0MB M1=3.5 M2=3 P=1.1e+2
    2206/4000  55.1% complete 09:59:05 ETA=    4.4m tpr=5.91e-01 ETF=10:03:30 mem:624.0MB M1=4.2 M2=0.27 P=4.2e+3
    2232/4000  55.8% complete 09:59:10 ETA=    4.4m tpr=5.99e-01 ETF=10:03:35 mem:624.0MB M1=4.2 M2=1.5 P=31
    2245/4000  56.1% complete 09:59:16 ETA=    4.4m tpr=6.06e-01 ETF=10:03:41 mem:624.1MB M1=4.2 M2=1.9 P=1.2e+3
    2258/4000  56.5% complete 09:59:21 ETA=    4.5m tpr=6.21e-01 ETF=10:03:51 mem:624.1MB M1=4.2 M2=2.3 P=5e+4
    2269/4000  56.7% complete 09:59:26 ETA=    4.6m tpr=6.36e-01 ETF=10:04:01 mem:626.1MB M1=4.2 M2=2.7 P=1.7e+5
    2282/4000  57.0% complete 09:59:31 ETA=    4.7m tpr=6.51e-01 ETF=10:04:11 mem:626.1MB M1=4.2 M2=3.6 P=31
    2300/4000  57.5% complete 09:59:36 ETA=    4.7m tpr=6.63e-01 ETF=10:04:18 mem:626.1MB M1=5 M2=0.32 P=2.6
    2329/4000  58.2% complete 09:59:41 ETA=    4.7m tpr=6.72e-01 ETF=10:04:22 mem:626.1MB M1=5 M2=1.3 P=1.7e+5
    2348/4000  58.7% complete 09:59:46 ETA=    4.6m tpr=6.73e-01 ETF=10:04:25 mem:626.2MB M1=5 M2=2.3 P=5e+4
    2365/4000  59.1% complete 09:59:52 ETA=    4.6m tpr=6.81e-01 ETF=10:04:30 mem:626.4MB M1=5 M2=3.3 P=1.2e+3
    2383/4000  59.6% complete 09:59:57 ETA=    4.7m tpr=6.90e-01 ETF=10:04:36 mem:626.4MB M1=5 M2=4.3 P=1.1e+2
    2400/4000  60.0% complete 10:00:02 ETA=    4.7m tpr=6.99e-01 ETF=10:04:42 mem:626.4MB M1=6 M2=0.37 P=2.6
    2423/4000  60.6% complete 10:00:07 ETA=    4.6m tpr=7.07e-01 ETF=10:04:46 mem:626.4MB M1=6 M2=1.6 P=1.1e+2
    2438/4000  61.0% complete 10:00:12 ETA=    4.6m tpr=7.10e-01 ETF=10:04:50 mem:628.5MB M1=6 M2=2.1 P=5e+4
    2454/4000  61.4% complete 10:00:18 ETA=    4.6m tpr=7.20e-01 ETF=10:04:56 mem:629.6MB M1=6 M2=3.3 P=3.6e+2
    2466/4000  61.6% complete 10:00:23 ETA=    4.7m tpr=7.29e-01 ETF=10:05:03 mem:629.6MB M1=6 M2=3.9 P=4.2e+3
    2477/4000  61.9% complete 10:00:28 ETA=    4.7m tpr=7.40e-01 ETF=10:05:10 mem:629.6MB M1=6 M2=4.5 P=1.5e+4
    2492/4000  62.3% complete 10:00:33 ETA=    4.7m tpr=7.51e-01 ETF=10:05:16 mem:630.7MB M1=6 M2=5.7 P=31
    2516/4000  62.9% complete 10:00:38 ETA=    4.7m tpr=7.59e-01 ETF=10:05:20 mem:630.7MB M1=7.2 M2=1.1 P=4.2e+3
    2537/4000  63.4% complete 10:00:44 ETA=    4.6m tpr=7.61e-01 ETF=10:05:22 mem:630.7MB M1=7.2 M2=2.6 P=1.5e+4
    2554/4000  63.9% complete 10:00:49 ETA=    4.6m tpr=7.65e-01 ETF=10:05:25 mem:630.8MB M1=7.2 M2=4 P=3.6e+2
    2570/4000  64.2% complete 10:00:54 ETA=    4.6m tpr=7.71e-01 ETF=10:05:29 mem:630.9MB M1=7.2 M2=5.4 P=2.6
    2590/4000  64.8% complete 10:00:59 ETA=    4.6m tpr=7.78e-01 ETF=10:05:33 mem:630.9MB M1=7.2 M2=6.8 P=2.6
    2622/4000  65.5% complete 10:01:04 ETA=    4.5m tpr=7.81e-01 ETF=10:05:33 mem:630.9MB M1=8.6 M2=2.2 P=31
    2636/4000  65.9% complete 10:01:09 ETA=    4.4m tpr=7.78e-01 ETF=10:05:34 mem:632.1MB M1=8.6 M2=3.1 P=4.2e+3
    2652/4000  66.3% complete 10:01:14 ETA=    4.4m tpr=7.87e-01 ETF=10:05:39 mem:632.1MB M1=8.6 M2=4.8 P=31
    2666/4000  66.7% complete 10:01:19 ETA=    4.7m tpr=8.37e-01 ETF=10:05:59 mem:634.6MB M1=8.6 M2=5.6 P=4.2e+3
    2691/4000  67.3% complete 10:01:25 ETA=    4.7m tpr=8.59e-01 ETF=10:06:06 mem:634.6MB M1=8.6 M2=8.2 P=9
    2733/4000  68.3% complete 10:01:31 ETA=    4.6m tpr=8.75e-01 ETF=10:06:08 mem:634.6MB M1=10 M2=3.7 P=1.1e+2
    2746/4000  68.7% complete 10:01:36 ETA=    4.6m tpr=8.83e-01 ETF=10:06:13 mem:634.6MB M1=10 M2=4.7 P=4.2e+3
    2757/4000  68.9% complete 10:01:42 ETA=    4.7m tpr=9.04e-01 ETF=10:06:22 mem:634.6MB M1=10 M2=5.7 P=1.5e+4
    2775/4000  69.4% complete 10:01:47 ETA=    4.8m tpr=9.37e-01 ETF=10:06:34 mem:634.8MB M1=10 M2=7.8 P=1.2e+3
    2795/4000  69.9% complete 10:01:52 ETA=    4.8m tpr=9.48e-01 ETF=10:06:38 mem:635.1MB M1=10 M2=9.8 P=1.2e+3
    2831/4000  70.8% complete 10:01:57 ETA=    4.7m tpr=9.71e-01 ETF=10:06:41 mem:636.3MB M1=12 M2=4.4 P=9
    2845/4000  71.1% complete 10:02:04 ETA=    4.7m tpr=9.67e-01 ETF=10:06:43 mem:636.3MB M1=12 M2=5.6 P=1.2e+3
    2861/4000  71.5% complete 10:02:09 ETA=    4.7m tpr=9.92e-01 ETF=10:06:52 mem:636.3MB M1=12 M2=8.1 P=9
    2875/4000  71.9% complete 10:02:16 ETA=    4.7m tpr=1.01e+00 ETF=10:07:01 mem:636.5MB M1=12 M2=9.3 P=1.2e+3
    2891/4000  72.3% complete 10:02:21 ETA=    4.8m tpr=1.03e+00 ETF=10:07:07 mem:636.5MB M1=12 M2=12 P=9
    2930/4000  73.2% complete 10:02:26 ETA=    4.7m tpr=1.05e+00 ETF=10:07:07 mem:636.7MB M1=15 M2=5.2 P=2.6
    2947/4000  73.7% complete 10:02:31 ETA=    4.6m tpr=1.05e+00 ETF=10:07:07 mem:636.7MB M1=15 M2=6.7 P=1.5e+4
    2965/4000  74.1% complete 10:02:37 ETA=    4.5m tpr=1.05e+00 ETF=10:07:09 mem:636.8MB M1=15 M2=9.7 P=1.2e+3
    2982/4000  74.5% complete 10:02:42 ETA=    4.5m tpr=1.06e+00 ETF=10:07:11 mem:637.3MB M1=15 M2=13 P=31

    [2021-11-01 10:02:46,515 DEBUG    MainProcess] --- Signaling stop to processes

    3010/4000  75.2% complete 10:02:47 ETA=    4.4m tpr=1.07e+00 ETF=10:07:12 mem:637.8MB M1=18 M2=2.7 P=2.6
    3035/4000  75.9% complete 10:02:52 ETA=    4.3m tpr=1.07e+00 ETF=10:07:11 mem:637.8MB M1=18 M2=6.3 P=1.2e+3
    3059/4000  76.5% complete 10:02:57 ETA=    4.2m tpr=1.07e+00 ETF=10:07:09 mem:637.8MB M1=18 M2=9.8 P=1.7e+5
    3078/4000  77.0% complete 10:03:02 ETA=    4.1m tpr=1.07e+00 ETF=10:07:08 mem:638.7MB M1=18 M2=13 P=5e+4
    3100/4000  77.5% complete 10:03:08 ETA=    4.0m tpr=1.08e+00 ETF=10:07:10 mem:638.9MB M1=21 M2=1.1 P=2.6
    3124/4000  78.1% complete 10:03:13 ETA=    3.9m tpr=1.08e+00 ETF=10:07:09 mem:639.1MB M1=21 M2=5.4 P=3.6e+2
    3148/4000  78.7% complete 10:03:18 ETA=    3.8m tpr=1.08e+00 ETF=10:07:08 mem:639.7MB M1=21 M2=9.7 P=5e+4
    3174/4000  79.3% complete 10:03:24 ETA=    3.7m tpr=1.07e+00 ETF=10:07:05 mem:640.1MB M1=21 M2=16 P=3.6e+2
    3197/4000  79.9% complete 10:03:29 ETA=    3.6m tpr=1.07e+00 ETF=10:07:04 mem:640.4MB M1=21 M2=20 P=1.5e+4
    3231/4000  80.8% complete 10:03:34 ETA=    3.4m tpr=1.07e+00 ETF=10:07:00 mem:640.9MB M1=26 M2=9 P=9
    3256/4000  81.4% complete 10:03:39 ETA=    3.3m tpr=1.05e+00 ETF=10:06:55 mem:640.9MB M1=26 M2=14 P=4.2e+3
    3273/4000  81.8% complete 10:03:44 ETA=    3.1m tpr=1.04e+00 ETF=10:06:53 mem:640.9MB M1=26 M2=19 P=1.1e+2
    3294/4000  82.3% complete 10:03:50 ETA=    3.0m tpr=1.03e+00 ETF=10:06:51 mem:641.4MB M1=26 M2=24 P=3.6e+2
    3321/4000  83.0% complete 10:03:55 ETA=    2.9m tpr=1.02e+00 ETF=10:06:48 mem:641.6MB M1=31 M2=7.7 P=9
    3348/4000  83.7% complete 10:04:00 ETA=    2.7m tpr=1.01e+00 ETF=10:06:45 mem:641.9MB M1=31 M2=14 P=5e+4
    3373/4000  84.3% complete 10:04:05 ETA=    2.6m tpr=1.01e+00 ETF=10:06:44 mem:641.9MB M1=31 M2=23 P=1.1e+2
    3394/4000  84.8% complete 10:04:10 ETA=    2.5m tpr=1.01e+00 ETF=10:06:43 mem:642.2MB M1=31 M2=29 P=3.6e+2
    3422/4000  85.5% complete 10:04:15 ETA=    2.4m tpr=1.00e+00 ETF=10:06:40 mem:642.5MB M1=37 M2=9.3 P=31
    3445/4000  86.1% complete 10:04:21 ETA=    2.3m tpr=9.93e-01 ETF=10:06:38 mem:642.7MB M1=37 M2=17 P=1.2e+3
    3464/4000  86.6% complete 10:04:26 ETA=    2.2m tpr=9.90e-01 ETF=10:06:39 mem:642.8MB M1=37 M2=24 P=3.6e+2
    3483/4000  87.1% complete 10:04:31 ETA=    2.1m tpr=9.96e-01 ETF=10:06:40 mem:642.8MB M1=37 M2=31 P=1.1e+2
    3509/4000  87.7% complete 10:04:37 ETA=    2.0m tpr=9.91e-01 ETF=10:06:38 mem:642.8MB M1=44 M2=2.3 P=1.7e+5
    3533/4000  88.3% complete 10:04:42 ETA=    1.9m tpr=9.81e-01 ETF=10:06:36 mem:642.8MB M1=44 M2=16 P=1.1e+2
    3550/4000  88.8% complete 10:04:47 ETA=    1.8m tpr=9.71e-01 ETF=10:06:36 mem:642.8MB M1=44 M2=24 P=2.6
    3568/4000  89.2% complete 10:04:52 ETA=    1.7m tpr=9.65e-01 ETF=10:06:37 mem:642.8MB M1=44 M2=29 P=5e+4
    3588/4000  89.7% complete 10:04:58 ETA=    1.7m tpr=9.64e-01 ETF=10:06:37 mem:643.1MB M1=44 M2=38 P=5e+4
    3622/4000  90.5% complete 10:05:03 ETA=    1.5m tpr=9.68e-01 ETF=10:06:34 mem:643.2MB M1=53 M2=13 P=313623/4000  90.6% complete 10:05:03 ETA=    1.5m tpr=9.68e-01 ETF=10:06:34 mem:643.2MB M1=53 M2=13 P=1.1e+2

    3640/4000  91.0% complete 10:05:08 ETA=    1.4m tpr=9.52e-01 ETF=10:06:33 mem:643.6MB M1=53 M2=24 P=2.6
    3657/4000  91.4% complete 10:05:13 ETA=    1.4m tpr=9.50e-01 ETF=10:06:34 mem:643.7MB M1=53 M2=29 P=1.5e+4
    3676/4000  91.9% complete 10:05:18 ETA=    1.3m tpr=9.53e-01 ETF=10:06:35 mem:643.7MB M1=53 M2=40 P=4.2e+3
    3706/4000  92.7% complete 10:05:23 ETA=    1.2m tpr=9.66e-01 ETF=10:06:34 mem:643.9MB M1=64 M2=3.2 P=4.2e+3
    3726/4000  93.2% complete 10:05:29 ETA=    1.1m tpr=9.50e-01 ETF=10:06:34 mem:644.5MB M1=64 M2=16 P=4.2e+3
    3746/4000  93.7% complete 10:05:34 ETA=    1.0m tpr=9.48e-01 ETF=10:06:34 mem:644.5MB M1=64 M2=29 P=4.2e+3
    3763/4000  94.1% complete 10:05:39 ETA=   55.8s tpr=9.42e-01 ETF=10:06:35 mem:644.5MB M1=64 M2=41 P=1.1e+2
    3786/4000  94.7% complete 10:05:45 ETA=   50.8s tpr=9.49e-01 ETF=10:06:35 mem:644.6MB M1=64 M2=54 P=4.2e+3
    3811/4000  95.3% complete 10:05:50 ETA=   45.5s tpr=9.63e-01 ETF=10:06:35 mem:645.0MB M1=76 M2=11 P=9
    3832/4000  95.8% complete 10:05:55 ETA=   39.9s tpr=9.51e-01 ETF=10:06:35 mem:645.0MB M1=76 M2=27 P=31
    3849/4000  96.2% complete 10:06:00 ETA=   35.5s tpr=9.42e-01 ETF=10:06:35 mem:645.1MB M1=76 M2=34 P=1.7e+5
    3875/4000  96.9% complete 10:06:05 ETA=   29.4s tpr=9.40e-01 ETF=10:06:34 mem:645.4MB M1=76 M2=57 P=1.2e+3
    3905/4000  97.6% complete 10:06:10 ETA=   22.2s tpr=9.36e-01 ETF=10:06:32 mem:645.5MB M1=91 M2=4.6 P=1.2e+3
    3930/4000  98.2% complete 10:06:15 ETA=   16.5s tpr=9.41e-01 ETF=10:06:32 mem:645.5MB M1=91 M2=32 P=2.6
    3931/4000  98.3% complete 10:06:15 ETA=   16.2s tpr=9.41e-01 ETF=10:06:31 mem:645.5MB M1=91 M2=32 P=9
    3954/4000  98.8% complete 10:06:20 ETA=   10.6s tpr=9.19e-01 ETF=10:06:31 mem:645.8MB M1=91 M2=50 P=3.6e+2
    3977/4000  99.4% complete 10:06:25 ETA=    5.2s tpr=9.06e-01 ETF=10:06:30 mem:645.8MB M1=91 M2=69 P=1.5e+4

    [2021-11-01 10:06:29,167 DEBUG    Process-3] --- Process-1 is finishing.

    process 1 free memory and return
    ****************************************************
    *                Process 1 finished:               *
    *  generator started at 2021-11-01T09:56:53.690194 *
    * generator finished at 2021-11-01T10:06:29.176751 *
    *                 total: 9m 35.49s                 *
    *         of which 9m 34.96s with binary_c         *
    *                 Ran 1001 systems                 *
    *       with a total probability of 0.160675       *
    *         This thread had 0 failing systems        *
    *       with a total failed probability of 0       *
    *   Skipped a total of 0 zero-probability systems  *
    *                                                  *
    ****************************************************

    process 1 queue put output_dict

    [2021-11-01 10:06:29,186 DEBUG    Process-3] --- Process-1 is finished.

    process 1 return

    [2021-11-01 10:06:29,342 DEBUG    Process-5] --- Process-3 is finishing.

    process 3 free memory and return
    ****************************************************
    *                Process 3 finished:               *
    *  generator started at 2021-11-01T09:56:53.706780 *
    * generator finished at 2021-11-01T10:06:29.345842 *
    *                 total: 9m 35.64s                 *
    *         of which 9m 35.06s with binary_c         *
    *                 Ran 1001 systems                 *
    *       with a total probability of 0.155662       *
    *         This thread had 0 failing systems        *
    *       with a total failed probability of 0       *
    *   Skipped a total of 0 zero-probability systems  *
    *                                                  *
    ****************************************************

    process 3 queue put output_dict

    [2021-11-01 10:06:29,350 DEBUG    Process-5] --- Process-3 is finished.

    process 3 return

    [2021-11-01 10:06:29,429 DEBUG    Process-2] --- Process-0 is finishing.

    process 0 free memory and return
    ****************************************************
    *                Process 0 finished:               *
    *  generator started at 2021-11-01T09:56:53.684890 *
    * generator finished at 2021-11-01T10:06:29.433207 *
    *                 total: 9m 35.75s                 *
    *         of which 9m 35.15s with binary_c         *
    *                 Ran 1025 systems                 *
    *       with a total probability of 0.162454       *
    *         This thread had 0 failing systems        *
    *       with a total failed probability of 0       *
    *   Skipped a total of 0 zero-probability systems  *
    *                                                  *
    ****************************************************

    process 0 queue put output_dict

    [2021-11-01 10:06:29,437 DEBUG    Process-2] --- Process-0 is finished.

    process 0 return

    [2021-11-01 10:06:29,449 DEBUG    Process-4] --- Process-2 is finishing.

    process 2 free memory and return
    ****************************************************
    *                Process 2 finished:               *
    *  generator started at 2021-11-01T09:56:53.694517 *
    * generator finished at 2021-11-01T10:06:29.453059 *
    *                 total: 9m 35.76s                 *
    *         of which 9m 35.25s with binary_c         *
    *                  Ran 973 systems                 *
    *       with a total probability of 0.169775       *
    *         This thread had 0 failing systems        *
    *       with a total failed probability of 0       *
    *   Skipped a total of 0 zero-probability systems  *
    *                                                  *
    ****************************************************

    process 2 queue put output_dict

    [2021-11-01 10:06:29,456 DEBUG    Process-4] --- Process-2 is finished.

    process 2 return
    ****************************************************************
    *     Population-25014bc73b334765a1c09a4e4a97ed66 finished!    *
    *              The total probability is 0.648566.              *
    *  It took a total of 9m 35.99s to run 4000 systems on 4 cores *
    *                   = 38m 23.97s of CPU time.                  *
    *                 Maximum memory use 856.406 MB                *
    ****************************************************************

    There were no errors found in this run.
    Done population run!

After the run is complete, some technical report on the run is returned. I stored that in `analytics`. As we can see below, this dictionary is like a status report of the evolution. Useful for e.g. debugging.

    [10]:

    print(analytics)

    {'population_name': '25014bc73b334765a1c09a4e4a97ed66', 'evolution_type': 'grid', 'failed_count': 0, 'failed_prob': 0, 'failed_systems_error_codes': [], 'errors_exceeded': False, 'errors_found': False, 'total_probability': 0.6485656144116352, 'total_count': 4000, 'start_timestamp': 1635760613.6602514, 'end_timestamp': 1635761189.652638, 'total_mass_run': 82563.09295167374, 'total_probability_weighted_mass_run': 0.6438124832773024, 'zero_prob_stars_skipped': 0}

    [11]:

    # make a plot of the luminosity distribution using Seaborn and Pandas
    import seaborn as sns
    import pandas as pd
    from binarycpython.utils.functions import pad_output_distribution

    # set up seaborn for use in the notebook
    sns.set(rc={'figure.figsize':(20,10)})
    sns.set_context("notebook",
                    font_scale=1.5,
                    rc={"lines.linewidth":2.5})

    titles = { 0 : "Primary",
               1 : "Secondary",
               2 : "Unresolved" }

    # choose to plot the
    # PRIMARY, SECONDARY or UNRESOLVED
    nstar = UNRESOLVED

    plots = {}

    # pad the distribution with zeros where data is missing
    for n in range(0,3):
        pad_output_distribution(population.grid_results['luminosity distribution'][n],
                                binwidth['luminosity'])
        plots[titles[n] + ' ZAMS luminosity distribution'] = population.grid_results['luminosity distribution'][n]

    # make pandas dataframe from our sorted dictionary of data
    plot_data = pd.DataFrame.from_dict(plots)

    # make the plot
    p = sns.lineplot(data=plot_data)
    p.set_xlabel("$\log_{10}$ ($L_\mathrm{ZAMS}$ / L$_{☉}$)")
    p.set_ylabel("Number of stars")
    p.set(yscale="log")

    [11]:

    [None]

You can see that the secondary stars are dimmer than the primaries - which you expect given they are lower in mass (by definition q=M2/M1<1).

Weirdly, in some places the primary distribution may exceed the unresolved distribution. This is a bit unphysical, but in this case is usually caused by limited resolution. If you increase the number of stars in the grid, this problem should go away (at a cost of more CPU time).

Things to try:

- Massive stars: can you see the effects of wind mass loss and rejuvenation in these stars?

- Alter the metallicity, does this make much of a difference?

- Change the binary fraction. Here we assume a 100% binary fraction, but a real population is a mixture of single and binary stars.

- How might you go about comparing these computed observations to real stars?

- What about evolved stars? Here we consider only the zero-age main sequnece. What about other main-sequence stars? What about stars in later phases of stellar evolution?

    [ ]: