diff --git a/.gitignore b/.gitignore
index aed2a66383de9b6d54dc10127ebd9981c138ad7f..4073a0198598387d8a284b82ecbe9ee00d9de6b0 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,3 @@
*.so
binary_c_python_api.o
-
+output/*
diff --git a/__init__.py b/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/__pycache__/defaults.cpython-36.pyc b/__pycache__/defaults.cpython-36.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..5e366ee6c1b3beae8fadebe63d1e7e413ff1ec67
Binary files /dev/null and b/__pycache__/defaults.cpython-36.pyc differ
diff --git a/defaults.py b/defaults.py
new file mode 100644
index 0000000000000000000000000000000000000000..f2e6906cd3f7e375bc915aaac1c6097b32cd6274
--- /dev/null
+++ b/defaults.py
@@ -0,0 +1,377 @@
+# File containing physics_defaults
+physics_defaults = {
+
+ # internal buffering and compression level
+ 'internal_buffering': 0,
+ 'internal_buffering_compression': 0,
+
+ # log filename
+ 'log_filename': '/dev/null',
+
+ 'metallicity': 0.02, # metallicity (default 0.02, solar)
+ 'max_evolution_time': 13700.0, # max evol time in Myr (default WMAP result)
+
+ # stellar evolution parameters
+ 'max_tpagb_core_mass': 1.38,
+ 'chandrasekhar_mass': 1.44,
+ 'max_neutron_star_mass': 1.8,
+ 'minimum_mass_for_carbon_ignition': 1.6,
+ 'minimum_mass_for_neon_ignition': 2.85,
+ 'AGB_core_algorithm': 0,
+ 'AGB_radius_algorithm': 0,
+ 'AGB_luminosity_algorithm': 0,
+ 'AGB_3dup_algorithm': 0,
+
+ # dredge up calibration (either automatic, or not)
+ 'delta_mcmin': 0.0,
+ 'lambda_min': 0.0,
+ 'minimum_envelope_mass_for_third_dredgeup': 0.5,
+
+ # minimum timestep (1yr = 1e-6 is the default)
+ 'minimum_timestep': 1e-6,
+ # maximum timestep (1 Gyr = 1e3 is the default)
+ 'maximum_timestep': 1e3,
+
+ # orbit
+ 'eccentricity': 0.0,
+
+ # tidally induced mass loss
+ 'CRAP_parameter': 0.0,
+
+ # tidal strength factor
+ 'tidal_strength_factor': 1.0,
+
+ # E2 tidal prescription. 0 = H02, 1 = Siess+2013
+ 'E2_prescription': 1,
+
+ # gravitational radiation model
+ # 0 = H02 model (Eggleton) for J and e
+ # 1 = H02 when R<RL for both stars for J and e
+ # 2 = None
+ # 3 = Landau and Lifshitz (1951) model for J (e is not changed)
+ # 4 = Landau and Lifshitz (1951) model for J when R<RL only (e is not changed)
+ 'gravitational_radiation_model': 0,
+
+ # magnetic braking multiplier
+ 'magnetic_braking_factor': 1.0,
+
+ ############################################################
+ ### Mass-loss prescriptions
+ ############################################################
+
+ # turn wind mass loss on or off
+ 'wind_mass_loss': 1,
+
+ # massive stars
+ 'wr_wind': 0, # default hurley et al wind
+ 'wr_wind_fac': 1.0, # factor applied to WR stars
+
+ # TPAGB wind details
+ 'tpagbwind': 0, # default to 0 (Karakas et al 2002)
+ 'superwind_mira_switchon': 500.0,
+ 'tpagb_reimers_eta': 1.0,
+
+ # eta for Reimers-like GB mass loss
+ 'gb_reimers_eta': 0.5,
+
+ # VW93 alterations
+ 'vw93_mira_shift': 0.0,
+ 'vw93_multiplier': 1.0,
+
+ # systemic wind angular momentum loss prescription
+ 'wind_angular_momentum_loss': 0,
+ 'lw': 1.0,
+
+ # enhanced mass loss due to rotation
+ # 0 = none = ROTATION_ML_NONE
+ # 1 = Langer+ formula (in mass-loss rate calculation,
+ # warning: can be unstable
+ # = ROTATION_ML_FORMULA)
+ # 2 = remove material in a decretion disc until J<Jcrit
+ # (ROTATION_ML_ANGMOM)
+ # 3 = 1 + 2 (not recommended!)
+ 'rotationally_enhanced_mass_loss': 2,
+ 'rotationally_enhanced_exponent': 1.0,
+
+ # timestep modulator
+ 'timestep_modulator': 1.0,
+
+ # initial rotation rates (0=automatic, >0 = in km/s)
+ 'vrot1': 0.0,
+ 'vrot2': 0.0,
+
+ ########################################
+ # Supernovae and kicks
+ ########################################
+
+ # Black hole masses:
+ # 0: H02=0
+ # 1: Belczynski
+ # 2: Spera+ 2015
+ # 3: Fryer 2012 (delayed)
+ # 4: Fryer 2012 (rapid)
+ 'BH_prescription': 2,
+ 'post_SN_orbit_method': 0,
+
+ 'wd_sigma': 0.0,
+ 'wd_kick_direction': 0,
+ 'wd_kick_pulse_number': 0,
+ 'wd_kick_when': 0,
+
+ # sn_kick_distribution and
+ # sn_kick_dispersion are only defined
+ # for SN types that leave a remnant
+ 'sn_kick_distribution_II': 1,
+ 'sn_kick_dispersion_II': 190.0,
+ 'sn_kick_distribution_IBC': 1,
+ 'sn_kick_dispersion_IBC': 190.0,
+ 'sn_kick_distribution_GRB_COLLAPSAR': 1,
+ 'sn_kick_dispersion_GRB_COLLAPSAR': 190.0,
+ 'sn_kick_distribution_ECAP': 1,
+ 'sn_kick_dispersion_ECAP': 190.0,
+ 'sn_kick_distribution_NS_NS': 0,
+ 'sn_kick_dispersion_NS_NS': 0.0,
+ 'sn_kick_distribution_TZ': 0,
+ 'sn_kick_dispersion_TZ': 0.0,
+ 'sn_kick_distribution_BH_BH': 0,
+ 'sn_kick_dispersion_BH_BH': 0.0,
+ 'sn_kick_distribution_BH_NS': 0,
+ 'sn_kick_dispersion_BH_NS': 0.0,
+ 'sn_kick_distribution_AIC_BH': 0,
+ 'sn_kick_dispersion_AIC_BH': 0.0,
+
+ 'sn_kick_companion_IA_He': 0,
+ 'sn_kick_companion_IA_ELD': 0,
+ 'sn_kick_companion_IA_CHAND': 0,
+ 'sn_kick_companion_AIC': 0,
+ 'sn_kick_companion_ECAP': 0,
+ 'sn_kick_companion_IA_He_Coal': 0,
+ 'sn_kick_companion_IA_CHAND_Coal': 0,
+ 'sn_kick_companion_NS_NS': 0,
+ 'sn_kick_companion_GRB_COLLAPSAR': 0,
+ 'sn_kick_companion_HeStarIa': 0,
+ 'sn_kick_companion_IBC': 0,
+ 'sn_kick_companion_II': 0,
+ 'sn_kick_companion_IIa': 0,
+ 'sn_kick_companion_WDKICK': 0,
+ 'sn_kick_companion_TZ': 0,
+ 'sn_kick_companion_AIC_BH': 0,
+ 'sn_kick_companion_BH_BH': 0,
+ 'sn_kick_companion_BH_NS': 0,
+
+ # evolution run splitting
+ 'evolution_splitting': 0,
+ 'evolution_splitting_sn_n': 10,
+ 'evolution_splitting_maxdepth': 1,
+
+ ########################################
+ #### Mass transfer
+ ########################################
+
+ # critical mass ratio for unstable RLOF
+ #
+ # qc = m (donor) / m (accretor) :
+ # if q>qc mass transfer is unstable
+ #
+ # H02 = Hurley et al. (2002)
+ # C14 = Claeys et al. (2014)
+
+ # non-degenerate accretors
+ 'qcrit_LMMS': 0.6944, # de Mink et al 2007 suggests 1.8, C14 suggest 1/1.44 = 0.694
+ 'qcrit_MS': 1.6, # C14 suggest 1.6
+ 'qcrit_HG': 4.0, # H02 sect. 2.6.1 gives 4.0
+ 'qcrit_GB': -1, # -1 is the H02 prescription for giants
+ 'qcrit_CHeB': 3.0,
+ 'qcrit_EAGB': -1, # -1 is the H02 prescription for giants
+ 'qcrit_TPAGB': -1, # -1 is the H02 prescription for giants
+ 'qcrit_HeMS': 3,
+ 'qcrit_HeHG': 0.784, # as in H02 2.6.1
+ 'qcrit_HeGB': 0.784, # as in H02 2.6.1
+ 'qcrit_HeWD': 3, # H02
+ 'qcrit_COWD': 3, # H02
+ 'qcrit_ONeWD': 3, # H02
+ 'qcrit_NS': 3, # H02
+ 'qcrit_BH': 3, # H02
+
+ # degenerate accretors
+ 'qcrit_degenerate_LMMS': 1.0, # C14
+ 'qcrit_degenerate_MS': 1.0, # C14
+ 'qcrit_degenerate_HG': 4.7619, # C14
+ 'qcrit_degenerate_GB': 1.15, # C14 (based on Hachisu)
+ 'qcrit_degenerate_CHeB': 3, # not used
+ 'qcrit_degenerate_EAGB': 1.15, # as GB
+ 'qcrit_degenerate_TPAGB': 1.15, # as GB
+ 'qcrit_degenerate_HeMS': 3,
+ 'qcrit_degenerate_HeHG': 4.7619, # C14
+ 'qcrit_degenerate_HeGB': 1.15, # C14
+ 'qcrit_degenerate_HeWD': 0.625, # C14
+ 'qcrit_degenerate_COWD': 0.625, # C14
+ 'qcrit_degenerate_ONeWD': 0.625, # C14
+ 'qcrit_degenerate_NS': 0.625, # C14
+ 'qcrit_degenerate_BH': 0.625, # C14
+
+ # disk wind for SNeIa
+ 'hachisu_disk_wind': 0, # 0
+ 'hachisu_qcrit': 1.15, # 1.15
+
+ # ELD accretion mass
+ 'mass_accretion_for_eld': 0.15, # 0.15, or 100(off)
+
+ # mergers have the critical angular momentum
+ # multiplied by this factor
+ 'merger_angular_momentum_factor': 1.0,
+
+ # RLOF rate method : 0=H02, 1=Adaptive R=RL, 3=Claeys et al 2014
+ 'RLOF_method': 3,
+ 'RLOF_mdot_factor': 1.0,
+
+ # RLOF time interpolation method
+ # 0 = binary_c (forward in time only), 1 = BSE (backwards allowed)
+ 'RLOF_interpolation_method': 0,
+
+ # Angular momentum in RLOF transfer model
+ # 0 : H02 (including disk treatment)
+ # 1 : Conservative
+ 'jorb_RLOF_transfer_model': 0,
+
+ # ang mom factor for non-conservative mass loss
+ # -2 : a wind from the secondary (accretor), i.e. gamma=Md/Ma (default in C14)
+ # -1 : donor (alternative in C14), i.e. gamma=Ma/Md
+ # >=0 : the specific angular momentum of the orbit multiplied by gamma
+ #
+ # (NB if Hachisu's disk wind is active, or loss is because of
+ # super-Eddington accretion or novae, material lost through the
+ # disk wind automatically gets gamma=-2 i.e. the Jaccretor)
+ 'nonconservative_angmom_gamma': -2,
+
+ # Hachisu's disc wind
+ 'hachisu_disk_wind': 0,
+ 'hachisu_qcrit': -1.0,
+
+ # donor rate limiters
+ 'donor_limit_thermal_multiplier': 1.0,
+ 'donor_limit_dynamical_multiplier': 1.0,
+
+ # RLOF assumes circular orbit (0) or at periastron (1)
+ 'rlperi': 0,
+
+ # general accretion limits
+ 'accretion_limit_eddington_multiplier': 1.0, # eddington limit factor
+ 'accretion_limit_dynamical_multiplier': 1.0, # dynamical limit factor
+ 'accretion_limit_thermal_multiplier': 1.0, # thermal limit on MS,HG and CHeB star factor
+ 'accretion_rate_novae_upper_limit': 1.03e-7,
+ 'accretion_rate_soft_xray_upper_limit': 2.71e-7,
+
+ # novae
+ 'nova_retention_method': 0,
+ 'nova_retention_fraction': 1e-3,
+ 'individual_novae': 0,
+ 'beta_reverse_nova': -1,
+ 'nova_faml_multiplier': 1.0,
+ 'nova_irradiation_multiplier': 0.0,
+
+ ########################################
+ # common envelope evolution
+ ########################################
+ 'comenv_prescription': 0, # 0=H02, 1=nelemans, 2=Nandez+Ivanova2016
+ 'alpha_ce': 1.0,
+ 'lambda_ce': -1, # -1 = automatically set
+ 'lambda_ionisation': 0.0,
+ 'lambda_enthalpy': 0.0,
+ 'comenv_splitmass': 0.0,
+ 'comenv_ms_accretion_mass': 0.0,
+ 'nelemans_minq': 0.0, # 0.0 min q for nelemans
+ 'nelemans_max_frac_j_change': 1.0, # 1.0
+ 'nelemans_gamma': 1.0, # 1.0
+ 'nelemans_n_comenvs': 1, # 1
+ 'comenv_merger_spin_method': 2,
+ 'comenv_ejection_spin_method': 1,
+ 'comenv_post_eccentricity': 0.0,
+ 'comenv_splitmass': 1.01,
+
+ # comenv accretion
+ 'comenv_ns_accretion_fraction': 0.0,
+ 'comenv_ns_accretion_mass': 0.0,
+
+ # #################################################
+ # circumbinary disc
+ # #################################################
+ # 'comenv_disc_mass_fraction' : 0.0,
+ # 'comenv_disc_angmom_fraction' : 0.0,
+ # 'cbdisc_gamma' : 1.6666666666,
+ # 'cbdisc_alpha' : 1e-6,
+ # 'cbdisc_kappa' : 1e-2,
+ # 'cbdisc_torquef' : 1.0,
+ # 'cbdisc_mass_loss_constant_rate' : 0.0,
+ # 'cbdisc_mass_loss_inner_viscous_accretion_method' : 1,
+ # 'cbdisc_mass_loss_inner_viscous_multiplier' : 1.0,
+ # 'cbdisc_mass_loss_inner_L2_cross_multiplier' : 0.0,
+ # 'cbdisc_mass_loss_ISM_ram_pressure_multiplier' : 0.0,
+ # 'cbdisc_mass_loss_ISM_pressure' : 3000.0,
+ # 'cbdisc_mass_loss_FUV_multiplier' : 0.0,
+ # 'cbdisc_mass_loss_Xray_multiplier' : 1.0,
+ # 'cbdisc_viscous_photoevaporation_coupling' : 1,
+ # 'cbdisc_inner_edge_stripping' : 1,
+ # 'cbdisc_outer_edge_stripping' : 1,
+ # 'cbdisc_minimum_luminosity' : 1e-4,
+ # 'cbdisc_minimum_mass' : 1e-6,
+ # 'cbdisc_eccentricity_pumping_method' : 1,
+ # 'cbdisc_resonance_multiplier' : 1.0,
+ # 'comenv_post_eccentricity' : 1e-5,
+
+ ##################################################
+ # wind accretion
+ ##################################################
+
+ # Bondi-Hoyle accretion multiplier
+ 'Bondi_Hoyle_accretion_factor': 1.5,
+ # Wind-RLOF method: 0=none, 1=q-dependent, 2=quadratic
+ # (See Abate et al. 2012,13,14 series of papers)
+ 'WRLOF_method': 0,
+
+ # pre-main sequence evolution
+ 'pre_main_sequence': 0,
+ 'pre_main_sequence_fit_lobes': 0,
+
+ ########################################
+ # Nucleosynthesis
+ ########################################
+
+ # lithium
+ 'lithium_hbb_multiplier': 1.0,
+ 'lithium_GB_post_1DUP': 0.0,
+ 'lithium_GB_post_Heflash': 0.0,
+
+ ########################################
+ # Yields vs time (GCE)
+ ########################################
+
+ 'yields_dt': 100000,
+ 'escape_velocity': 1e9, # if wind v < this, ignore the yield
+ # and assume it is lost to the IGM
+ 'escape_fraction': 0.0, # assume all yield is kept in the population
+
+ # minimum sep/per for RLOF on the ZAMS
+ 'minimum_separation_for_instant_RLOF': 0,
+ 'minimum_orbital_period_for_instant_RLOF': 0,
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/test_david.py b/test_david.py
index 59c22675547d5ea27854be53e79981dffdaccede..ba36dec3ac419778da32216117253a201d0fc458 100644
--- a/test_david.py
+++ b/test_david.py
@@ -1,8 +1,11 @@
#!/usr/bin/python3
import os
-
import binary_c
+from defaults import physics_defaults
+
+
+
############################################################
# Test script to run a binary using the binary_c Python
@@ -73,9 +76,32 @@ def run_simple_loop_binary():
with open("test_output.txt", "a") as myfile:
myfile.write(' '.join(value_array[:3])+'\n')
- print(mass_1)
- print(initial_m1)
- print(time)
+ print('mass1:', mass_1)
+ print('zams_mass1:' ,initial_m1)
+ print('time:', time)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
# print("Current binary_c object class functions: ")
@@ -96,5 +122,44 @@ def run_simple_loop_binary():
# run_test_binary()
# Test grid-like
-run_simple_loop_binary()
+# run_simple_loop_binary()
+
+def create_arg_string(arg_dict):
+ """
+ Function that creates the arg string
+ """
+
+ arg_string = ''
+ for key in arg_dict.keys():
+ arg_string += f'{key} {arg_dict[key]} '
+ arg_string = arg_string.strip()
+ return arg_string
+
+def example_with_loading_default_args():
+ """
+ Example function loading the default physics args for a binary_c system. Got
+ it from the binary_grid2 perl module
+ """
+
+ # Load args
+ physics_args = physics_defaults.copy()
+
+ # Manually set M_1, M_2, orbital_period and separation values:
+ physics_args['M_1'] = 20
+ physics_args['M_2'] = 15
+ physics_args['separation'] = 0 # 0 = ignored, use period
+ physics_args['orbital_period'] = 4530.0
+
+ arg_string = create_arg_string(physics_args)
+
+ arg_string = f'binary_c {arg_string}'
+
+ output = binary_c.run_binary(arg_string)
+ print (output)
+
+
+
+example_with_loading_default_args()
+
+quit()