Example use-case: predicting GAIA-observable populations of star systems

In this notebook we show how to convolve data to predict gaia-observable populations for the present time.

Several ingredients are necessary here:

population-synthesis results that contain counts of star systems at particular times
a Milky-Way galaxy star formation rate history model
a method to sample the positions (/distances to sun) and on-sky angles of the star systems
a model to model binary orbit properties like projected distance. We will make use of Dyad.
a model to assign detection probabilities for the gaia survey. We will make use of GaiaUnlimited

We will work towards a single epoch observation in this notebook. Further sophistications like multi-epoch observation of the same systems to identify their binary character is not covered here (yet). Moreover, we will use a single-metallicity population.

The convolution broadly is done as follows:

calculate the formation yields of star systems that would exit at present day
Calculate detection probabilities for the rest based on their position (either randomly assigned or motivated by a spatially-defined SFH) in the Milkyway and their system properties. For this we make use of GaiaUnlimited
Use this information to predict observable populations for the GAIA survey.

Some things to consider:

The total stellar mass in a galaxy is large, so if one tries to model too many types of stellar systems, the output will be large too.
To avoid being overloaded with data, you should either filter out systems of lesser interest before the convolution

[ ]:

import os
import json
import time
import copy
import h5py
import pkg_resources

import astropy.units as u
import numpy as np
import astropy.constants as const
import pandas as pd

from syntheticstellarpopconvolve import convolve, default_convolution_config, default_convolution_instruction
from syntheticstellarpopconvolve.general_functions import temp_dir, generate_boilerplate_outputfile
from syntheticstellarpopconvolve.convolution_by_sampling import (
    select_dict_entries_with_new_indices,
)

from syntheticstellarpopconvolve.usecase_notebook_utils.usecase_lisa_utils import get_mass_norm, sample_distances_simple, get_period, sample_distances_interpolated

TMP_DIR = temp_dir("code", "notebook_example_GAIA_distribution", clean_path=True)

# The flag below allows the user to run this notebook without the full data or starformation rate.
FULL_VERSION = os.getenv("EXAMPLE_USECASE_GAIA_VERSION", False)

post-convolution function

The post-convolution function is where the magic happens.

The convolution is configured to sample systems. That means we must use the sampled_indices field to fetch the system properties, and we can assign some key properties ourselves. Most notably, we need to assign:

distance relative to the sun
on-sky right-ascension and declination angles

TODO: consider using agama for this instead

Because we are dealing with binary systems we must also consider the binary orbital properties:

mean anomaly
projected separation
inclination of the binary system

[2]:

import functools
from syntheticstellarpopconvolve.usecase_notebook_utils.usecase_lisa_utils import precompute_radial_cdf, sample_distances_and_angles_interpolated

Hr = 4  # Radial scale length
inverse_cdf = precompute_radial_cdf(Hr)
bound_sample_distances_and_angles_interpolated = functools.partial(
    sample_distances_and_angles_interpolated,
    inverse_cdf=inverse_cdf
)

def post_convolution_function(
    config, sfr_dict, data_dict, time_bin_info_dict, convolution_results, convolution_instruction
):

    # unpack some data
    system_indices = convolution_results["indices"]
    local_indices = np.arange(len(system_indices))

    # sample distances and angles
    d, ra, dec = bound_sample_distances_and_angles_interpolated(NBin=len(local_indices_unmerged_systems_in_waveband), Hr=Hr, inverse_cdf=inverse_cdf)
    d, ra, dec = sample_distances_and_angles_interpolated(NBin=len(local_indices_unmerged_systems_in_waveband), Hr=Hr)

    # load systems as binaries
    # sample in-orbit angles, mean anomalies, and inclinations.

    # Filter out systems that are

    pass

[ ]: