Using COSMIC to evolve binaries
COSMIC can evolve binaries for several different use cases. Below you’ll find examples to run a single binary system, multiple binary systems or a grid of binaries.
single binary
Below is the process to initialize and evolve a binary that could have formed a GW150914-like binary. First, import the modules in COSMIC that initialize and evolve the binary.
In [1]: from cosmic.sample.initialbinarytable import InitialBinaryTable
In [2]: from cosmic.evolve import Evolve
To initialize a single binary, populate the InitialBinaries method in the InitialBinaryTable class. Each initialized binary requires the following parameters:
m1 : ZAMS mass of the primary star in
m2 : ZAMS mass of the secondary star in
porb : initial orbital period in days
ecc : initial eccentricity
tphysf : total evolution time of the binary in Myr
kstar1 : initial primary stellar type, following the BSE convention
kstar2 : initial secondary stellar type, following the BSE convention
metallicity : metallicity of the population (e.g. )
In [3]: single_binary = InitialBinaryTable.InitialBinaries(m1=85.543645, m2=84.99784, porb=446.795757, ecc=0.448872, tphysf=13700.0, kstar1=1, kstar2=1, metallicity=0.002)
In [4]: print(single_binary)
kstar_1 kstar_2 mass_1 mass_2 porb ecc metallicity tphysf mass0_1 mass0_2 rad_1 rad_2 lum_1 ... B_2 bacc_1 bacc_2 tacc_1 tacc_2 epoch_1 epoch_2 tms_1 tms_2 bhspin_1 bhspin_2 tphys binfrac
0 1.0 1.0 85.543645 84.99784 446.795757 0.448872 0.002 13700.0 85.543645 84.99784 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
[1 rows x 38 columns]
The flags for the various binary evolution prescriptions used in BSE also need to be set. Each flag is saved in the BSEDict dictionary. Note that the BSEDict only needs to be specified the first time a binary is evolved with COSMIC or if you need to change the binary evolution prescriptions.
If you are unfamiliar with these prescriptions, it is highly advised to either run the defaults from the COSMIC install which are consistent with Breivik+2020
In [5]: BSEDict = {'xi': 1.0, 'bhflag': 1, 'neta': 0.5, 'windflag': 3, 'wdflag': 1, 'alpha1': 1.0, 'pts1': 0.001, 'pts3': 0.02, 'pts2': 0.01, 'epsnov': 0.001, 'hewind': 0.5, 'ck': 1000, 'bwind': 0.0, 'lambdaf': 0.0, 'mxns': 3.0, 'beta': -1.0, 'tflag': 1, 'acc2': 1.5, 'grflag' : 1, 'remnantflag': 4, 'ceflag': 0, 'eddfac': 1.0, 'ifflag': 0, 'bconst': 3000, 'sigma': 265.0, 'gamma': -2.0, 'pisn': 45.0, 'natal_kick_array' : [[-100.0,-100.0,-100.0,-100.0,0.0], [-100.0,-100.0,-100.0,-100.0,0.0]], 'bhsigmafrac' : 1.0, 'polar_kick_angle' : 90, 'qcrit_array' : [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], 'cekickflag' : 2, 'cehestarflag' : 0, 'cemergeflag' : 0, 'ecsn' : 2.25, 'ecsn_mlow' : 1.6, 'aic' : 1, 'ussn' : 0, 'sigmadiv' :-20.0, 'qcflag' : 1, 'eddlimflag' : 0, 'fprimc_array' : [2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0], 'bhspinflag' : 0, 'bhspinmag' : 0.0, 'rejuv_fac' : 1.0, 'rejuvflag' : 0, 'htpmb' : 1, 'ST_cr' : 1, 'ST_tide' : 1, 'bdecayfac' : 1, 'rembar_massloss' : 0.5, 'kickflag' : 0, 'zsun' : 0.014, 'bhms_coll_flag' : 0, 'don_lim' : -1, 'acc_lim' : -1}
Once the binary is initialized and the BSE model is set, the system is evolved with the the Evolve class, which calls the evolv2.f subroutine in the BSE source code.
In [6]: bpp, bcm, initC, kick_info = Evolve.evolve(initialbinarytable=single_binary, BSEDict=BSEDict)
For every evolved binary system, BSE generates two arrays, which are stored as pandas DataFrames in COSMIC:
bpp - contains binary parameters at important stages in the binary’s evolution, including stellar evolutionary phase changes or mass transfer episodes.
bcm - contains several binary parameters at user specified time steps during the binary’s evolution. The default setting in COSMIC is to output the final stage of the binary at the evolution time specified by the user.
You can see the different parameters included in each DataFrame using the columns attribute of the DataFrame:
In [7]: print(bpp.columns)
Index(['tphys', 'mass_1', 'mass_2', 'kstar_1', 'kstar_2', 'sep', 'porb', 'ecc',
'RRLO_1', 'RRLO_2', 'evol_type', 'aj_1', 'aj_2', 'tms_1', 'tms_2',
'massc_1', 'massc_2', 'rad_1', 'rad_2', 'mass0_1', 'mass0_2', 'lum_1',
'lum_2', 'teff_1', 'teff_2', 'radc_1', 'radc_2', 'menv_1', 'menv_2',
'renv_1', 'renv_2', 'omega_spin_1', 'omega_spin_2', 'B_1', 'B_2',
'bacc_1', 'bacc_2', 'tacc_1', 'tacc_2', 'epoch_1', 'epoch_2',
'bhspin_1', 'bhspin_2', 'bin_num'],
dtype='object')
In [8]: print(bcm.columns)
Index(['tphys', 'kstar_1', 'mass0_1', 'mass_1', 'lum_1', 'rad_1', 'teff_1',
'massc_1', 'radc_1', 'menv_1', 'renv_1', 'epoch_1', 'omega_spin_1',
'deltam_1', 'RRLO_1', 'kstar_2', 'mass0_2', 'mass_2', 'lum_2', 'rad_2',
'teff_2', 'massc_2', 'radc_2', 'menv_2', 'renv_2', 'epoch_2',
'omega_spin_2', 'deltam_2', 'RRLO_2', 'porb', 'sep', 'ecc', 'B_1',
'B_2', 'SN_1', 'SN_2', 'bin_state', 'merger_type', 'bin_num'],
dtype='object')
The units are broadly consistent with BSE and are described in Describing the output of COSMIC/BSE: Columns names/Values/Units.
The evol_type column in bpp indicates the evolutionary change that occurred for each line. The meaning of each number is described here, Evolve Type.
Each of the parameters in bpp or bcm can be accessed in the usual way for DataFrames:
In [9]: bpp.mass_1
Out[9]:
0 85.543645
0 72.720332
0 72.589218
0 71.959504
0 32.352601
0 32.176152
0 25.488585
0 24.988585
0 24.988590
0 24.989628
0 24.990676
0 24.990676
0 24.990676
0 24.990676
0 24.990687
0 24.990687
0 24.990687
Name: mass_1, dtype: float64
In [10]: bpp = bpp[['mass_1', 'mass_2', 'kstar_1', 'kstar_2', 'sep', 'evol_type']]
You can use the utils.convert_kstar_evol_type
function to convert the
kstar_1
, kstar_2
, and evol_type
columns from integers to strings
that describe each int:
In [11]: from cosmic.utils import convert_kstar_evol_type
In [12]: convert_kstar_evol_type(bpp)
Out[12]:
mass_1 mass_2 kstar_1 kstar_2 sep evol_type
0 85.543645 84.997840 MS, > 0.7 M⊙ MS, > 0.7 M⊙ 1363.435508 initial state
0 72.720332 72.376321 Hertzsprung Gap MS, > 0.7 M⊙ 1602.531512 kstar change
0 72.589218 72.377845 Hertzsprung Gap MS, > 0.7 M⊙ 883.827396 begin Roche lobe overflow
0 71.959504 72.806393 Core Helium Burning MS, > 0.7 M⊙ 882.511619 kstar change
0 32.352601 110.573279 Core Helium Burning MS, > 0.7 M⊙ 1614.251471 end Roche lobe overlow
0 32.176152 110.618802 Naked Helium Star MS MS, > 0.7 M⊙ 1614.063847 kstar change
0 25.488585 106.891756 Naked Helium Star MS MS, > 0.7 M⊙ 1741.069066 supernova of primary
0 24.988585 106.891756 Black Hole MS, > 0.7 M⊙ 1747.695130 kstar change
0 24.988590 88.885767 Black Hole Hertzsprung Gap 2023.280971 kstar change
0 24.989628 88.681512 Black Hole Hertzsprung Gap 2018.828817 begin Roche lobe overflow
0 24.990676 88.469775 Black Hole Core Helium Burning 2017.805217 kstar change
0 24.990676 88.469775 Black Hole Core Helium Burning 2017.805217 begin common envelope
0 24.990676 41.981621 Black Hole Naked Helium Star MS 320.660165 end common envelope
0 24.990676 41.981621 Black Hole Naked Helium Star MS 320.660165 end Roche lobe overlow
0 24.990687 31.554344 Black Hole Naked Helium Star MS 379.806762 supernova of secondary
0 24.990687 31.054344 Black Hole Black Hole 383.225670 kstar change
0 24.990687 31.054344 Black Hole Black Hole 382.989460 max evolution time
Note that utils.convert_kstar_evol_type
is only applicable to the bpp
array.
You can also use the built in plotting function to see how the system evolves:
In [13]: from cosmic.plotting import evolve_and_plot
In [14]: single_binary = InitialBinaryTable.InitialBinaries(m1=85.543645, m2=84.99784, porb=446.795757, ecc=0.448872, tphysf=13700.0, kstar1=1, kstar2=1, metallicity=0.002)
In [15]: BSEDict = {'xi': 1.0, 'bhflag': 1, 'neta': 0.5, 'windflag': 3, 'wdflag': 1, 'alpha1': 1.0, 'pts1': 0.001, 'pts3': 0.02, 'pts2': 0.01, 'epsnov': 0.001, 'hewind': 0.5, 'ck': 1000, 'bwind': 0.0, 'lambdaf': 0.0, 'mxns': 3.0, 'beta': -1.0, 'tflag': 1, 'acc2': 1.5, 'grflag' : 1, 'remnantflag': 4, 'ceflag': 0, 'eddfac': 1.0, 'ifflag': 0, 'bconst': 3000, 'sigma': 265.0, 'gamma': -2.0, 'pisn': 45.0, 'natal_kick_array' : [[-100.0,-100.0,-100.0,-100.0,0.0], [-100.0,-100.0,-100.0,-100.0,0.0]], 'bhsigmafrac' : 1.0, 'polar_kick_angle' : 90, 'qcrit_array' : [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], 'cekickflag' : 2, 'cehestarflag' : 0, 'cemergeflag' : 0, 'ecsn' : 2.25, 'ecsn_mlow' : 1.6, 'aic' : 1, 'ussn' : 0, 'sigmadiv' :-20.0, 'qcflag' : 1, 'eddlimflag' : 0, 'fprimc_array' : [2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0], 'bhspinflag' : 0, 'bhspinmag' : 0.0, 'rejuv_fac' : 1.0, 'rejuvflag' : 0, 'htpmb' : 1, 'ST_cr' : 1, 'ST_tide' : 1, 'bdecayfac' : 1, 'rembar_massloss' : 0.5, 'kickflag' : 0, 'zsun' : 0.014, 'bhms_coll_flag' : 0, 'don_lim' : -1, 'acc_lim' : -1}
In [16]: fig = evolve_and_plot(single_binary, t_min=None, t_max=None, BSEDict=BSEDict, sys_obs={})
(Source code, png, hires.png, pdf)
In this case, all the action happens in the first few Myr, so let’s specify a t_max:
In [17]: fig = evolve_and_plot(initC, t_min=None, t_max=6.0, BSEDict={}, sys_obs={})
(Source code, png, hires.png, pdf)
multiple binaries
Multiple systems can also be initialized and evolved; below is an example for systems that could form GW150914 and GW170817 - like binaries.
In [18]: binary_set = InitialBinaryTable.InitialBinaries(m1=[85.543645, 11.171469], m2=[84.99784, 6.67305], porb=[446.795757, 170.758343], ecc=[0.448872, 0.370], tphysf=[13700.0, 13700.0], kstar1=[1, 1], kstar2=[1, 1], metallicity=[0.002, 0.02])
In [19]: print(binary_set)
kstar_1 kstar_2 mass_1 mass_2 porb ecc metallicity tphysf mass0_1 mass0_2 rad_1 rad_2 lum_1 ... B_2 bacc_1 bacc_2 tacc_1 tacc_2 epoch_1 epoch_2 tms_1 tms_2 bhspin_1 bhspin_2 tphys binfrac
0 1.0 1.0 85.543645 84.99784 446.795757 0.448872 0.002 13700.0 85.543645 84.99784 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
1 1.0 1.0 11.171469 6.67305 170.758343 0.370000 0.020 13700.0 11.171469 6.67305 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
[2 rows x 38 columns]
In [20]: import numpy as np
In [21]: np.random.seed(5)
In [22]: bpp, bcm, initC, kick_info = Evolve.evolve(initialbinarytable=binary_set, BSEDict=BSEDict)
Note that the BSEDict did not be reinitialized since the BSE model did not change.
As before, bpp, bcm, and initC are returned as pandas DataFrames which assign an index to each binary system we evolve. We can access each binary as follows:
In [23]: print(bpp.loc[0])
tphys mass_1 mass_2 kstar_1 kstar_2 sep porb ecc RRLO_1 RRLO_2 evol_type ... B_1 B_2 bacc_1 bacc_2 tacc_1 tacc_2 epoch_1 epoch_2 bhspin_1 bhspin_2 bin_num
0 0.000000 85.543645 84.997840 1.0 1.0 1363.435508 446.795757 0.448872 3.878791e-02 3.877174e-02 1.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.000000 0.000000 0.0 0.0 0
0 3.717070 72.719783 72.376037 2.0 1.0 1602.540712 617.246046 0.448872 1.074695e-01 1.299162e-01 2.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 -0.130507 -0.129952 0.0 0.0 0
0 3.718368 72.588666 72.377556 2.0 1.0 883.832613 252.926180 0.000000 1.000132e+00 1.268762e-01 3.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 -0.130507 -0.129922 0.0 0.0 0
0 3.720007 71.958903 72.806115 4.0 1.0 882.510342 252.534112 0.000000 1.686913e+01 1.504802e-01 2.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 -0.130507 -0.101292 0.0 0.0 0
0 3.739911 32.346795 110.583971 4.0 1.0 1616.153256 629.844167 0.000000 3.282858e-01 2.753324e-02 4.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 -0.130507 1.466049 0.0 0.0 0
0 3.741053 32.175552 110.600277 7.0 1.0 1617.101339 630.740428 0.000000 4.252227e-03 2.751049e-02 2.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 3.722851 1.466122 0.0 0.0 0
0 4.071364 25.486398 106.874794 7.0 1.0 1744.364154 733.918026 0.000000 3.606448e-03 3.002358e-02 15.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 3.698581 1.447327 0.0 0.0 0
0 4.071364 24.986398 106.874794 14.0 1.0 1751.003726 739.510366 0.003792 2.311299e-07 2.993589e-02 2.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 4.071364 1.447327 0.0 0.0 0
0 4.894907 24.986403 88.691971 14.0 2.0 2030.288753 994.420218 0.003792 1.895385e-07 4.575119e-02 2.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 4.071364 1.289560 0.0 0.0 0
0 4.896430 24.987444 88.487242 14.0 2.0 2025.833989 992.038312 0.000000 1.891237e-07 1.001885e+00 3.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 4.071364 1.289560 0.0 0.0 0
0 4.897383 24.988491 88.276317 14.0 4.0 2024.854520 992.236963 0.000000 1.891001e-07 6.931196e+00 2.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 4.071364 1.289560 0.0 0.0 0
0 4.897383 24.988491 88.276317 14.0 4.0 2024.854520 992.236963 0.000000 1.891001e-07 6.931196e+00 7.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 4.071364 1.289560 0.0 0.0 0
0 4.897383 24.988491 41.858050 14.0 7.0 320.766753 81.436052 0.000000 1.891001e-07 6.931196e+00 8.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 4.071364 1.289560 0.0 0.0 0
0 4.897383 24.988491 41.858050 14.0 7.0 320.766753 81.436052 0.000000 9.852501e-07 1.645015e-02 4.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 4.071364 4.897383 0.0 0.0 0
0 5.206247 24.988501 31.476750 14.0 7.0 379.755264 114.139832 0.000000 7.807651e-07 1.258480e-02 16.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 4.071364 4.872584 0.0 0.0 0
0 5.206247 24.988501 30.976750 14.0 14.0 383.178626 116.202332 0.008934 7.739647e-07 8.697561e-07 2.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 4.071364 5.206247 0.0 0.0 0
0 205.190000 24.988501 30.976750 14.0 14.0 382.959826 116.102817 0.008934 7.744068e-07 8.702530e-07 10.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 4.071364 5.206247 0.0 0.0 0
[17 rows x 44 columns]
In [24]: print(bcm.loc[0])
tphys kstar_1 mass0_1 mass_1 lum_1 rad_1 teff_1 massc_1 radc_1 menv_1 renv_1 ... RRLO_2 porb sep ecc B_1 B_2 SN_1 SN_2 bin_state merger_type bin_num
0 0.00 1.0 85.543645 85.543645 9.416137e+05 11.060285 54306.854457 0.000000 0.000000 1.000000e-10 1.000000e-10 ... 3.877174e-02 446.795757 1363.435508 0.448872 0.0 0.0 0.0 0.0 0 -001 0
0 0.01 1.0 85.527367 85.527367 9.423653e+05 11.069603 54294.820398 0.000000 0.000000 1.000000e-10 1.000000e-10 ... 3.879707e-02 446.965185 1363.693996 0.448872 0.0 0.0 0.0 0.0 0 -001 0
0 0.02 1.0 85.511049 85.511049 9.431226e+05 11.078928 54282.866439 0.000000 0.000000 1.000000e-10 1.000000e-10 ... 3.882240e-02 447.135124 1363.953214 0.448872 0.0 0.0 0.0 0.0 0 -001 0
0 0.03 1.0 85.494691 85.494691 9.438858e+05 11.088262 54270.989379 0.000000 0.000000 1.000000e-10 1.000000e-10 ... 3.884773e-02 447.305577 1364.213166 0.448872 0.0 0.0 0.0 0.0 0 -001 0
0 0.04 1.0 85.478292 85.478292 9.446548e+05 11.097605 54259.186024 0.000000 0.000000 1.000000e-10 1.000000e-10 ... 3.887306e-02 447.476547 1364.473858 0.448872 0.0 0.0 0.0 0.0 0 -001 0
.. ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
0 205.15 14.0 25.486398 24.988501 1.000000e-10 0.000106 1781.215957 24.988501 0.000106 1.000000e-10 1.000000e-10 ... 8.702530e-07 116.102817 382.959826 0.008934 0.0 0.0 1.0 1.0 0 -001 0
0 205.16 14.0 25.486398 24.988501 1.000000e-10 0.000106 1781.215957 24.988501 0.000106 1.000000e-10 1.000000e-10 ... 8.702530e-07 116.102817 382.959826 0.008934 0.0 0.0 1.0 1.0 0 -001 0
0 205.17 14.0 25.486398 24.988501 1.000000e-10 0.000106 1781.215957 24.988501 0.000106 1.000000e-10 1.000000e-10 ... 8.702530e-07 116.102817 382.959826 0.008934 0.0 0.0 1.0 1.0 0 -001 0
0 205.18 14.0 25.486398 24.988501 1.000000e-10 0.000106 1781.215957 24.988501 0.000106 1.000000e-10 1.000000e-10 ... 8.702530e-07 116.102817 382.959826 0.008934 0.0 0.0 1.0 1.0 0 -001 0
0 205.19 14.0 25.486398 24.988501 1.000000e-10 0.000106 1781.215957 24.988501 0.000106 1.000000e-10 1.000000e-10 ... 8.702530e-07 116.102817 382.959826 0.008934 0.0 0.0 1.0 1.0 0 -001 0
[20520 rows x 39 columns]
In [25]: print(initC.loc[0])
kstar_1 1.000000
kstar_2 1.000000
mass_1 85.543645
mass_2 84.997840
porb 446.795757
...
fprimc_11 0.095238
fprimc_12 0.095238
fprimc_13 0.095238
fprimc_14 0.095238
fprimc_15 0.095238
Name: 0, Length: 136, dtype: float64
In [26]: print(bpp.loc[1])
tphys mass_1 mass_2 kstar_1 kstar_2 sep porb ecc RRLO_1 RRLO_2 evol_type ... B_1 B_2 bacc_1 bacc_2 tacc_1 tacc_2 epoch_1 epoch_2 bhspin_1 bhspin_2 bin_num
1 0.000000 11.171469 6.673050 1.0 1.0 338.356712 170.758343 0.370000 0.049210 0.045915 1.0 ... 0.0 0.000000e+00 0.0 0.0 0.0 0.0 0.000000 0.000000 0.0 0.0 1
1 19.427278 10.767909 6.664658 2.0 1.0 346.350714 178.922774 0.369999 0.115026 0.053609 2.0 ... 0.0 0.000000e+00 0.0 0.0 0.0 0.0 -0.969072 -0.029132 0.0 0.0 1
1 19.461438 10.765777 6.664702 2.0 1.0 218.223971 89.489313 0.000000 1.000479 0.053622 3.0 ... 0.0 0.000000e+00 0.0 0.0 0.0 0.0 -0.976313 -0.028841 0.0 0.0 1
1 19.476912 10.544417 6.884988 3.0 1.0 211.992707 85.686473 0.000000 2.771636 0.055156 2.0 ... 0.0 0.000000e+00 0.0 0.0 0.0 0.0 -1.094470 1.906517 0.0 0.0 1
1 19.476912 10.544417 6.884988 3.0 1.0 211.992707 85.686473 0.000000 2.771636 0.055156 7.0 ... 0.0 0.000000e+00 0.0 0.0 0.0 0.0 -1.094470 1.906517 0.0 0.0 1
1 19.476912 2.415442 6.884988 7.0 1.0 20.150243 3.437487 0.000000 2.771636 0.055156 8.0 ... 0.0 0.000000e+00 0.0 0.0 0.0 0.0 -1.094470 1.906517 0.0 0.0 1
1 19.476912 2.415442 6.884988 7.0 1.0 20.150243 3.437487 0.000000 0.064435 0.421709 4.0 ... 0.0 0.000000e+00 0.0 0.0 0.0 0.0 19.476912 1.906517 0.0 0.0 1
1 22.952681 2.226294 6.887974 8.0 1.0 20.614122 3.593011 0.000000 0.060956 0.424422 2.0 ... 0.0 0.000000e+00 0.0 0.0 0.0 0.0 19.165010 1.925371 0.0 0.0 1
1 23.231406 2.177806 6.896735 8.0 1.0 20.656641 3.612014 0.000000 1.000105 0.423629 3.0 ... 0.0 0.000000e+00 0.0 0.0 0.0 0.0 19.165010 1.984764 0.0 0.0 1
1 23.254641 1.889676 7.181936 9.0 1.0 22.458960 4.095570 0.000000 25.240301 0.382091 2.0 ... 0.0 0.000000e+00 0.0 0.0 0.0 0.0 19.165010 4.538299 0.0 0.0 1
1 23.256935 1.371085 7.700471 9.0 1.0 31.223370 6.713530 0.000000 0.274471 0.263521 4.0 ... 0.0 0.000000e+00 0.0 0.0 0.0 0.0 19.165010 8.211818 0.0 0.0 1
1 23.258694 1.370344 7.700740 12.0 1.0 31.223233 6.713660 0.000000 0.000146 0.263511 2.0 ... 0.0 0.000000e+00 0.0 0.0 0.0 0.0 23.258694 8.212927 0.0 0.0 1
1 47.641575 1.370365 7.628335 12.0 2.0 32.650896 7.208168 0.000000 0.000139 0.483936 2.0 ... 0.0 0.000000e+00 0.0 0.0 0.0 0.0 23.258694 7.614577 0.0 0.0 1
1 47.673758 1.370368 7.628122 12.0 2.0 32.703941 7.225825 0.000000 0.000139 1.000649 3.0 ... 0.0 0.000000e+00 0.0 0.0 0.0 0.0 23.258694 7.612219 0.0 0.0 1
1 47.673758 1.370368 7.628122 12.0 2.0 32.703941 7.225825 0.000000 0.000139 1.000649 7.0 ... 0.0 0.000000e+00 0.0 0.0 0.0 0.0 23.258694 7.612219 0.0 0.0 1
1 47.673758 1.370344 8.188112 15.0 5.0 0.000000 0.000000 0.000000 0.000139 1.000649 8.0 ... 0.0 0.000000e+00 0.0 0.0 0.0 0.0 23.258694 7.612219 0.0 0.0 1
1 47.785623 0.000000 7.998227 15.0 5.0 0.000000 0.000000 -1.000000 -1.000000 0.000100 16.0 ... 0.0 0.000000e+00 0.0 0.0 0.0 0.0 23.258694 21.502403 0.0 0.0 1
1 47.785623 0.000000 1.277584 15.0 13.0 0.000000 0.000000 -1.000000 -1.000000 0.000100 2.0 ... 0.0 0.000000e+00 0.0 0.0 0.0 0.0 23.258694 47.785623 0.0 0.0 1
1 247.770000 0.000000 1.277584 15.0 13.0 0.000000 0.000000 -1.000000 -1.000000 0.000100 10.0 ... 0.0 5.688661e+12 0.0 0.0 0.0 0.0 23.258694 47.785623 0.0 0.0 1
[19 rows x 44 columns]
The plotting function can also take in multiple binaries. Let’s plot both the GW150914-like progenitor evolution and the GW170817-like progenitor evolutions. For the GW170817-like progenitor, we expect most of the evolution to take place in the first ~60 Myr.
In [27]: fig = evolve_and_plot(binary_set, t_min=None, t_max=[6.0, 60.0], BSEDict=BSEDict, sys_obs={})
grid of binaries
Sometimes it is helpful to run a grid of initial binaries to explore how changing a single paramter affects the evolved binary. Here we evolve the same system that produces a GW150914-like binary, but run over several initial orbital periods spaced evenly in log space.
In [28]: n_grid = 10
In [29]: binary_grid = InitialBinaryTable.InitialBinaries(m1=np.ones(n_grid)*100.0, m2=np.ones(n_grid)*85.0, porb=np.logspace(3,5,n_grid), ecc=np.ones(n_grid)*0.65, tphysf=np.ones(n_grid)*13700.0, kstar1=np.ones(n_grid), kstar2=np.ones(n_grid), metallicity=np.ones(n_grid)*0.005)
In [30]: BSEDict = {'xi': 1.0, 'bhflag': 1, 'neta': 0.5, 'windflag': 3, 'wdflag': 1, 'alpha1': 1.0, 'pts1': 0.001, 'pts3': 0.02, 'pts2': 0.01, 'epsnov': 0.001, 'hewind': 0.5, 'ck': 1000, 'bwind': 0.0, 'lambdaf': 0.0, 'mxns': 3.0, 'beta': -1.0, 'tflag': 1, 'acc2': 1.5, 'grflag' : 1, 'remnantflag': 4, 'ceflag': 0, 'eddfac': 1.0, 'ifflag': 0, 'bconst': 3000, 'sigma': 265.0, 'gamma': -2.0, 'pisn': 45.0, 'natal_kick_array' : [[-100.0,-100.0,-100.0,-100.0,0.0], [-100.0,-100.0,-100.0,-100.0,0.0]], 'bhsigmafrac' : 1.0, 'polar_kick_angle' : 90, 'qcrit_array' : [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], 'cekickflag' : 2, 'cehestarflag' : 0, 'cemergeflag' : 0, 'ecsn' : 2.25, 'ecsn_mlow' : 1.6, 'aic' : 1, 'ussn' : 0, 'sigmadiv' :-20.0, 'qcflag' : 1, 'eddlimflag' : 0, 'fprimc_array' : [2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0], 'bhspinflag' : 0, 'bhspinmag' : 0.0, 'rejuv_fac' : 1.0, 'rejuvflag' : 0, 'htpmb' : 1, 'ST_cr' : 1, 'ST_tide' : 1, 'bdecayfac' : 1, 'rembar_massloss' : 0.5, 'kickflag' : 0, 'zsun' : 0.014, 'bhms_coll_flag' : 0, 'don_lim' : -1, 'acc_lim' : -1}
In [31]: print(binary_grid)
kstar_1 kstar_2 mass_1 mass_2 porb ecc metallicity tphysf mass0_1 mass0_2 rad_1 rad_2 lum_1 lum_2 ... B_1 B_2 bacc_1 bacc_2 tacc_1 tacc_2 epoch_1 epoch_2 tms_1 tms_2 bhspin_1 bhspin_2 tphys binfrac
0 1.0 1.0 100.0 85.0 1000.000000 0.65 0.005 13700.0 100.0 85.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
1 1.0 1.0 100.0 85.0 1668.100537 0.65 0.005 13700.0 100.0 85.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
2 1.0 1.0 100.0 85.0 2782.559402 0.65 0.005 13700.0 100.0 85.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
3 1.0 1.0 100.0 85.0 4641.588834 0.65 0.005 13700.0 100.0 85.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
4 1.0 1.0 100.0 85.0 7742.636827 0.65 0.005 13700.0 100.0 85.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
5 1.0 1.0 100.0 85.0 12915.496650 0.65 0.005 13700.0 100.0 85.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
6 1.0 1.0 100.0 85.0 21544.346900 0.65 0.005 13700.0 100.0 85.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
7 1.0 1.0 100.0 85.0 35938.136638 0.65 0.005 13700.0 100.0 85.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
8 1.0 1.0 100.0 85.0 59948.425032 0.65 0.005 13700.0 100.0 85.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
9 1.0 1.0 100.0 85.0 100000.000000 0.65 0.005 13700.0 100.0 85.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
[10 rows x 38 columns]
In [32]: bpp, bcm, initC, kick_info = Evolve.evolve(initialbinarytable=binary_grid, BSEDict=BSEDict)
In [33]: print(bpp)
tphys mass_1 mass_2 kstar_1 kstar_2 sep porb ecc RRLO_1 RRLO_2 evol_type ... B_1 B_2 bacc_1 bacc_2 tacc_1 tacc_2 epoch_1 epoch_2 bhspin_1 bhspin_2 bin_num
0 0.000000 100.000000 85.000000 1.0 1.0 2397.042879 1.000000e+03 0.650000 4.214342e-02 4.122696e-02 1.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.000000 0.000000 0.0 0.0 0
0 3.561836 74.433283 67.900444 2.0 1.0 3115.630394 1.689418e+03 0.650000 1.241381e-01 1.262321e-01 2.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 -0.192142 -0.163778 0.0 0.0 0
0 3.563725 74.195874 67.901480 2.0 1.0 1092.094610 3.508888e+02 0.000000 1.001714e+00 1.263648e-01 3.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 -0.192142 -0.163756 0.0 0.0 0
0 3.564957 73.978352 67.965270 4.0 1.0 1091.881014 3.509758e+02 0.000000 3.925547e+00 1.263083e-01 2.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 -0.192142 -0.160638 0.0 0.0 0
0 3.656147 34.358643 98.037335 4.0 1.0 1408.149978 5.322408e+02 0.000000 6.132831e-01 4.403118e-02 4.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 -0.192142 1.156534 0.0 0.0 0
.. ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
9 3.930319 28.153730 35.520384 14.0 4.0 146504.333172 8.144537e+05 0.645026 6.394120e-09 9.156138e-02 2.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 3.930319 -0.176003 0.0 0.0 9
9 3.960762 28.159125 31.484483 14.0 7.0 156308.745072 9.273958e+05 0.644909 5.824845e-09 9.079288e-05 2.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 3.930319 3.758602 0.0 0.0 9
9 4.100892 28.159125 24.732495 14.0 7.0 176275.664624 1.179417e+06 0.644909 5.824845e-09 9.079288e-05 16.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 3.930319 3.722109 0.0 0.0 9
9 4.100892 28.159125 24.232495 14.0 14.0 178591.634598 1.208462e+06 0.646513 4.824793e-09 4.446878e-09 2.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 3.930319 4.100892 0.0 0.0 9
9 13700.000000 28.159125 24.232495 14.0 14.0 178489.776222 1.207428e+06 0.646513 4.827546e-09 4.449415e-09 10.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 3.930319 4.100892 0.0 0.0 9
[140 rows x 44 columns]
In [34]: print(bcm)
tphys kstar_1 mass0_1 mass_1 lum_1 rad_1 teff_1 massc_1 radc_1 menv_1 renv_1 ... RRLO_2 porb sep ecc B_1 B_2 SN_1 SN_2 bin_state merger_type bin_num
0 0.0 1.0 100.000000 100.000000 1.240906e+06 13.897608 51907.996879 0.000000 0.000000 1.000000e-10 1.000000e-10 ... 4.122696e-02 1.000000e+03 2397.042879 0.650000 0.0 0.0 0.0 0.0 0 -001 0
0 13700.0 14.0 22.440552 21.943153 1.000000e-10 0.000093 1900.803264 21.943153 0.000093 1.000000e-10 1.000000e-10 ... 1.865774e-06 2.458504e+01 123.469509 0.011920 0.0 0.0 1.0 1.0 0 -001 0
1 0.0 1.0 100.000000 100.000000 1.240906e+06 13.897608 51907.996879 0.000000 0.000000 1.000000e-10 1.000000e-10 ... 2.931117e-02 1.668101e+03 3371.506109 0.650000 0.0 0.0 0.0 0.0 0 -001 1
1 13700.0 14.0 23.651593 23.154215 1.000000e-10 0.000098 1850.425698 23.154215 0.000098 1.000000e-10 1.000000e-10 ... 1.788960e-06 2.619483e+01 129.953548 0.011609 0.0 0.0 1.0 1.0 0 -001 1
2 0.0 1.0 100.000000 100.000000 1.240906e+06 13.897608 51907.996879 0.000000 0.000000 1.000000e-10 1.000000e-10 ... 2.083939e-02 2.782559e+03 4742.115190 0.650000 0.0 0.0 0.0 0.0 0 -001 2
2 13700.0 14.0 25.319576 25.079267 1.000000e-10 0.000106 1777.989786 25.079267 0.000106 1.000000e-10 1.000000e-10 ... 1.029859e-07 2.555575e+03 2943.446752 0.009533 0.0 0.0 1.0 1.0 0 -001 2
3 0.0 1.0 100.000000 100.000000 1.240906e+06 13.897608 51907.996879 0.000000 0.000000 1.000000e-10 1.000000e-10 ... 1.481620e-02 4.641589e+03 6669.914200 0.650000 0.0 0.0 0.0 0.0 0 -001 3
3 13700.0 14.0 27.066154 26.858148 1.000000e-10 0.000114 1718.100835 26.858148 0.000114 1.000000e-10 1.000000e-10 ... 8.746688e-08 3.389138e+03 3613.904203 0.009059 0.0 0.0 1.0 1.0 0 -001 3
4 0.0 1.0 100.000000 100.000000 1.240906e+06 13.897608 51907.996879 0.000000 0.000000 1.000000e-10 1.000000e-10 ... 1.053389e-02 7.742637e+03 9381.416025 0.650000 0.0 0.0 0.0 0.0 0 -001 4
4 13700.0 14.0 24.723822 24.223941 1.000000e-10 0.000103 1809.107130 24.223941 0.000103 1.000000e-10 1.000000e-10 ... 1.280214e-06 4.279488e+01 181.050323 0.017606 0.0 0.0 1.0 1.0 0 -001 4
5 0.0 1.0 100.000000 100.000000 1.240906e+06 13.897608 51907.996879 0.000000 0.000000 1.000000e-10 1.000000e-10 ... 7.489289e-03 1.291550e+04 13195.217208 0.650000 0.0 0.0 0.0 0.0 0 -001 5
5 13700.0 14.0 29.517261 29.070916 1.000000e-10 0.000123 1651.419177 29.070916 0.000123 1.000000e-10 1.000000e-10 ... 2.158945e-08 6.521852e+04 25850.016841 0.473375 0.0 0.0 1.0 1.0 0 -001 5
6 0.0 1.0 100.000000 100.000000 1.240906e+06 13.897608 51907.996879 0.000000 0.000000 1.000000e-10 1.000000e-10 ... 5.324667e-03 2.154435e+04 18559.432469 0.650000 0.0 0.0 0.0 0.0 0 -001 6
6 13700.0 14.0 28.990549 28.531297 1.000000e-10 0.000121 1666.962868 28.531297 0.000121 1.000000e-10 1.000000e-10 ... 1.356464e-08 1.951836e+05 53258.251279 0.604085 0.0 0.0 1.0 1.0 0 -001 6
7 0.0 1.0 100.000000 100.000000 1.240906e+06 13.897608 51907.996879 0.000000 0.000000 1.000000e-10 1.000000e-10 ... 3.785684e-03 3.593814e+04 26104.347365 0.650000 0.0 0.0 0.0 0.0 0 -001 7
7 13700.0 14.0 28.809440 28.333276 1.000000e-10 0.000120 1672.777925 28.333276 0.000120 1.000000e-10 1.000000e-10 ... 9.238461e-09 3.844055e+05 83492.641891 0.631604 0.0 0.0 1.0 1.0 0 -001 7
8 0.0 1.0 100.000000 100.000000 1.240906e+06 13.897608 51907.996879 0.000000 0.000000 1.000000e-10 1.000000e-10 ... 2.691511e-03 5.994843e+04 36716.475706 0.650000 0.0 0.0 0.0 0.0 0 -001 8
8 13700.0 14.0 28.708663 28.219714 1.000000e-10 0.000120 1676.140336 28.219714 0.000120 1.000000e-10 1.000000e-10 ... 6.306465e-09 6.873549e+05 122635.915245 0.637013 0.0 0.0 1.0 1.0 0 -001 8
9 0.0 1.0 100.000000 100.000000 1.240906e+06 13.897608 51907.996879 0.000000 0.000000 1.000000e-10 1.000000e-10 ... 1.913586e-03 1.000000e+05 51642.723315 0.650000 0.0 0.0 0.0 0.0 0 -001 9
9 13700.0 14.0 28.653730 28.159125 1.000000e-10 0.000119 1677.942604 28.159125 0.000119 1.000000e-10 1.000000e-10 ... 4.449415e-09 1.207428e+06 178489.776222 0.646513 0.0 0.0 1.0 1.0 0 -001 9
[20 rows x 39 columns]
dynamically set time resolution for bcm array
COSMIC has the ability to set time resolution of the bcm array depending on the current state of the evolution. Below we demonstrate three scenarios, setting dtp only during mass transfer, setting dtp to the same resolution for all of the evolution except for after the system merges or is disrupted, and finally an example of setting dtp only during the HMB stage of the evolution.
First, print all time steps during mass transfer
In [35]: single_binary = InitialBinaryTable.InitialBinaries(m1=7.806106, m2=5.381412, porb=2858.942021, ecc=0.601408, tphysf=13700.0, kstar1=1, kstar2=1, metallicity=0.02)
In [36]: BSEDict = {'xi': 1.0, 'bhflag': 1, 'neta': 0.5, 'windflag': 3, 'wdflag': 1, 'alpha1': 1.0, 'pts1': 0.001, 'pts3': 0.02, 'pts2': 0.01, 'epsnov': 0.001, 'hewind': 0.5, 'ck': 1000, 'bwind': 0.0, 'lambdaf': 0.0, 'mxns': 3.0, 'beta': -1.0, 'tflag': 1, 'acc2': 1.5, 'grflag' : 1, 'remnantflag': 4, 'ceflag': 0, 'eddfac': 1.0, 'ifflag': 0, 'bconst': 3000, 'sigma': 265.0, 'gamma': -2.0, 'pisn': 45.0, 'natal_kick_array' : [[-100.0,-100.0,-100.0,-100.0,0.0], [-100.0,-100.0,-100.0,-100.0,0.0]], 'bhsigmafrac' : 1.0, 'polar_kick_angle' : 90, 'qcrit_array' : [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], 'cekickflag' : 2, 'cehestarflag' : 0, 'cemergeflag' : 0, 'ecsn' : 2.25, 'ecsn_mlow' : 1.6, 'aic' : 1, 'ussn' : 0, 'sigmadiv' :-20.0, 'qcflag' : 1, 'eddlimflag' : 0, 'fprimc_array' : [2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0], 'bhspinflag' : 0, 'bhspinmag' : 0.0, 'rejuv_fac' : 1.0, 'rejuvflag' : 0, 'htpmb' : 1, 'ST_cr' : 1, 'ST_tide' : 1, 'bdecayfac' : 1, 'rembar_massloss' : 0.5, 'kickflag' : 0, 'zsun' : 0.014, 'bhms_coll_flag' : 0, 'don_lim' : -1, 'acc_lim' : -1}
In [37]: bpp, bcm, initC, kick_info = Evolve.evolve(initialbinarytable=single_binary, BSEDict=BSEDict, timestep_conditions =[['RRLO_1>=1', 'dtp=0.0'], ['RRLO_2>=1', 'dtp=0.0']])
In [38]: print(bcm[['tphys', 'kstar_1', 'kstar_2', 'mass_1', 'mass_2', 'RRLO_1', 'RRLO_2']])
tphys kstar_1 kstar_2 mass_1 mass_2 RRLO_1 RRLO_2
0 0.000000 1.0 1.0 7.806106 5.381412 0.010953 0.010459
0 43.565604 5.0 1.0 7.346069 5.396914 1.000222 0.008730
0 43.565604 8.0 1.0 2.037972 5.396914 0.520417 0.110930
0 43.579670 8.0 1.0 2.032299 5.399465 1.000951 0.110963
0 43.579674 8.0 1.0 2.032297 5.399466 1.001137 0.110963
.. ... ... ... ... ... ... ...
0 43.610746 9.0 1.0 1.323823 6.103425 1.704861 0.064712
0 43.611735 9.0 1.0 1.322435 6.104787 0.033876 0.064677
0 80.634129 12.0 2.0 1.321769 6.083292 0.000090 1.000175
0 80.634129 15.0 5.0 0.000000 5.225974 -1.000000 0.000100
0 13700.000000 15.0 13.0 0.000000 1.277584 -1.000000 0.000100
[107 rows x 7 columns]
Second, pick a certain resolution for the bcm array until the system mergers or is disrutped and then only print the final state
In [39]: bpp, bcm, initC, kick_info = Evolve.evolve(initialbinarytable=single_binary, BSEDict=BSEDict, timestep_conditions =[['binstate=0', 'dtp=1.0']])
In [40]: print(bcm[['tphys', 'kstar_1', 'kstar_2', 'mass_1', 'mass_2', 'bin_state']])
tphys kstar_1 kstar_2 mass_1 mass_2 bin_state
0 0.0 1.0 1.0 7.806106 5.381412 0
0 1.0 1.0 1.0 7.805049 5.381334 0
0 2.0 1.0 1.0 7.803970 5.381256 0
0 3.0 1.0 1.0 7.802864 5.381177 0
0 4.0 1.0 1.0 7.801731 5.381097 0
.. ... ... ... ... ... ...
0 78.0 12.0 1.0 1.321764 6.085939 0
0 79.0 12.0 1.0 1.321764 6.085078 0
0 80.0 12.0 1.0 1.321764 6.084196 0
0 81.0 15.0 13.0 0.000000 1.277584 1
0 13700.0 15.0 13.0 0.000000 1.277584 1
[83 rows x 6 columns]
Finally, we show how to print a fine resolution only during the HMXB stage of the evolution.
In [41]: single_binary = InitialBinaryTable.InitialBinaries(m1=85.543645, m2=84.99784, porb=446.795757, ecc=0.448872, tphysf=13700.0, kstar1=1, kstar2=1, metallicity=0.002)
In [42]: BSEDict = {'xi': 1.0, 'bhflag': 1, 'neta': 0.5, 'windflag': 3, 'wdflag': 1, 'alpha1': 1.0, 'pts1': 0.001, 'pts3': 0.02, 'pts2': 0.01, 'epsnov': 0.001, 'hewind': 0.5, 'ck': 1000, 'bwind': 0.0, 'lambdaf': 0.0, 'mxns': 3.0, 'beta': -1.0, 'tflag': 1, 'acc2': 1.5, 'grflag' : 1, 'remnantflag': 4, 'ceflag': 0, 'eddfac': 1.0, 'ifflag': 0, 'bconst': 3000, 'sigma': 265.0, 'gamma': -2.0, 'pisn': 45.0, 'natal_kick_array' : [[-100.0,-100.0,-100.0,-100.0,0.0], [-100.0,-100.0,-100.0,-100.0,0.0]], 'bhsigmafrac' : 1.0, 'polar_kick_angle' : 90, 'qcrit_array' : [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], 'cekickflag' : 2, 'cehestarflag' : 0, 'cemergeflag' : 0, 'ecsn' : 2.25, 'ecsn_mlow' : 1.6, 'aic' : 1, 'ussn' : 0, 'sigmadiv' :-20.0, 'qcflag' : 1, 'eddlimflag' : 0, 'fprimc_array' : [2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0], 'bhspinflag' : 0, 'bhspinmag' : 0.0, 'rejuv_fac' : 1.0, 'rejuvflag' : 0, 'htpmb' : 1, 'ST_cr' : 1, 'ST_tide' : 1, 'bdecayfac' : 1, 'rembar_massloss' : 0.5, 'kickflag' : 0, 'zsun' : 0.014, 'bhms_coll_flag' : 0, 'don_lim' : -1, 'acc_lim' : -1}
In [43]: bpp, bcm, initC, kick_info = Evolve.evolve(initialbinarytable=single_binary, BSEDict=BSEDict, timestep_conditions =[['kstar_1=14', 'kstar_2<10','dtp=0.1'], ['kstar_2=14', 'kstar_1<10','dtp=0.1']])
In [44]: print(bcm[['tphys', 'kstar_1', 'kstar_2', 'mass_1', 'mass_2', 'bin_state']])
tphys kstar_1 kstar_2 mass_1 mass_2 bin_state
0 0.000000 1.0 1.0 85.543645 84.997840 0
0 4.071374 14.0 1.0 24.988585 106.891756 0
0 4.171374 14.0 1.0 24.988586 105.634889 0
0 4.271374 14.0 1.0 24.988586 104.335873 0
0 4.371374 14.0 1.0 24.988586 103.016130 0
0 4.471374 14.0 1.0 24.988586 101.707633 0
0 4.571374 14.0 1.0 24.988586 100.454871 0
0 4.671374 14.0 1.0 24.988586 99.314847 0
0 4.771374 14.0 1.0 24.988586 98.352198 0
0 4.871374 14.0 1.0 24.988590 89.040292 0
0 4.971374 14.0 7.0 24.990678 39.441122 0
0 5.071374 14.0 7.0 24.990681 35.988161 0
0 5.171374 14.0 7.0 24.990685 32.642016 0
0 5.271374 14.0 14.0 24.990687 31.054448 0
0 13700.000000 14.0 14.0 24.990687 31.054448 0
restarting a binary
COSMIC allows you to restart a binary from any point in its evolution from a COSMIC generated bpp array. Below we provide an example of the same evolutionary track started from the beginning and three different points in the evolution, once sometime between the beginning and the first object going supernova, once between the first and second supernova, and finally after both supernova
single_binary = InitialBinaryTable.InitialBinaries(m1=25.543645, m2=20.99784, porb=446.795757, ecc=0.448872, tphysf=13700.0, kstar1=1, kstar2=1, metallicity=0.002)
BSEDict = {'xi': 1.0, 'bhflag': 1, 'neta': 0.5, 'windflag': 3, 'wdflag': 1, 'alpha1': 1.0, 'pts1': 0.001, 'pts3': 0.02, 'pts2': 0.01, 'epsnov': 0.001, 'hewind': 0.5, 'ck': 1000, 'bwind': 0.0, 'lambdaf': 0.0, 'mxns': 3.0, 'beta': -1.0, 'tflag': 1, 'acc2': 1.5, 'remnantflag': 3, 'ceflag': 0, 'eddfac': 1.0, 'ifflag': 0, 'bconst': 3000, 'sigma': 265.0, 'gamma': -2.0, 'pisn': 45.0, 'natal_kick_array' : [[-100.0,-100.0,-100.0,-100.0,0.0], [-100.0,-100.0,-100.0,-100.0,0.0]], 'bhsigmafrac' : 1.0, 'polar_kick_angle' : 90, 'qcrit_array' : [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], 'cekickflag' : 2, 'cehestarflag' : 0, 'cemergeflag' : 0, 'ecsn' : 2.5, 'ecsn_mlow' : 1.4, 'aic' : 1, 'ussn' : 0, 'sigmadiv' :-20.0, 'qcflag' : 1, 'eddlimflag' : 0, 'fprimc_array' : [2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0], 'bhspinflag' : 0, 'bhspinmag' : 0.0, 'rejuv_fac' : 1.0, 'rejuvflag' : 0, 'htpmb' : 1, 'ST_cr' : 1, 'ST_tide' : 0, 'bdecayfac' : 1, 'randomseed' : -1235453, 'grflag' : 1, 'rembar_massloss' : 0.5, 'kickflag' : 0, 'zsun' : 0.014, 'grflag' : 1, 'bhms_coll_flag' : 0, 'don_lim' : -1, 'acc_lim' : -1}
for i in [3, 7, 11]:
bpp, bcm, initC, kick_info = Evolve.evolve(initialbinarytable=single_binary, BSEDict=BSEDict)
for column in bpp.columns:
initC = initC.assign(**{column:bpp.iloc[i][column]})
bpp_mid, bcm_mid, initC_mid, kick_info = Evolve.evolve(initialbinarytable=initC)
if i == 3:
print("From beginning")
print(bpp)
print("Started in middle at Index {0}".format(i))
print(bpp_mid)