[Theory paper](https://arxiv.org/pdf/1710.09839.pdf)
[GAIA paper](https://www.aanda.org/articles/aa/full_html/2016/11/aa29272-16/aa29272-16.html)

Note about units: The default unit for length is parsec, for mass is a solar mass, and for time is years.


```python
import sympy
sympy.init_printing()
import numpy
import pylab
import time
from sklearn.manifold import TSNE
```

Progress bar


```python
def log_progress(sequence, every=None, size=None, name='Items'):
    from ipywidgets import IntProgress, HTML, VBox
    from IPython.display import display

    is_iterator = False
    if size is None:
        try:
            size = len(sequence)
        except TypeError:
            is_iterator = True
    if size is not None:
        if every is None:
            if size <= 200:
                every = 1
            else:
                every = int(size / 200)     # every 0.5%
    else:
        assert every is not None, 'sequence is iterator, set every'

    if is_iterator:
        progress = IntProgress(min=0, max=1, value=1)
        progress.bar_style = 'info'
    else:
        progress = IntProgress(min=0, max=size, value=0)
    label = HTML()
    box = VBox(children=[label, progress])
    display(box)

    index = 0
    try:
        for index, record in enumerate(sequence, 1):
            if index == 1 or index % every == 0:
                if is_iterator:
                    label.value = '{name}: {index} / ?'.format(
                        name=name,
                        index=index
                    )
                else:
                    progress.value = index
                    label.value = u'{name}: {index} / {size}'.format(
                        name=name,
                        index=index,
                        size=size
                    )
            yield record
    except:
        progress.bar_style = 'danger'
        raise
    else:
        progress.bar_style = 'success'
        progress.value = index
        label.value = "{name}: {index}".format(
            name=name,
            index=str(index or '?')
        )
```

[Salpeter initial mass function](https://en.wikipedia.org/wiki/Initial_mass_function)


```python
M_star = sympy.Symbol('M_*', positive=True)
M_star_min = sympy.Symbol('M_{*min}', positive=True)
M_star_max = sympy.Symbol('M_{*max}', positive=True)
a_salpeter = sympy.Symbol('a_{salpeter}', positive=True)
C_salpeter = sympy.Symbol('C_{salpeter}', positive=True)
salpeter_imf = C_salpeter*M_star**(-a_salpeter)
salpeter_imf
```


```python
temp = salpeter_imf
temp = sympy.integrate(temp, M_star, conds='none')
C_salpeter_expr = sympy.solve(temp.subs(M_star, M_star_max) - temp.subs(M_star, M_star_min)-1, C_salpeter)[0]
C_salpeter_expr
```


```python
temp = salpeter_imf.subs(C_salpeter, C_salpeter_expr)
temp = sympy.integrate(temp, M_star, conds='none')
temp = temp - temp.subs(M_star, M_star_min)
cumulative_salpeter_imf = temp.simplify()
cumulative_salpeter_imf
```

What is the probability that a star is born with a mass in excess of 20 $M_{\odot}$?


```python
1-cumulative_salpeter_imf.subs({
    M_star: 20,
    M_star_min:0.08,
    M_star_max:100,
    a_salpeter:2.35
})
```

Distribution function for masses above 20 $M_{\odot}$


```python
xi = sympy.Symbol('xi', positive=True)
f = sympy.Symbol('f', positive=True)
temp = sympy.Eq(f, cumulative_salpeter_imf)
temp = sympy.solve(temp, M_star)[0].simplify()
temp = sympy.expand_power_base(temp.subs(M_star_max, xi*M_star_min), force=True).simplify()
temp = sympy.expand_power_base(temp, force=True).simplify()
temp = temp.subs(xi, M_star_max/M_star_min)
temp.subs({a_salpeter: 2.35,
           M_star_min: 20,
           M_star_max: 100
          })
```

This is a simplified version of the IMFs used in the theory paper (equation 1 and 2 in theory paper)


```python
def randomise_salpeter_bh_progenitor_mass(l):
    
    f = numpy.random.rand(l)
    
    return 20.0/(1.0-0.886*f)**0.741
```

Mass ratio (equation 3)


```python
def randomise_mass_ratio(M1_list):
    
    q_min = 0.08/M1_list
    f_list = numpy.random.rand(len(M1_list))
    return q_min + f_list*(1-q_min)
```

Semi major axis distribution (equation 6)


```python
def randomise_semi_major_axis(M1_list):
    
    f_list = numpy.random.rand(len(M1_list))
    A_min = 1e-8
    A_max = 1e-2
    #log10A_list = 8.0*(f_list - 1)
    return A_min*(A_max/A_min)**f_list
```

Position in the galaxy (equation 11)


```python
def randomise_height_above_galactic_disc(M1_list):
    
    f_list = numpy.random.rand(len(M1_list))
    return -250*numpy.log(1-f_list)
```


```python
def randomise_distance_from_galactic_centre(M1_list):
    
    f_list = numpy.random.rand(len(M1_list))
    return -3500*numpy.log(1-f_list)
```

Angle around the galactic rotation axis, relative to the line connecting the galactic centre to the sun


```python
def randomise_angle_gc_sol(M1_list):
    
    return numpy.pi*numpy.random.rand(len(M1_list))
```

Relation between mass and radius, using the empiric relation from [Demircan and Kahraman (1991)](http://adsabs.harvard.edu/abs/1991Ap%26SS.181..313D), page 319


```python
def mass2radius_demircan_kahraman(m):
    
    return 1.6*m**0.83*2.3e-8
```

Calculation of the ratio between the Roche lobe radius and the radius of the donor, according to the mass ratio $q$, using (equation 7)


```python
def roche_lobe_radius_ratio(q):
    
    return (0.6*q**-0.67+numpy.log(1+q**-0.33))/(0.49*q**-0.67)
```

Terminal semi major axis for binaries with comparable masses (equation 20)


```python
def calc_current_semi_major_axis_ratio_massive(q_list, k_list):
    
    a = 1.0
    b = 0.5
    c1 = a*(1-b)-2
    c2 = -a*(1-b)/b-2
    return (k_list+b*(1-k_list)+q_list)/(1+q_list)*k_list**c1*(b*(1-k_list)/q_list+1)**c2
```

Terminal semi major axis for binaries with large mass ratio (equation 23)


```python
def calc_current_semi_major_axis_ratio_wimpy(q_list, k_list, alphalambda=1.0):
    
    rl_list = 1.0/ roche_lobe_radius_ratio(q_list)
    
    return k_list/(2.0*(1.0-k_list)/(alphalambda*rl_list*q_list)+1.0)
```

Terminal semi major axis for both cases


```python
def calc_current_semi_major_axis(q_list, k_list, isma_list, alphalambda=1.0):
    
    return isma_list*numpy.where(
        q_list>0.5,
        calc_current_semi_major_axis_ratio_massive(q_list, k_list),
        calc_current_semi_major_axis_ratio_wimpy(q_list, k_list, alphalambda=alphalambda)
                               )
```

Relation between progenitor and black hole mass (equation 19)


```python
def calc_terminal_black_hole_mass(M1_list):
    
    return 8.0*numpy.log(M1_list - 19)/numpy.log(81.0)+2
```

Stellar classification

Based on [this table](http://www.pas.rochester.edu/~emamajek/EEM_dwarf_UBVIJHK_colors_Teff.txt)

Parse table


```python
with open('EEM_dwarf_UBVIJHK_colors_Teff.txt') as f:
    data = f.readlines()
table_lines = []
important_flag = False

for line in data:
    if '#SpT' in line:
        important_flag = not important_flag
        continue
    if important_flag:
        table_lines.append(line)

table_mass = []
table_vmag = []
for line in table_lines:
    mass_candidate_text = line.split()[19]
    mv_candidate_text = line.split()[4]
    if '...' in mass_candidate_text or '...' in mv_candidate_text:
        continue
    table_mass.append(float(mass_candidate_text))
    table_vmag.append(float(mv_candidate_text))
table_mass = numpy.array(table_mass)
table_vmag = numpy.array(table_vmag)
```


```python
def stellar_mass_to_vmag(mass):
    return numpy.interp(-mass, -table_mass, table_vmag)
%matplotlib inline
pylab.semilogx(table_mass, table_vmag,'.')
pylab.semilogx(table_mass, stellar_mass_to_vmag(table_mass));
```


```python
def absolute_to_apparent_magnitude(absolute, d_pc, extinction=False):
    if extinction:
        return absolute+5*(numpy.log10(d_pc)-1)+d_pc/1e3
    return absolute+5*(numpy.log10(d_pc)-1)
```

Gaia Sensitivity (section 8 in Gaia paper)


```python
def calc_gaia_z12p09(gmag):
    
    capped_gmag = numpy.clip(gmag, 12, 200)
    return numpy.clip(10**(0.4*(capped_gmag-15)),10**(0.4*(12.09-15)),1e30)

def calc_sigma_pomega(gmag):
    
    z = calc_gaia_z12p09(gmag)
    return numpy.sqrt(-1.63+680.8*z+32.7*z**2)

def calc_sigma_G(gmag):
    
    z = calc_gaia_z12p09(gmag)
    return 1.2e-3*numpy.sqrt(0.05*z**2+1.9*z+0.0002)
```

Synthetic stellar population


```python
ssp = {'primary mass':randomise_salpeter_bh_progenitor_mass(int(1e6))}
ssp['initial mass ratio'] = randomise_mass_ratio(ssp['primary mass'])
ssp['initial companion mass'] = ssp['primary mass']*ssp['initial mass ratio']
ssp['initial semi major axis'] = randomise_semi_major_axis(ssp['primary mass'])
ssp['z'] = randomise_height_above_galactic_disc(ssp['primary mass'])
ssp['r'] = randomise_distance_from_galactic_centre(ssp['primary mass'])
ssp['angle'] = randomise_angle_gc_sol(ssp['primary mass'])
ssp['black hole mass'] = calc_terminal_black_hole_mass(ssp['primary mass'])
#ssp['black hole mass'] = 0.2*ssp['primary mass']
ssp['terminal companion mass'] = (ssp['initial companion mass'] + 
                                  numpy.where(ssp['initial mass ratio']>0.5,
                                              0.5*(ssp['primary mass']-ssp['black hole mass']),
                                              numpy.zeros_like(ssp['primary mass'])))
ssp['terminal mass ratio'] = ssp['black hole mass']/ssp['terminal companion mass']
ssp['age'] = numpy.random.rand(len(ssp['primary mass'])) # Age / galactic age
ssp['lifetime'] = ssp['terminal companion mass']**-2.5 - ssp['primary mass']**-2.5 # Normalised by the age of the galaxy
ssp['companion mass'] = ssp['primary mass']*ssp['initial mass ratio']
ssp['companion radius'] = mass2radius_demircan_kahraman(ssp['companion mass'])
ssp['initial primary radius'] = mass2radius_demircan_kahraman(ssp['primary mass'])
ssp['initial roche radius'] = roche_lobe_radius_ratio(ssp['initial mass ratio'])*ssp['initial primary radius']
ssp['terminal roche radius'] = roche_lobe_radius_ratio(1.0/ssp['terminal mass ratio'])*ssp['companion radius']
ssp['terminal semi major axis'] = calc_current_semi_major_axis(
    ssp['initial mass ratio'],
    ssp['black hole mass']/ssp['primary mass'],
    ssp['initial semi major axis'],
    alphalambda=1.0
)
ssp['terminal period'] = 9.4e7*(ssp['terminal semi major axis']**1.5*
                                (ssp['companion mass']+ssp['black hole mass'])**(-0.5))
ssp['companion absolute magnitude'] = stellar_mass_to_vmag(ssp['primary mass']*ssp['initial mass ratio'])
ssp['distance from earth'] = numpy.sqrt(ssp['z']**2+8000**2+ssp['r']**2-2*8000*ssp['r']*numpy.cos(ssp['angle']))
ssp['companion apparent magnitude'] = absolute_to_apparent_magnitude(ssp['companion absolute magnitude'],
                                                                     ssp['distance from earth'],
                                                                    extinction=True)
ssp['sigma pomega'] = calc_sigma_pomega(ssp['companion apparent magnitude'])
ssp['parallax'] = 1e6/ssp['distance from earth'] # in uas
ssp['sigma G'] = calc_sigma_G(ssp['companion apparent magnitude'])
ssp['statistical weight'] = numpy.clip(ssp['lifetime'], 0, 1)
ssp_no_kick = ssp
```

Count how many will be visible by Gaia


```python
filter_list = [ssp['terminal semi major axis']>ssp['terminal roche radius'],
               ssp['terminal period']<5,
               ssp['terminal period']>0.137,
               #ssp['initial semi major axis']<3.2e-5,
               ssp['black hole mass']>3.0/(1.0-ssp['sigma pomega']/ssp['parallax']-ssp['sigma G']-ssp['sigma pomega']*ssp['distance from earth']/2e11/ssp['terminal semi major axis']),
               ssp['companion apparent magnitude']<20,
               ssp['terminal semi major axis']>10*(1+ssp['initial mass ratio']/0.2)*ssp['sigma pomega']*ssp['distance from earth']*5e-12
              ]
mask = numpy.all(filter_list,axis=0)
gaia_no_kick_mask = mask
ssp['statistical weight'][mask].sum()
```

Statistics of detected black holes - luminous companion binaries


```python
%matplotlib inline
import matplotlib
font = {'size':22}
matplotlib.rc('font',**font)

pylab.hist(ssp['black hole mass'][mask],20, weights=ssp['statistical weight'][mask])
pylab.xlabel('Black hole mass [$M_{\odot}$]');
pylab.ylabel('Number')
pass
```


```python
%matplotlib inline
import matplotlib
font = {'size':22}
matplotlib.rc('font',**font)

pylab.hist(ssp['initial semi major axis'][mask],20, weights=ssp['statistical weight'][mask])
#pylab.xlabel('Black hole mass [$M_{\odot}$]');
pylab.ylabel('Number')
pylab.xscale('log')
pass
```


```python
pylab.hist(ssp['terminal period'][mask],20, weights=ssp['statistical weight'][mask])
pylab.xlabel('Terminal period [years]');
pylab.ylabel('Number')
pass
```


```python
pylab.hist(ssp['distance from earth'][mask]/1000,20, weights=ssp['statistical weight'][mask])
pylab.xlabel('Distance from earth [kpc]');
pylab.ylabel('Number')
pass
```


```python
pylab.hist(ssp['companion mass'][mask],20, weights=ssp['statistical weight'][mask])
pylab.xscale('log')
pylab.xlabel(r'Companion mass [$M_{\odot}$]');
pylab.ylabel('Number')
pass
```


```python
pylab.hist(ssp['terminal semi major axis'][mask]/ssp['distance from earth'][mask]/4.8e-9,50,
          weights=ssp['statistical weight'][mask])
pylab.xscale('log')
pylab.xlabel(r'Angular separation [mas]')
pylab.ylabel('Number')
pass
```


```python
pylab.scatter(ssp['companion mass'][mask],
              ssp['terminal semi major axis'][mask]/ssp['distance from earth'][mask]/4.8e-9
             )
pylab.xscale('log')
pylab.yscale('log')
pylab.xlabel('Companion mass[$M_{\odot}$]')
pylab.ylabel('Angular separation [mas]')
pass
```

Count how many will be seen by [Hipparcos](https://arxiv.org/pdf/1202.2064.pdf)


```python
filter_list = [ssp['terminal semi major axis']>ssp['terminal roche radius'],
               ssp['terminal period']<5,
               ssp['terminal period']>0.137,
               ssp['initial semi major axis']<3.2e-5,
               ssp['black hole mass']>
               3.0/(1.0-100*ssp['sigma pomega']/ssp['parallax']-
                    100*ssp['sigma pomega']*ssp['distance from earth']/2e11/ssp['terminal semi major axis']),
               ssp['companion apparent magnitude']<12.4,
               ssp['terminal semi major axis']>
               10*(1+ssp['initial mass ratio']/0.2)*100*ssp['sigma pomega']*ssp['distance from earth']*5e-12
              ]
mask = numpy.all(filter_list, axis=0)
ssp['statistical weight'][mask].sum()
```

# Inefficient Common Envelope


```python
ssp = {'primary mass':randomise_salpeter_bh_progenitor_mass(int(1e6))}
ssp['initial mass ratio'] = randomise_mass_ratio(ssp['primary mass'])
ssp['initial companion mass'] = ssp['primary mass']*ssp['initial mass ratio']
ssp['initial semi major axis'] = randomise_semi_major_axis(ssp['primary mass'])
ssp['z'] = randomise_height_above_galactic_disc(ssp['primary mass'])
ssp['r'] = randomise_distance_from_galactic_centre(ssp['primary mass'])
ssp['angle'] = randomise_angle_gc_sol(ssp['primary mass'])
ssp['black hole mass'] = calc_terminal_black_hole_mass(ssp['primary mass'])
#ssp['black hole mass'] = 0.2*ssp['primary mass']
ssp['terminal companion mass'] = (ssp['initial companion mass'] + 
                                  numpy.where(ssp['initial mass ratio']>0.5,
                                              0.5*(ssp['primary mass']-ssp['black hole mass']),
                                              numpy.zeros_like(ssp['primary mass'])))
ssp['terminal mass ratio'] = ssp['black hole mass']/ssp['terminal companion mass']
ssp['age'] = numpy.random.rand(len(ssp['primary mass'])) # Age / galactic age
ssp['lifetime'] = ssp['terminal companion mass']**-2.5 - ssp['primary mass']**-2.5 # Normalised by the age of the galaxy
ssp['companion mass'] = ssp['primary mass']*ssp['initial mass ratio']
ssp['companion radius'] = mass2radius_demircan_kahraman(ssp['companion mass'])
ssp['initial primary radius'] = mass2radius_demircan_kahraman(ssp['primary mass'])
ssp['initial roche radius'] = roche_lobe_radius_ratio(ssp['initial mass ratio'])*ssp['initial primary radius']
ssp['terminal roche radius'] = roche_lobe_radius_ratio(1.0/ssp['terminal mass ratio'])*ssp['companion radius']
ssp['terminal semi major axis'] = calc_current_semi_major_axis(
    ssp['initial mass ratio'],
    ssp['black hole mass']/ssp['primary mass'],
    ssp['initial semi major axis'],
    alphalambda=0.1
)
ssp['terminal period'] = 9.4e7*(ssp['terminal semi major axis']**1.5*
                                (ssp['companion mass']+ssp['black hole mass'])**(-0.5))
ssp['companion absolute magnitude'] = stellar_mass_to_vmag(ssp['primary mass']*ssp['initial mass ratio'])
ssp['distance from earth'] = numpy.sqrt(ssp['z']**2+8000**2+ssp['r']**2-2*8000*ssp['r']*numpy.cos(ssp['angle']))
ssp['companion apparent magnitude'] = absolute_to_apparent_magnitude(ssp['companion absolute magnitude'],
                                                                     ssp['distance from earth'],
                                                                    extinction=True)
ssp['sigma pomega'] = calc_sigma_pomega(ssp['companion apparent magnitude'])
ssp['parallax'] = 1e6/ssp['distance from earth'] # in uas
ssp['sigma G'] = calc_sigma_G(ssp['companion apparent magnitude'])
ssp['statistical weight'] = numpy.clip(ssp['lifetime'], 0, 1)
ssp_inefficient = ssp
```

Count how many will be detected by gaia


```python
filter_list = [ssp['terminal semi major axis']>ssp['terminal roche radius'],
               ssp['terminal period']<5,
               ssp['terminal period']>0.137,
               ssp['black hole mass']>3.0/(1.0-ssp['sigma pomega']/ssp['parallax']-ssp['sigma G']-ssp['sigma pomega']*ssp['distance from earth']/2e11/ssp['terminal semi major axis']),
               ssp['companion apparent magnitude']<20,
               ssp['terminal semi major axis']>10*(1+ssp['initial mass ratio']/0.2)*ssp['sigma pomega']*ssp['distance from earth']*5e-12
              ]
mask = numpy.all(filter_list,axis=0)
gaia_inefficient_mask = mask
ssp['statistical weight'][mask].sum()
```

Count how many will be detected by hipparcos


```python
filter_list = [ssp['terminal semi major axis']>ssp['terminal roche radius'],
               ssp['terminal period']<5,
               ssp['terminal period']>0.137,
               ssp['black hole mass']>
               3.0/(1.0-100*ssp['sigma pomega']/ssp['parallax']-
                    100*ssp['sigma pomega']*ssp['distance from earth']/2e11/ssp['terminal semi major axis']),
               ssp['companion apparent magnitude']<12.4,
               ssp['terminal semi major axis']>
               10*(1+ssp['initial mass ratio']/0.2)*100*ssp['sigma pomega']*ssp['distance from earth']*5e-12
              ]
mask = numpy.all(filter_list, axis=0)
ssp['statistical weight'][mask].sum()
```

# Natal Kicks

## Same velocity as neutron stars


```python
def kick_semi_major_axis_ratio(vk_list, vi_list, chi, mu_list):
    
    y = vk_list/vi_list
    return 2-chi*(1+2*mu_list*y+y**2)
```


```python
def kick_eccentricity(vk_list, vi_list, chi, mu_list, phi_list):
    
    sma_ratio = kick_semi_major_axis_ratio(vk_list,
                                          vi_list,
                                          chi,
                                          mu_list)
    mu = mu_list
    y = vk_list/vi_list
    temp = (1+mu*y)**2-(mu**2-1)*y**2*numpy.sin(phi_list)**2
    
    return numpy.sqrt(1-sma_ratio*temp)
```


```python
ssp = {'primary mass':randomise_salpeter_bh_progenitor_mass(int(1e6))}
ssp['initial mass ratio'] = randomise_mass_ratio(ssp['primary mass'])
ssp['initial companion mass'] = ssp['primary mass']*ssp['initial mass ratio']
ssp['initial semi major axis'] = randomise_semi_major_axis(ssp['primary mass'])
ssp['z'] = randomise_height_above_galactic_disc(ssp['primary mass'])
ssp['r'] = randomise_distance_from_galactic_centre(ssp['primary mass'])
ssp['angle'] = randomise_angle_gc_sol(ssp['primary mass'])
ssp['black hole mass'] = calc_terminal_black_hole_mass(ssp['primary mass'])
ssp['age'] = numpy.random.rand(len(ssp['primary mass'])) # Age / galactic age
ssp['companion mass'] = ssp['primary mass']*ssp['initial mass ratio']
ssp['companion radius'] = mass2radius_demircan_kahraman(ssp['companion mass'])
ssp['initial primary radius'] = mass2radius_demircan_kahraman(ssp['primary mass'])
ssp['initial roche radius'] = roche_lobe_radius_ratio(ssp['initial mass ratio'])*ssp['initial primary radius']
ssp['terminal companion mass'] = (ssp['initial companion mass'] + 
                                  numpy.where(ssp['initial mass ratio']>0.5,
                                              0.5*(ssp['primary mass']-ssp['black hole mass']),
                                              numpy.zeros_like(ssp['primary mass'])))
ssp['pre kick semi major axis'] = calc_current_semi_major_axis(
    ssp['initial mass ratio'],
    ssp['black hole mass']/ssp['primary mass'],
    ssp['initial semi major axis'],
    alphalambda=1.0
)
ssp['terminal mass ratio'] = ssp['black hole mass']/ssp['terminal companion mass']
ssp['lifetime'] = ssp['terminal companion mass']**-2.5 - ssp['primary mass']**-2.5 # Normalised by the age of the galaxy
ssp['terminal roche radius'] = roche_lobe_radius_ratio(1.0/ssp['terminal mass ratio'])*ssp['companion radius']

# Kicks
ssp['kick mu'] = numpy.random.rand(len(ssp['primary mass']))*2-1
ssp['kick phi'] = numpy.random.rand(len(ssp['primary mass']))*numpy.pi
ssp['kick v_k'] = numpy.absolute(numpy.random.normal(0, 265, len(ssp['primary mass']))) # km/s
ssp['kick v_i'] = 0.07*numpy.sqrt((ssp['primary mass'] + ssp['companion mass'])/ssp['pre kick semi major axis']) # km/s
ssp['kick chi'] = (ssp['primary mass']+ssp['initial companion mass'])/(ssp['black hole mass'] + ssp['terminal companion mass'])

ssp['terminal semi major axis'] = ssp['pre kick semi major axis']/(
    kick_semi_major_axis_ratio(ssp['kick v_k'], 
                               ssp['kick v_i'], 
                               ssp['kick chi'],
                               ssp['kick mu'])
)

ssp['terminal eccentricity'] = kick_eccentricity(ssp['kick v_k'],
                                                 ssp['kick v_i'],
                                                 ssp['kick chi'], 
                                                 ssp['kick mu'], 
                                                 ssp['kick phi'])

ssp['terminal period'] = 9.4e7*(ssp['terminal semi major axis']**1.5*
                                ssp['primary mass']**-0.5*
                                (0.2+ssp['initial mass ratio'])**-0.5)
ssp['companion absolute magnitude'] = stellar_mass_to_vmag(ssp['primary mass']*ssp['initial mass ratio'])
ssp['distance from earth'] = numpy.sqrt(ssp['z']**2+8000**2+ssp['r']**2-2*8000*ssp['r']*numpy.cos(ssp['angle']))
ssp['companion apparent magnitude'] = absolute_to_apparent_magnitude(ssp['companion absolute magnitude'],
                                                                     ssp['distance from earth'],
                                                                    extinction=True)
ssp['sigma pomega'] = calc_sigma_pomega(ssp['companion apparent magnitude'])
ssp['parallax'] = 1e6/ssp['distance from earth'] # in uas
ssp['sigma G'] = calc_sigma_G(ssp['companion apparent magnitude'])
ssp['statistical weight'] = numpy.clip(ssp['lifetime'],0,1)
ssp_kick_velocity=ssp
```

    /usr/local/lib/python2.7/dist-packages/ipykernel_launcher.py:49: RuntimeWarning: invalid value encountered in power


Count how many will be seen by hipparcos


```python
filter_list = [ssp['terminal semi major axis']>0,
               ssp['terminal eccentricity'] <0.9,
               ssp['lifetime'] > ssp['age'],
               ssp['terminal semi major axis']>ssp['terminal roche radius'],
               ssp['terminal period']<5,
               ssp['terminal period']>0.137,
               ssp['black hole mass']>
               3.0/(1.0-100*ssp['sigma pomega']/ssp['parallax']-
                    100*ssp['sigma pomega']*ssp['distance from earth']/2e11/ssp['terminal semi major axis']),
               ssp['companion apparent magnitude']<12.4,
               ssp['terminal semi major axis']>
               10*(1+ssp['initial mass ratio']/0.2)*100*ssp['sigma pomega']*ssp['distance from earth']*5e-12
              ]
mask = numpy.all(filter_list,axis=0)
ssp['statistical weight'][mask].sum()
```

Count how many will be seen by Gaia


```python
filter_list = [
    ssp['terminal semi major axis'] > 0,
    ssp['terminal eccentricity'] <0.9,
    ssp['terminal semi major axis']>ssp['terminal roche radius'],
    ssp['terminal period']<5,
    ssp['terminal period']>0.137,
    ssp['black hole mass']>
    3.0/(1.0-ssp['sigma pomega']/ssp['parallax']-ssp['sigma G']-
         ssp['sigma pomega']*ssp['distance from earth']/2e11/ssp['terminal semi major axis']),
    ssp['companion apparent magnitude']<20,
    ssp['terminal semi major axis']>
    10*(1+ssp['initial mass ratio']/0.2)*ssp['sigma pomega']*ssp['distance from earth']*5e-12
]
mask = numpy.all(filter_list, axis=0)
gaia_kick_velocity_mask = mask
ssp['statistical weight'][mask].sum()
```


```python
%matplotlib inline
import matplotlib
font = {'size':22}
matplotlib.rc('font',**font)

pylab.hist(ssp['black hole mass'][mask],10, weights=ssp['statistical weight'][mask])
pylab.xlabel('Black hole mass [$M_{\odot}$]');
pylab.ylabel('Number')
pass
```


```python
pylab.hist(ssp['distance from earth'][mask]/1000,10, weights=ssp['statistical weight'][mask])
pylab.xlabel('Distance from earth [kpc]');
pylab.ylabel('Number')
pass
```


```python
pylab.hist(ssp['companion mass'][mask],10, weights=ssp['statistical weight'][mask])
pylab.xscale('log')
pylab.xlabel(r'Companion mass [$M_{\odot}$]');
pylab.ylabel('Number')
pass
```


```python
pylab.hist(ssp['terminal semi major axis'][mask]/ssp['distance from earth'][mask]/4.8e-9,10,
          weights=ssp['statistical weight'][mask])
pylab.xscale('log')
pylab.xlabel(r'Angular separation [mas]')
pylab.ylabel('Number')
pass
```


```python
pylab.scatter(ssp['companion mass'][mask],
              ssp['terminal semi major axis'][mask]/ssp['distance from earth'][mask]/4.8e-9
             )
pylab.xscale('log')
pylab.yscale('log')
pylab.xlabel('Companion mass[$M_{\odot}$]')
pylab.ylabel('Angular separation [mas]')
pass
```

## Momentum conservatoin


```python
ssp = {'primary mass':randomise_salpeter_bh_progenitor_mass(int(1e6))}
ssp['initial mass ratio'] = randomise_mass_ratio(ssp['primary mass'])
ssp['initial companion mass'] = ssp['primary mass']*ssp['initial mass ratio']
ssp['initial semi major axis'] = randomise_semi_major_axis(ssp['primary mass'])
ssp['z'] = randomise_height_above_galactic_disc(ssp['primary mass'])
ssp['r'] = randomise_distance_from_galactic_centre(ssp['primary mass'])
ssp['angle'] = randomise_angle_gc_sol(ssp['primary mass'])
ssp['black hole mass'] = calc_terminal_black_hole_mass(ssp['primary mass'])
ssp['age'] = numpy.random.rand(len(ssp['primary mass'])) # Age / galactic age
ssp['companion mass'] = ssp['primary mass']*ssp['initial mass ratio']
ssp['companion radius'] = mass2radius_demircan_kahraman(ssp['companion mass'])
ssp['initial primary radius'] = mass2radius_demircan_kahraman(ssp['primary mass'])
ssp['initial roche radius'] = roche_lobe_radius_ratio(ssp['initial mass ratio'])*ssp['initial primary radius']
ssp['terminal companion mass'] = (ssp['initial companion mass'] + 
                                  numpy.where(ssp['initial mass ratio']>0.5,
                                              0.5*(ssp['primary mass']-ssp['black hole mass']),
                                              numpy.zeros_like(ssp['primary mass'])))
ssp['pre kick semi major axis'] = calc_current_semi_major_axis(
    ssp['initial mass ratio'],
    ssp['black hole mass']/ssp['primary mass'],
    ssp['initial semi major axis'],
    alphalambda=1.0
)
ssp['terminal mass ratio'] = ssp['black hole mass']/ssp['terminal companion mass']
ssp['lifetime'] = ssp['terminal companion mass']**-2.5 - ssp['primary mass']**-2.5 # Normalised by the age of the galaxy
ssp['terminal roche radius'] = roche_lobe_radius_ratio(1.0/ssp['terminal mass ratio'])*ssp['companion radius']

# Kicks
ssp['kick mu'] = numpy.random.rand(len(ssp['primary mass']))*2-1
ssp['kick phi'] = numpy.random.rand(len(ssp['primary mass']))*numpy.pi
ssp['kick v_k'] = numpy.absolute(numpy.random.normal(0, 265, len(ssp['primary mass'])))/ssp['black hole mass'] # km/s
ssp['kick v_i'] = 0.07*numpy.sqrt((ssp['primary mass'] + ssp['companion mass'])/ssp['pre kick semi major axis']) # km/s
ssp['kick chi'] = (ssp['primary mass']+ssp['initial companion mass'])/(ssp['black hole mass'] + ssp['terminal companion mass'])

ssp['terminal semi major axis'] = ssp['pre kick semi major axis']/(
    kick_semi_major_axis_ratio(ssp['kick v_k'], 
                               ssp['kick v_i'], 
                               ssp['kick chi'],
                               ssp['kick mu'])
)

ssp['terminal eccentricity'] = kick_eccentricity(ssp['kick v_k'],
                                                 ssp['kick v_i'],
                                                 ssp['kick chi'], 
                                                 ssp['kick mu'], 
                                                 ssp['kick phi'])

ssp['terminal period'] = 9.4e7*(ssp['terminal semi major axis']**1.5*
                                ssp['primary mass']**-0.5*
                                (0.2+ssp['initial mass ratio'])**-0.5)
ssp['companion absolute magnitude'] = stellar_mass_to_vmag(ssp['primary mass']*ssp['initial mass ratio'])
ssp['distance from earth'] = numpy.sqrt(ssp['z']**2+8000**2+ssp['r']**2-2*8000*ssp['r']*numpy.cos(ssp['angle']))
ssp['companion apparent magnitude'] = absolute_to_apparent_magnitude(ssp['companion absolute magnitude'],
                                                                     ssp['distance from earth'],
                                                                    extinction=True)
ssp['sigma pomega'] = calc_sigma_pomega(ssp['companion apparent magnitude'])
ssp['parallax'] = 1e6/ssp['distance from earth'] # in uas
ssp['sigma G'] = calc_sigma_G(ssp['companion apparent magnitude'])
ssp['statistical weight'] = numpy.clip(ssp['lifetime'],0,1)
ssp_kick_momentum=ssp
```

    /usr/local/lib/python2.7/dist-packages/ipykernel_launcher.py:49: RuntimeWarning: invalid value encountered in power


Hipparcos prediction


```python
filter_list = [
    ssp['terminal semi major axis'] > 0,
    ssp['terminal eccentricity'] <0.9,
    ssp['lifetime'] > ssp['age'],
    ssp['terminal semi major axis']>ssp['terminal roche radius'],
    ssp['terminal period']<5,
    ssp['terminal period']>0.137,
    ssp['black hole mass']>
    3.0/(1.0-100*ssp['sigma pomega']/ssp['parallax']-
         100*ssp['sigma pomega']*ssp['distance from earth']/2e11/ssp['terminal semi major axis']),
    ssp['companion apparent magnitude']<12.4,
    ssp['terminal semi major axis']>
    10*(1+ssp['initial mass ratio']/0.2)*100*ssp['sigma pomega']*ssp['distance from earth']*5e-12
]
mask = numpy.all(filter_list, axis=0)
ssp['statistical weight'][mask].sum()
```

Gaia prediction


```python
filter_list = [
    ssp['terminal semi major axis'] > 0,
    ssp['terminal eccentricity'] <0.9,
    ssp['terminal semi major axis']>ssp['terminal roche radius'],
    ssp['terminal period']<5,
    ssp['terminal period']>0.137,
    ssp['black hole mass']>
    3.0/(1.0-ssp['sigma pomega']/ssp['parallax']-ssp['sigma G']-
         ssp['sigma pomega']*ssp['distance from earth']/2e11/ssp['terminal semi major axis']),
    ssp['companion apparent magnitude']<20,
    ssp['terminal semi major axis']>
    10*(1+ssp['initial mass ratio']/0.2)*ssp['sigma pomega']*ssp['distance from earth']*5e-12
]
mask = numpy.all(filter_list, axis=0)
gaia_kick_momentum_mask = mask
ssp['statistical weight'][mask].sum()
```


```python
%matplotlib inline
import matplotlib
font = {'size':22}
matplotlib.rc('font',**font)

pylab.hist(ssp['black hole mass'][mask],10, weights=ssp['statistical weight'][mask])
pylab.xlabel('Black hole mass [$M_{\odot}$]');
pylab.ylabel('Number')
pass
```


```python
pylab.hist(ssp['distance from earth'][mask]/1000,10, weights=ssp['statistical weight'][mask])
pylab.xlabel('Distance from earth [kpc]');
pylab.ylabel('Number')
pass
```


```python
pylab.hist(ssp['companion mass'][mask],10, weights=ssp['statistical weight'][mask])
pylab.xscale('log')
pylab.xlabel(r'Companion mass [$M_{\odot}$]');
pylab.ylabel('Number')
pass
```


```python
pylab.hist(ssp['terminal semi major axis'][mask]/ssp['distance from earth'][mask]/4.8e-9,10,
          weights=ssp['statistical weight'][mask])
pylab.xscale('log')
pylab.xlabel(r'Angular separation [mas]')
pylab.ylabel('Number')
pass
```


```python
pylab.scatter(ssp['companion mass'][mask],
              ssp['terminal semi major axis'][mask]/ssp['distance from earth'][mask]/4.8e-9
             )
pylab.xscale('log')
pylab.yscale('log')
pylab.xlabel('Companion mass[$M_{\odot}$]')
pylab.ylabel('Angular separation [mas]')
pass
```

# Figures for the paper


```python
%matplotlib inline
font = {'size':22}
matplotlib.rc('font',**font)

for ssp, mask in zip([ssp_no_kick, ssp_inefficient, ssp_kick_momentum, ssp_kick_velocity],
                    [gaia_no_kick_mask, gaia_inefficient_mask, gaia_kick_momentum_mask, gaia_kick_velocity_mask]):
    pylab.hist(ssp['terminal semi major axis'][mask]/
               ssp['distance from earth'][mask]/4.8e-9,10,
               weights=ssp['statistical weight'][mask]*100,
               histtype='step')
pylab.yscale('log')
pylab.xlabel(r'Angular separation [mas]')
pylab.ylabel('Number')
pass
```


```python
%matplotlib inline
font = {'size':22}
matplotlib.rc('font',**font)

for ssp, mask in zip([ssp_no_kick, ssp_inefficient, ssp_kick_momentum, ssp_kick_velocity],
                    [gaia_no_kick_mask, gaia_inefficient_mask, gaia_kick_momentum_mask, gaia_kick_velocity_mask]):
    pylab.hist(ssp['black hole mass'][mask],
               weights=ssp['statistical weight'][mask]*100,
               histtype='step',
              linewidth=4)
pylab.yscale('log')
pylab.xlabel(r'Black hole mass [$M_{\odot}$]')
pylab.ylabel('Number')
pylab.legend(['fid','al01','nkm','nkv'],bbox_to_anchor=(1,1))
pylab.tight_layout()
pass
```


```python
%matplotlib inline
font = {'size':22}
matplotlib.rc('font',**font)

for ssp, mask in zip([ssp_no_kick, ssp_inefficient, ssp_kick_momentum, ssp_kick_velocity],
                    [gaia_no_kick_mask, gaia_inefficient_mask, gaia_kick_momentum_mask, gaia_kick_velocity_mask]):
    pylab.hist(ssp['terminal companion mass'][mask],
               weights=ssp['statistical weight'][mask]*100,
               histtype='step',
              linewidth=4)
pylab.yscale('log')
pylab.xscale('log')
pylab.xlabel(r'Companion mass [$M_{\odot}$]')
pylab.legend(['fid','al01','nkm','nkv'],loc='best',bbox_to_anchor=(1,1))
pylab.ylabel('Number')
pass
```


```python
%matplotlib inline
font = {'size':22}
matplotlib.rc('font',**font)

for ssp, mask in zip([ssp_no_kick, ssp_inefficient, ssp_kick_momentum, ssp_kick_velocity],
                    [gaia_no_kick_mask, gaia_inefficient_mask, gaia_kick_momentum_mask, gaia_kick_velocity_mask]):
    pylab.hist(ssp['terminal period'][mask],
               weights=ssp['statistical weight'][mask]*100,
               histtype='step',
              linewidth=4)
pylab.yscale('log')
pylab.xlabel(r'Period [years]')
pylab.ylabel('Number')
pylab.legend(['fid','al01','nkm','nkv'],loc='best',bbox_to_anchor=(1,1))
pass
```


```python
%matplotlib inline
font = {'size':22}
matplotlib.rc('font',**font)

for ssp, mask in zip([ssp_no_kick, ssp_inefficient, ssp_kick_momentum, ssp_kick_velocity],
                    [gaia_no_kick_mask, gaia_inefficient_mask, gaia_kick_momentum_mask, gaia_kick_velocity_mask]):
    pylab.hist(ssp['distance from earth'][mask]/1000,
               weights=ssp['statistical weight'][mask]*100,
               histtype='step',
              linewidth=4)
pylab.yscale('log')
pylab.xlabel(r'Distance from earth [kpc]')
pylab.ylabel('Number')
pylab.legend(['fid','al01','nkm','nkv'],loc='best',bbox_to_anchor=(1,1))
pass
```

# Machine Learning Playground


```python
X = numpy.vstack((
    numpy.concatenate((ssp_no_kick['black hole mass'][gaia_no_kick_mask],
                       ssp_kick_momentum['black hole mass'][gaia_kick_momentum_mask],
                       ssp_kick_velocity['black hole mass'][gaia_kick_velocity_mask])),
        numpy.concatenate((ssp_no_kick['companion mass'][gaia_no_kick_mask],
                       ssp_kick_momentum['companion mass'][gaia_kick_momentum_mask],
                       ssp_kick_velocity['companion mass'][gaia_kick_velocity_mask])),
        numpy.concatenate((ssp_no_kick['terminal semi major axis'][gaia_no_kick_mask],
                       ssp_kick_momentum['terminal semi major axis'][gaia_kick_momentum_mask],
                       ssp_kick_velocity['terminal semi major axis'][gaia_kick_velocity_mask])),
    numpy.concatenate((ssp_no_kick['distance from earth'][gaia_no_kick_mask],
                       ssp_kick_momentum['distance from earth'][gaia_kick_momentum_mask],
                       ssp_kick_velocity['distance from earth'][gaia_kick_velocity_mask])))).T
    #ssp_no_kick['terminal semi major axis'][gaia_no_kick_mask],
    #ssp_no_kick['distance from earth'][gaia_no_kick_mask])).T
X_embedded = TSNE(n_components=2).fit_transform(X)
n1 = len(ssp_no_kick['black hole mass'][gaia_no_kick_mask])
n2 = len(ssp_kick_momentum['black hole mass'][gaia_kick_momentum_mask])
n3 = len(ssp_kick_velocity['black hole mass'][gaia_kick_velocity_mask])
%matplotlib inline
pylab.scatter(X_embedded.T[0][:n1],
              X_embedded.T[1][:n1])
pylab.scatter(X_embedded.T[0][n1:n1+n2],
              X_embedded.T[1][n1:n1+n2])
pylab.scatter(X_embedded.T[0][n1+n2:n1+n2+n3],
              X_embedded.T[1][n1+n2:n1+n2+n3])
```

All models occupy the same region in a reduced dimenion space. This mean they are indistinguishable.


```python

```

6e561beffabb388f846329749cdeaabe658f964a
ipynb
Jupyter Notebook
stellar_binary_monte_carlo.ipynb
bolverk/gaia_black_hole_binaries
966c35cf122e46923fcb8c0b4d9a6bbb66174751
MIT
2018-06-30T01:15:10.000Z
2018-06-30T01:15:10.000Z
stellar_binary_monte_carlo.ipynb
bolverk/gaia_black_hole_binaries
966c35cf122e46923fcb8c0b4d9a6bbb66174751
MIT
2017-11-09T18:46:04.000Z
2018-02-16T16:53:52.000Z
stellar_binary_monte_carlo.ipynb
bolverk/gaia_black_hole_binaries
966c35cf122e46923fcb8c0b4d9a6bbb66174751
MIT
Qwen/Qwen-72B
1. YES
2. YES