scottcoughlin2014
commented
5 years ago This will get you the bpp and bcm array for an example event with the most recent version of the code.
from cosmic.sample.initialbinarytable import InitialBinaryTable
from cosmic.evolve import Evolve
import pandas as pd
import numpy as np
from cosmic import _evolvebin
BSEDict = {'neta': 0.5, 'bwind': 0.0, 'hewind': 1.0, 'alpha1': 1.0, 'lambdaf': 1.0, 'ceflag': 0,
'tflag': 1, 'ifflag': 0, 'wdflag': 0, 'bhflag': 1, 'nsflag': 3, 'mxns': 3.0, 'pts1': 0.05,
'pts2': 0.01, 'pts3': 0.02, 'sigma': 265.0, 'beta': -1.0, 'xi': 0.5, 'acc2': 1.5, 'epsnov': 0.001,
'eddfac': 1.0, 'gamma': -2.0, 'bconst': -3000, 'CK': -1000, 'merger': 0, 'windflag': 3,
'B_0': [0.0, 0.0], 'bacc': [0.0, 0.0], 'tacc': [0.0, 0.0],
'bkick': [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 'massc': [0.0, 0.0],
'opsin': [0.0, 0.0], 'epoch': [0.0, 0.0], 'tphys': 0.0, 'ppsn': 1}
single_binary = InitialBinaryTable.SingleBinary(m1=3.518899, m2=3.508906, porb=3.811531, ecc=0.020918,
tphysf=5016.353066,
kstar1=1.0, kstar2=1.0, metallicity=0.02)
bpp, bcm, initC = Evolve.evolve(initialbinarytable=single_binary, BSEDict=BSEDict)
initC['randomseed'] = 243848
bpp_columns = ['tphys', 'mass_1', 'mass_2', 'kstar_1', 'kstar_2' ,
'sep', 'ecc', 'RROL_1', 'RROL_2', 'evol_type', 'bin_num']
bcm_columns = ['tphys', 'kstar_1', 'mass0_1', 'mass_1', 'lumin_1', 'rad_1',
'teff_1', 'massc_1', 'radc_1', 'menv_1', 'renv_1', 'epoch_1',
'ospin_1', 'deltam_1', 'RROL_1', 'kstar_2', 'mass0_2', 'mass_2',
'lumin_2', 'rad_2', 'teff_2', 'massc_2', 'radc_2', 'menv_2',
'renv_2', 'epoch_2', 'ospin_2', 'deltam_2', 'RROL_2',
'porb', 'sep', 'ecc', 'B_0_1', 'B_0_2', 'SN_1',
'SN_2', 'bin_state', 'merger_type', 'bin_num']
f = np.array(initC)
f = f[0]
[bpp, bcm] = _evolvebin.evolv2(f[0], f[1], f[2], f[3], f[4], f[5], f[6], f[7], f[8], f[9],
f[10], f[11], f[12], f[13], f[14], f[15], f[16], f[17], f[18], f[19],
f[20], f[21], f[22], f[23], f[24], f[25], f[26], f[27], f[28], f[29],
f[30], f[31], f[32], f[33], f[34], f[35], f[36])
bpp = bpp[:np.argwhere(bpp[:,0] == -1)[0][0]]
bcm = bcm[:np.argwhere(bcm[:,0] == -1)[0][0]]
bpp_bin_numbers = np.atleast_2d(np.array([f[37]] * len(bpp))).T
bcm_bin_numbers = np.atleast_2d(np.array([f[37]] * len(bcm))).T
bpp = np.hstack((bpp, bpp_bin_numbers))
bcm = np.hstack((bcm, bcm_bin_numbers))
bpp = pd.DataFrame(bpp,
columns=bpp_columns,
index=bpp[:, -1].astype(int))
bcm = pd.DataFrame(bcm,
columns=bcm_columns,
index=bcm[:, -1].astype(int))
bcm.merger_type = bcm.merger_type.astype(int).astype(str).apply(lambda x: x.zfill(4))
bcm.bin_state = bcm.bin_state.astype(int)
bpp.bin_num = bpp.bin_num.astype(int)
bcm.bin_num = bcm.bin_num.astype(int)
print(bcm)
print(bpp)
katiebreivik
In hrdiag.f line 762ish:
We need to figure out what the actual remnant should be for a failed ECSN. I think it might just be a He-star that cools over time. I don't see why it would be a He-WD, so we need to decide how to handle this since it is borking out in evolv2.f at line 3389:
Since this happens BEFORE tphys >= tphysf, the binary is not tracked and is lost from the population. Should we set it to kw=7? (I don't think this makes sense because it's not a He MS star) Should we create a kw=16 for failed ECSN?
Let's get some input from folks who know better :-)