Doxygen

CMC

cmc.c

Defines

_MAIN_

Functions

int main(int argc, char *argv[])

The main function.

Monte Carlo (MC) methods calculate the dynamical evolution of a collisional system of N stars in the Fokker-Planck approximation, which applies when the evolution of the cluster is dominated by two-body relaxation, and the relaxation time is much larger than the dynamical time. In practice, further assumptions of spherical symmetry and dynamical equilibrium have to be made. The Henon MC technique (Henon 1971), which is based on orbit averaging, represents a balanced compromise between realism and speed. The MC method allows for a star-by-star realization of the cluster, with its N particles representing the N stars in the cluster. Integration is done on the relaxation timescale, and the total computational cost scales as O(NlogN) (Henon 1971).

Our code here is based on the Henon-type MC cluster evolution code CMC (Cluster Monte Carlo), developed over many years by Joshi et al. (2000, 2001), Fregeau et al. (2003), Fregeau & Rasio (2007), Chatterjee et al. (2010), and Umbreit et al. (2012). CMC includes a detailed treatment of strong binary star interactions and physical stellar collisions (Fregeau & Rasio 2007), as well as an implementation of single and binary star evolution (Chatterjee et al. 2010) and the capability of handling the dynamics around a central massive black hole (Umbreit et al. 2012).

Code Overview:

The principal data structure is an array of structures of type star_t. One or more of these stars can be binaries. Binaries are stored in a separate array of structures of type binary_t. If an element in the star array is a binary, it’s binind property/variable is assigned a non-zero value. Moreover, the value of this variable is assigned in such a way that it indicates the index of the binary array which holds the properties of this binary.

CMC can be run both serially as well as parallely.

Parallelization Strategy:

There are various routines which have varying dependencies and accesses between elements of the data structures. To minmize inter-process communication required by these routines as much as possible, we partition the data such that the number of stars held by each processor is a multiple of MIN_CHUNK_SIZE (input parameter, defaults to 20, refer http://adsabs.harvard.edu/abs/2013ApJS..204…15P for a justification for the choice of this value, and detailed explanation). The remainder of stars which is < MIN_CHUNK_SIZE after division goes to the last processor.

Most routines do computations on a star by star basis, and access only the properties of the star being processed. However, the values of gravitational potential, masses of the stars and radial positions of the stars are accessed in a complex, data-dependent manner throughout the code. Hence, we store these values in separate arrays and duplicate/copy these across all processors. Since these properties change during a timestep, we synchronize these arrays at the end of each time step.

Like mentioned above, most routines perform computations on a star by star basis. The only routine that differs from this pattern is the sorting routine that sorts all the stars by their radial positions. We use a parallel sorting algorithm called Sample Sort (see http://adsabs.harvard.edu/abs/2013ApJS..204…15P for a detailed discussion).

The parallel sorting routine does the sorting and automatically redistributes the data to the processors, however, there is no control over the number of stars that will end up on each processor. So, after the sorting takes place, we exchange data between processors to adhere to the data partitioning scheme mentioned above.

OPTIONS: -q —quiet : do not print diagnostic info to stdout -d —debug : turn on debugging -V —version : print version info -h —help : display this help text -s —streams :Run with multiple random streams. To mimic the parallel version with the given number of processors

Parameters

argc – Input argument count
argv[] – Input argument list: USAGE: <path>/cmc [options…] <input_file> <output_file_prefix> mpirun -np <procs> <path>/cmc <input_file> <output_file_prefix>

Returns

indicates how the program exited. 0 implies normal exit, and nonzero value implies abnormal termination.

cmc_bhlosscone.c

Functions

void create_rwalk_file(char *fname)

?

Parameters: fname – file name

void write_rwalk_data(char *fname, long index, double Trel, double dt, double l2_scale, double n_steps, double beta, double n_local, double W, double P_orb, double n_orb)

?

Parameters

fname – ?
index – star index
Trel – ?
dt – ?
l2_scale – ?
n_steps – ?
beta – ?
n_local – ?
W – ?
P_orb – ?
n_orb – ?

void bh_rand_walk(long index, double v[4], double vcm[4], double beta, double dt)

This is the random walk procedure as outlined by Freitag & Benz (2002). Change of notation: beta here is theta in the paper.

Parameters

index – star index
v[4] – ?
vcm[4] – ?
beta – ?
dt – ?

void get_3d_velocities(double *w, double vr, double vt)

?

Parameters

w – ?
vr – ?
vt – ?

void do_random_step(double *w, double beta, double delta)

?

Parameters

w – ?
beta – ?
delta – ?

double check_angle_w_w_new(double *w, double *w_new, double delta)

calculate the angle between w and w_new and compare to deltat

Parameters

w – ?
w_new – ?
delta – ?

Returns

?

double calc_P_orb(long index)

calculate star’s radial orbital period

Parameters: index – star index
Returns: star’s radial orbital period

double calc_p_orb_f(double x, void *params)

integrand for calc_P_orb

Parameters

x – ?
params – ?

Returns

?

double calc_p_orb_f2(double x, void *params)

integrand for calc_P_orb

Parameters

x – ?
params – ?

Returns

?

double calc_p_orb_gc(double x, void *params)

integrand for calc_P_orb, using Gauss-Chebyshev for regularizing the integrand near the endpoints

Parameters

x – ?
params – ?

Returns

?

struct Interval get_r_interval(double r)

?

Parameters: r – ?
Returns: ?

cmc_binbin.c

Functions

void bb_calcunits(fb_obj_t *obj[2], fb_units_t *bb_units)

calculate the units used

Parameters

obj[2] – ?
bb_units – ?

fb_ret_t binbin(double *t, long k, long kp, double W, double bmax, fb_hier_t *hier, gsl_rng *rng)

the main attraction

Parameters

t – ?
k – index of star 1
kp – index of star 2
W – ?
bmax – ?
hier – ?
rng – gsl rng

Returns

?

cmc_binsingle.c

Functions

void bs_calcunits(fb_obj_t *obj[2], fb_units_t *bs_units)

calculate the units used

Parameters

obj[2] – ?
bs_units – ?

fb_ret_t binsingle(double *t, long ksin, long kbin, double W, double bmax, fb_hier_t *hier, gsl_rng *rng)

?

Parameters

t – ?
ksin – ?
kbin – ?
W – ?
bmax – ?
hier – ?
rng – ?

Returns

?

cmc_bse_utils.c

Functions

void update_bse_from_sse(bse_t *bvars, sse_t *svars, int bmember)

?

Parameters

bvars – ?
svars – ?
bmember – ?

void update_sse_from_bse(bse_t *bvars, sse_t *svars, int bmember)

?

Parameters

bvars – ?
svars – ?
bmember – ?

void get_sse_from_star(sse_t *svars, star_t *star)

?

Parameters

svars – ?
star – ?

void update_star_from_sse(star_t *star, sse_t svars)

?

Parameters

star – ?
svars – ?

void get_bse_from_binary(bse_t *bvars, binary_t *star)

?

Parameters

bvars – ?
star – ?

void update_binary_from_bse(binary_t *star, bse_t *bvars)

?

Parameters

star – ?
bvars – ?

void compress_binary(star_t *bincom, binary_t *bin)

Breaks a binary if it contains a k=15 “star”, i.e. a massless remnant.

Use this function only if STELLAR_EVOLUTION is turned on (Spares one an extra if-statement)!

Parameters

bincom – center-of-mass properties
bin – binary properties

void cmc_bse_comenv(binary_t *tempbinary, double cmc_l_unit, double RbloodySUN, double *zpars, double *vs, int *fb)

?

Parameters

tempbinary – ?
cmc_l_unit – ?
RbloodySUN – ?
zpars – ?
vs – ?
fb – ?
ecsnp – ?
ecsn_mlow – ?
ST_tide – ?

cmc_core.c

Functions

static inline int is_member(int k, int *set, int len)

?

Parameters

k – ?
set – ?
len – ?

Returns

?

struct densities density_estimators(int n_points, int *startypes, int len)

Calculates density estimators including only stars with certain types. (mpi version of density_estimators)

Parameters

n_points – the number of points (stars) over which the local density is estimated (Casertano & Hut ‘85 suggest n_points=6)
startypes – Array containing the star types that should be considered.
len – The length of the startype array.

Returns

The array with the estimated local densities. The array is malloc’ed so you have to free it yourself.

struct core_t core_properties(struct densities rhoj)

calculate core quantities using density weighted averages (note that in Casertano & Hut (1985) only rho and rc are analyzed and tested) (mpi version of core_properties)

Parameters: rhoj – ?
Returns: ?

void append_core_header(FILE *cfile, char *tag, int core_number)

Appends the header information for a certain core to the core file.

No newline character will be added, so you have to do it yourself once all core headers are written.

Parameters

cfile – The core file.
tag – Shorthand for the type of core.
core_number – The position of the core within the file (used to calculate the column numbers)

void write_core_data(FILE *cfile, struct core_t core)

Writes out all core data to a file.

No newline character will be added, so you have to do it yourself once all core data is written.

Parameters

cfile – The core file (open file handle).
core – The struct that holds all core data.

struct core_t no_remnants_core(int n_points)

?

Parameters: n_points – ?
Returns: ?

cmc_dynamics.c

Functions

void dynamics_apply(double dt, gsl_rng *rng)

core of the code: applies relaxation, does single-single collisions and binary interactions

Parameters

dt – timestep
rng – gsl rng

cmc_dynamics_helper.c

Functions

void zero_star(long j)

Zero’es out all star variables corresponding to given index.

Parameters: j – star index

void zero_binary(long j)

Zero’es out all variables of the binary corresponding to the given index.

Parameters: j – index of binary

void print_interaction_status(char status_text[])

prints the specified interaction status text

Parameters: status_text[] – status text

void print_interaction_error(void): prints error message during interaction

void destroy_obj(long i)

destroy an object (can be star or binary)

Parameters: i – Index of star to be destroyed

void destroy_obj_new(long i)

destroy an object (can be star or binary) - pure serial version

Parameters: i – index of object/star to be destroyed

void destroy_binary(long i)

destroy a binary

Parameters: i – index of binary (in the binary array) to be destroyed

long create_star(int idx, int dyn_0_se_1)

create a new star, returning its index

Parameters

idx – index of star that creates new star
dyn_0_se_1 – if 0, created by dynamics, if 1 created by stellar evolution

Returns

index of new star

long create_binary(int idx, int dyn_0_se_1)

create a new binary, returning its index. A new binary creation is done as follows. Since binaries get destroyed, so, first we scan the existing binary array and look for holes i.e. which were left behind by destroyed stars. If one is found, we insert a new binary there, if not we insert it at the end.

Parameters

idx – index of the star that is creating the binary
dyn_0_se_1 – 0 if created by dynamics, 1 if created by stellar evolution

Returns

index of new binary

void sort_three_masses(long sq, long *k1, long *k2, long *k3)

Sort by mass: k1, k2, k3 in order of most massive to least massive. (parallel version of sort_three_masses)

Parameters

sq – index of the first of the three consecutive stars
k1 – variable to store the most massive star’s index
k2 – variable to store the 2nd most massive star’s index
k3 – variable to store the least massive star’s index

void calc_3bb_encounter_dyns(long k1, long k2, long k3, double v1[4], double v2[4], double v3[4], double (*vrel12)[4], double (*vrel3)[4], gsl_rng *rng)

Function for calculating quantities relevant for three-body binary formation Generates random direction for vt for all stars (so then have 3D veloc. for each) Then, for the two possible binaries that could form (1,2 or 1,3) we calculate: (1) COM veloc of candidate binary pair (2) v3: veloc of 3rd star relative to binary In COM frame of the system of three stars, calculate total momentum and energy.

Parameters

k1 – index of 1st star
k2 – index of 2nd star
k3 – index of 3rd star
v1[4] – ?
v2[4] – ?
v3[4] – ?
vrel12 – ?
vrel3 – ?
rng – gsl rng

double get_eta(double eta_min, long k1, long k2, long k3, double vrel12[4], double vrel3[4])

?

Parameters

eta_min – ?
k1 – index of star 1
k2 – index of star 2
k3 – index of star 3
vrel12[4] – ?
vrel3[4] – ?

Returns

?

void make_threebodybinary(double P_3bb, long k1, long k2, long k3, long form_binary, double eta_min, double ave_local_mass, double n_local, double sigma_local, double v1[4], double v2[4], double v3[4], double vrel12[4], double vrel3[4], double delta_E_3bb, gsl_rng *rng)

?

Parameters

P_3bb – ?
k1 – index of star 1
k2 – index of star 2
k3 – index of star 3
form_binary – ?
eta_min – ?
ave_local_mass – ?
n_local – ?
sigma_local – ?
v1[4] – ?
v2[4] – ?
v3[4] – ?
vrel12[4] – ?
vrel3[4] – ?
delta_E_3bb – ?
rng – gsl rng

void calc_sigma_local(long k, long p, long N_LIMIT, double *ave_local_mass, double *sigma_local)

?

Parameters

k – ?
p – ?
N_LIMIT – ?
ave_local_mass – ?
sigma_local – ?

long star_get_id_new(void)

generate unique star id’s (original version). Returns successive ids after N.

Returns: unique id for newly created star

long star_get_merger_id_new(long id1, long id2)

generate unique star id’s. The serial version of this function star_get_id_new generated new ids incrementally starting with N. But, we cant do this in the parallel version since this would need additional communication to keep the latest assigned id synchronized across processors, and also will result in a race condition. So, we use the following formula new_id = N + min(id1N, id2N) + id1/N + id2/N. Here / implies integer division, and % is the modulo.

Parameters

id1 – id of first star
id2 – id of second star

Returns

new id

double calc_n_local(long k, long p, long N_LIMIT)

calculate local density by averaging

Parameters

k – star index
p – number of stars to average over
N_LIMIT – typically the total number of stars

Returns

local density for the given star

double calc_Ai_local(long k, long kp, long p, double W, long N_LIMIT)

?

Parameters

k – index of star 1
kp – index of star 2
p – ?
W – ?
N_LIMIT – ?

Returns

?

void calc_encounter_dyns(long k, long kp, double v[4], double vp[4], double w[4], double *W, double *rcm, double vcm[4], gsl_rng *rng, int setY)

?

Parameters

k – index of star 1
kp – index of star 2
v[4] – ?
vp[4] – ?
w[4] – ?
W – ?
rcm – ?
vcm[4] – ?
rng – ?
setY – ?

void set_star_EJ(long k)

Sets E and J values of the star based on vr and vt values, current position and potential.

Parameters: k – index of star

void set_star_news(long k)

copies some variables of the star to the new variables

Parameters: k – index of star

void set_star_olds(long k)

copies some variables of the star to the old variables

Parameters: k – index of star

void copy_globals_to_locals(long k)

copies the global/duplicated array values to the local star structure

Parameters: k – index of star

double binint_get_mass(long k, long kp, long id)

find masses of merging stars from binary interaction components

Parameters

k – index of first star
kp – index of second star
id – id of star

Returns

masses of merging stars

double binint_get_spins(long k, long kp, long id)

find spins of merging black holes from binary interaction components

Parameters

k – index of first star
kp – index of second star
id – id of star

Returns

spin of merging black holes

long binint_get_startype(long k, long kp, long id)

Sourav: find stellar types of merging stars from binary interaction components.

Parameters

k – index of star 1
kp – index of star 2
id – star id to be matched

Returns

?

double binint_get_createtime(long k, long kp, long id)

Sourav: Sourav: toy rejuvenation- finding creation times of binary interaction components.

Parameters

k – index of star 1
kp – index of star 2
id – star id to be matched

Returns

?

double binint_get_lifetime(long k, long kp, long id)

Sourav: Sourav: toy rejuvenation- finding MS lifetimes of binary interaction components.

Parameters

k – index of star 1
kp – index of star 2
id – star id to be matched

Returns

?

long binint_get_indices(long k, long kp, long id, int *bi)

return star index of star with id “id” from binary interaction components, along with which member

Parameters

k – index of star 1
kp – index of star 2
id – id of star to be matched
bi – “bi=0,1” if a binary, else single

Returns

return star index of star with id “id” from binary interaction components, along with which member

void binint_log_obj(fb_obj_t *obj, fb_units_t units)

?

Parameters

obj – ?
units – ?

void binint_log_status(fb_ret_t retval, double vesc)

?

Parameters: retval – ?

void binint_log_collision(const char interaction_type[], long id, double mass, double r, fb_obj_t obj, long k, long kp, long startype)

?

Parameters

interaction_type[] – ?
id – ?
mass – ?
r – ?
obj – ?
k – index of star 1
kp – index of star 2
startype – star type

void binint_do(long k, long kp, double rperi, double w[4], double W, double rcm, double vcm[4], gsl_rng *rng)

do binary interaction (bin-bin or bin-single)

Parameters

k – index of 1st star
kp – index of 2nd star
rperi – ?
w[4] – ?
W – ?
rcm – ?
vcm[4] – ?
rng – gsl rng

double simul_relax_new(void)

Parallel version of simul_relax_new.

Returns: relaxation timestep

double simul_relax(gsl_rng *rng)

simulate relaxation to get timestep (original serial version). Calculates timestep for each star using some average quantities taken around the star, and then returns the minimum of these.

Parameters: rng – gsl rng
Returns: relaxation timestep

int destroy_bbh(double m1, double m2, double a, double e, double nlocal, double sigma, struct rng_t113_state *rng_st): Here we need to compare the inspiral time for binary black holes to the typical timescale for another encounter. Note that I’m not being careful with units, since we only need to compare timescales.

void break_wide_binaries(struct rng_t113_state *rng_st): Since the binary interactions are done in a vaccuum, it is possible for them to produce pathologically wide binaries, which must be broken by hand, else they shorten the timestep to a crawl.

void calc_sigma_r(long p, long N_LIMIT, double *sig_r_or_mave, double *sig_sigma, long *sig_n, int r_0_mave_1)

Computes the local average velocity dispersion value for each star (parallel version of calc_sigma_r).

Parameters

p – the window to be averaged over
N_LIMIT – total number of stars in the processor
sig_r_or_mave – gets filled up with either the radial positions or average masses
sig_sigma – sigma array
sig_n – n value of sigma structure
r_0_mave_1 – if 0, sig_r_or_mave is filled with r values, if not, average mass values

double calc_average_mass_sqr(long index, long N_LIMIT)

calculates sliding averages of mass^2 around given index

Parameters

index – index of star around which average is needed
N_LIMIT – total number of stars

Returns

average of mass^2

double floateq(double a, double b)

generic routine for testing for approximate equality of floating point numbers

Parameters

a – first floating point number
b – second floating point number

Returns

?

double sigma_tc_nd(double n, double m1, double r1, double m2, double vinf)

tidal capture (including merger) cross section; inputs are assumed to be in (self-consistent) code units

Parameters

n – ?
m1 – ?
r1 – ?
m2 – ?
vinf – ?

Returns

?

double sigma_tc_nn(double na, double ma, double ra, double nb, double mb, double rb, double vinf)

tidal capture (including merger) cross section; inputs are assumed to be in (self-consistent) code units

Parameters

na – ?
ma – ?
ra – ?
nb – ?
mb – ?
rb – ?
vinf – ?

Returns

?

double Tl(int order, double polytropicindex, double eta)

T_l function for use in tidal capture calculations.

Parameters

order – ?
polytropicindex – ?
eta – ?

Returns

?

double Etide(double rperi, double Mosc, double Rosc, double nosc, double Mpert)

tidal energy of Mpert acting on Mosc, in the polytropic approximation

Parameters

rperi – ?
Mosc – ?
Rosc – ?
nosc – ?
Mpert – ?

Returns

?

void vt_add_kick(double *vt, double vs1, double vs2, struct rng_t113_state *rng_st)

Add kicks to vt.

Parameters

vt – transverse velocity
vs1 – ?
vs2 – ?
rng_st – random state

void binary_bh_merger(long k, long kb, long knew, int kprev0, int kprev1, struct rng_t113_state *rng_st)

add GW recoil kicks and mass loss for mergers of BBHs

These are all based on fits to NR simulations, though the functional forms I’ve taken from Section V of Gerosa and Kesden, PRD, 93, 12, 124066 (2016)

cmc_evolution_thr.c

Functions

orbit_rs_t calc_orbit_rs(long si, double E, double J)

function to calculate r_p, r_a, dQ/dr|_r_p, and dQ/dr|_r_a for an orbit

Parameters

si – index of star
E – energy
J – angular momentum

Returns

orbit structure

double get_Tbb(central_t central)

computes the binary-binary interaction timestep

Parameters: central – properties of central stars of the cluster
Returns: bin-bin timestep

double get_Tbs(central_t central)

computes the binary-single interaction timestep

Parameters: central – properties of central stars of the cluster
Returns: bin-sin timestep

double GetTimeStep(gsl_rng *rng)

Computes the timesteps for relaxation, collision, binary interactions, stellar evolution etc, and returns the minimum of them as the timestep for the current iteration.

Parameters: rng – gsl rng
Returns: the timestep

void tidally_strip_stars(void): removes tidally-stripped stars

void remove_star(long j, double phi_rtidal, double phi_zero)

removes star

Parameters

j – index of star
phi_rtidal – potential at tidal radius
phi_zero – potential at zero

void count_esc_bhs(long j)

Meagan: output file for escaping BHs.

Parameters: j – star index

void remove_star_center(long j)

sents star’s position to infinity, sets mass to a very small number, and sets up vr and vt for future calculations.

Parameters: j – index of star

void get_positions_loop(struct get_pos_str *get_pos_dat)

Requires indexed (sorted in increasing r) stars with potential computed in star[].phi and N_MAX set. Uses sE[], sJ[] and sr[] from previous iteration. Returns positions and velocities in srnew[], svrnew[], and svtnew[]. Returns Max r for all stars, or -1 on error.

Parameters: get_pos_dat – ?

double get_positions()

return maximum stellar radius, get r_p and r_a for all objects

Returns: maximum stellar radius

cmc_fits.c

Functions

void load_binary_data(int d_k, int s_k)

Copies the properties corresponding to index s_k of the cfd structure to the index d_k of the binary array.

Parameters

d_k – destination index in the binary array where the data needs to be copied to
s_k – source index in the cfd structure from where data needs to be copied

void load_fits_file_data(void): Populates the star and binary arrays from the data obtained from the cfd structure. Each processor just reads and stores the corresponding data from the cfd structure it is responsible for, and ignores the rest.

cmc_funcs.c

Warning

doxygenfile: Cannot find file “cmc_funcs.c

cmc_init.c

Warning

doxygenfile: Cannot find file “cmc_init.c

cmc_io.c

Defines

PRINT_PARSED(DOC)

CHECK_PARSED(A)

CHECK_PARSED(A, DEFAULT)

CHECK_PARSED(A, DEFAULT, DOC)

Functions

void print_version(FILE *stream)

print the version

Parameters: stream – stream to be printed to

void cmc_print_usage(FILE *stream, char *argv[])

print the usage

Parameters

stream – stream to be printed to
argv[] – input arg list

void print_results(void): Writes output to stdout and all files required for post-simulation analysis.

void print_2Dsnapshot(void): prints a 2D snapshot

void print_bh_snapshot(void): prints BH snapshots

void PrintLogOutput(void): Calculates some quantities required for output and prints out log files. These are mostly files that need data only from the root node to be written out.

void PrintFileOutput(void): prints out some more files

void print_bh_summary(): Meagan: extra output for bhs.

void print_esc_bh_summary(): Meagan - extra output for bhs.

void find_nstars_within_r(double r, long *ns, long *nb)

finds number of single and binary stars within a given radial position

Parameters

r – radial position
ns – variable to hold/return number of single stars
nb – variable to hold/return number of binary stars

void PrintParaFileOutput(void): This is the function which actually writes the parallel file buffers into the file in parallel using MPI IO.

void mpi_para_file_write(char *wrbuf, long long *len, long long *prev_cum_offset, MPI_File *fh)

Flushes out data in parallel present in the char buffer to the corresponding file using MPI-IO.

Parameters

wrbuf – write buffer containing the data to be flushed out
len – length/size of the data in the buffer
prev_cum_offset – offset of the file where the data needs to be written
fh – MPI-IO File handle

void parser(int argc, char *argv[], gsl_rng *r)

Parsing of Input Parameters / Memory allocation / File I/O.

Parameters

argc – Input argument count from main()
argv[] – Input argument list from main()
r – gsl rng

void close_buffers(void): close file buffers/pointers

void close_root_buffers(void): Closes some of the file pointers - of files which require writing only by the root node.

void close_node_buffers(void): Closes some of the file pointers - of files which require writing by the all nodes.

void mpi_close_node_buffers(void): Closes the MPI file pointers - of files which require writing only by the all nodes.

void trap_sigs(void): traps signals

void print_initial_binaries(void): print out initial binary paramaeters to data file

void print_conversion_script(void): print handy script for converting output files to physical units

char *sprint_star_dyn(long k, char string[MAX_STRING_LENGTH])

routines for printing star/binary info in a unified log format

Parameters

k – index of star
string[MAX_STRING_LENGTH] – ?

Returns

?

char *sprint_bin_dyn(long k, char string[MAX_STRING_LENGTH])

?

Parameters

k – index of star
string[MAX_STRING_LENGTH] – ?

Returns

?

void parse_snapshot_windows(char *param_string)

?

Parameters: param_string – ?

void print_snapshot_windows(void): ?

int valid_snapshot_window_units(void)

?

Returns: ?

void write_snapshot(char *filename, int bh_only, char *tablename)

writes out snapshot to the given file

Parameters

filename – name of the file
bh_only – if bh_only>0 this’ll print only BHs.

void print_denprof_snapshot(char *infile)

smaller snapshot outputting limited data.

Parameters: infile – file name

void get_star_data(int argc, char *argv[], gsl_rng *rng)

Does a miscellaneous set of things, primarily copying the data read from the input file into the star and binary data structures. Also allocates duplicate array for the parallel version. Also, initializes a few global variables and the random number generator.

Parameters

argc – input arg count
argv[] – input arg list
rng – random number generator

void mpiInitGlobArrays(): Allocate global arrays based on total number of stars.

void load_dynamical_friction_data()

void load_tidal_tensor()

void save_global_vars(restart_struct_t *rest)

void load_global_vars(restart_struct_t *rest)

void save_restart_file()

void load_restart_file()

struct restart_struct_t

Public Members

double s_Eescaped

double s_Jescaped

double s_Eintescaped

double s_Ebescaped

double s_TidalMassLoss

double s_Etidal

double s_Prev_Dt

long long s_mpi_logfile_len

long long s_mpi_escfile_len

long long s_mpi_binintfile_len

long long s_mpi_collisionfile_len

long long s_mpi_tidalcapturefile_len

long long s_mpi_semergedisruptfile_len

long long s_mpi_removestarfile_len

long long s_mpi_relaxationfile_len

long long s_mpi_pulsarfile_len

long long s_mpi_morepulsarfile_len

long long s_mpi_triplefile_len

long long s_mpi_bhmergerfile_len

long long s_mpi_logfile_ofst_total

long long s_mpi_escfile_ofst_total

long long s_mpi_binaryfile_ofst_total

long long s_mpi_binintfile_ofst_total

long long s_mpi_collisionfile_ofst_total

long long s_mpi_tidalcapturefile_ofst_total

long long s_mpi_semergedisruptfile_ofst_total

long long s_mpi_removestarfile_ofst_total

long long s_mpi_relaxationfile_ofst_total

long long s_mpi_pulsarfile_ofst_total

long long s_mpi_morepulsarfile_ofst_total

long long s_mpi_triplefile_ofst_total

long long s_mpi_bhmergerfile_ofst_total

double s_OldTidalMassLoss

double s_TidalMassLoss_old

long s_N_bb

long s_N_bs

double s_E_bb

double s_E_bs

long s_Echeck

int s_se_file_counter

long s_snap_num

long s_bh_snap_num

long s_StepCount

long s_tcount

double s_TotalTime

long s_newstarid

double s_cenma_m

double s_cenma_m_new

double s_cenma_e

double s_cenma_e_new

cmc_mpi.c

Functions

void mpiFindIndicesEvenGen(long N, int i, int *mpiStart, int *mpiEnd)

Function to find start and end indices for each processor for a loop over N. Makes sure the number of elements each proc works on is an even number. Is a generic function that does not use myid, instead it should be passed through parameter i.

Parameters

N – data size for which data partitioning scheme needs to be found
i – parameters for processor i that is needed
mpiStart – start index of processor i in the global dataset
mpiEnd – end index of processor i in the global dataset

void mpiFindIndicesSpecial(long N, int *mpiStart, int *mpiEnd)

Function specially made for our code as dynamics_apply() takes 2 stars at a time. This function divides particles in pairs if the no.of stars is even, but if the no.of stars is odd, it gives one star to the last processor (since only this will be skipped by the serial code, and hence will help comparison) and divides the rest in pairs.

Parameters

N – data size for which data partitioning scheme needs to be found
mpiStart – array for start indices
mpiEnd – array for end indices

void mpiFindIndicesSpecialGen(long N, int i, int *mpiStart, int *mpiEnd)

Same as mpiFindIndicesSpecial, but does not use myid, instead allows it to be passed through parameter i.

Parameters

N – data size for which data partitioning scheme needs to be found
i – parameters for processor i that is needed
mpiStart – start index of processor i in the global dataset
mpiEnd – end index of processor i in the global dataset

void mpiFindIndices(long N, int blkSize, int i, int *mpiStart, int *mpiEnd)

Function to find start and end indices for the given processor i. If blkSize is 1, N is divided equally among processors. If blkSize > 0, N*blkSize stars are divided among procs in such a way that each contain multiples of blkSize.

Parameters

N – the number of blocks of size blkSize to be divided among processors
blkSize – size of blocks
i – processor id
mpiStart – indices
mpiEnd –

void mpiFindIndicesCustom(long N, int blkSize, int i, int *mpiStart, int *mpiEnd)

Populates the mpiStart and mpiEnd values given the initial data size N, and the factor blkSize which data in each processor is required to be a multiple of (except the last processor). Followin is the data partitioning scheme - Each processor gets data that is a multiple of blkSize, and the remainder after division goes to the last processor. For more details, refer to: http://adsabs.harvard.edu/abs/2013ApJS..204…15P.

Parameters

N – data size that needs to be partitioned
blkSize – the factor blkSize which data in each processor is required to be a multiple of (except the last one).
mpiStart – start index of the data elements in the global array belonging to processor i.
mpiEnd – array end index of the data elements in the global array belonging to processor i.

void mpiFindDispAndLenCustom(long N, int blkSize, int *mpiDisp, int *mpiLen)

finds displacement and count for all processors

Parameters

N – data size
blkSize – block size which each chunk has to be a multiple of
mpiDisp – displacement array
mpiLen – count array

MPI_Comm inv_comm_create(int procs, MPI_Comm old_comm)

Creates a new communicator with inverse order of processes.

Parameters

procs – number of processors
old_comm – old communicator

Returns

new communicator with inverse order of ranks

cmc_nr.c

Defines

FUNC(j, k, E, J)

Parallel version of FUNC macro, requires global indices as input after index transformation.

Parameters

j – ?
k – ?
E – energy
J – angular momentum

Returns

?

Functions

long FindZero_r_serial(long kmin, long kmax, double r)

binary search on r (pure serial version): given the array star[].r and the two indices kmin and kmax, with the conditions 1) array is monotonic in its indices, 2) kmin<kmax, find and return the index k, such that star[k].r<r<star[k+1].r

Parameters

kmin – min start index for bisection
kmax – max start index for bisection
r – target value

Returns

index k, such that star[k].r<r<star[k+1].r

long FindZero_r(long kmin, long kmax, double r)

binary search on r: given the array star[].r and the two indices kmin and kmax, with the conditions 1) array is monotonic in its indices, 2) kmin<kmax, find and return the index k, such that star[k].r<r<star[k+1].r

Parameters

kmin – min start index for bisection
kmax – max start index for bisection
r – target value

Returns

index k, such that star[k].r<r<star[k+1].r

double sigma_r(double r)

binary search on sigma_array.r given the array sigma_array.r[] and the two indices kmin and kmax, with the conditions 1) array is monotonic in its indices, 2) kmin<kmax, find the index k, such that sigma_array.r[k]<r<sigma_array.r[k+1]

Parameters: r – target value
Returns: sigma value at index k, such that sigma_array.r[k]<r<sigma_array.r[k+1]

long FindZero_Q(long j, long kmin, long kmax, double E, double J)

another binary search, except FUNC(k) may be decreasing rather than increasing

Parameters

j – star index
kmin – min index for bisection
kmax – max index for bisection
E – energy
J – angular momentum

Returns

index k, such that sigma_array.r[k]<r<sigma_array.r[k+1]

double calc_pot_within_interval(double r, void *p)

?

Parameters

r – ?
p – ?

Returns

?

double calc_pot_in_interval(double r, long k)

?

Parameters

r – ?
k – ?

Returns

?

long find_zero_Q(long j, long kmin, long kmax, long double E, long double J)

another binary search, except FUNC(k) may be decreasing rather than increasing

Parameters

j – star index
kmin – min index for bisection
kmax – max index for bisection
E – energy
J – angular momentum

Returns

index k, such that sigma_array.r[k]<r<sigma_array.r[k+1]

double calc_vr_within_interval(double r, void *p)

Calculate the square of vr !

Parameters

r – ?
p – ?

Returns

?

double calc_vr_in_interval(double r, long index, long k, double E, double J)

Calculates the square of vr.

Parameters

r – ?
index – ?
k – ?
E – energy
J – angular momentum

Returns

?

double calc_Q_within_interval(double r, void *p)

?

Parameters

r – ?
p – ?

Returns

?

double calc_vr(double r, long index, double E, double J)

Calculates the square of vr!

Parameters

r – radial position
index – star index
E – energy
J – angular momentum

Returns

square of vr

double find_root_vr(long index, long k, double E, double J)

Calculates the root of vr with the additional constraint that when substituting it back into calc_vr_within_interval leads to vr>= 0.0 .

Parameters

index – ?
k – ?
E – energy
J – angular momentum

Returns

?

struct calc_vr_params

Public Members

long index

long k

double E

double J

cmc_orbit.c

Functions

long get_positive_Q_index(long index, double E, double J)

?

Parameters

index – ?
E – energy
J – angular momentum

Returns

?

long symmetric_positive_Q_search(long index, long kmid, double E, double J, long kmin, long kmax)

Calculates Q at star position kmid and follows N_Q_TRACE Q-evaluations symmetrically to both sides of kmid. If there are no positive Q-values and Q has been decreasing monotonously a circular orbit is returned (-1). If there was a positive Q-value its index k is returned. Otherwise the next N_Q_TRACE points are evaluated. If the loop ends without result -2 is returned.

Parameters

index – ?
kmid – ?
E – ?
J – ?
kmin – ?
kmax – ?

Returns

?

long get_positive_Q_index_new_J(long index, double E, double J, orbit_rs_t orbit_old)

?

Parameters

index – ?
E – ?
J – ?
orbit_old – ?

Returns

?

struct star_coords get_position(long index, double E, double J, double old_r, orbit_rs_t orbit)

?

Parameters

index – ?
E – ?
J – ?
old_r – ?
orbit – ?

Returns

?

void debug_msg_orbit_new_inside_sqrt(double actual_rmin, double rmin, double inside_sqrt, long k, int ismin)

?

Parameters

actual_rmin – ?
rmin – ?
inside_sqrt – ?
k – ?
ismin – ?

void debug_msg_supp_mpi_info(char *msg, int proc, long g_si, long si, int mpiBegin, int mpiEnd)

Print supplementary information related to MPI parallelization.

Parameters

msg – Message context string
proc – Current processor
g_si – Global particle index
si – Local particle index
mpiBegin – Start of memory partition in global address space
mpiEnd – End of memory partition in global address space

void debug_msg_orbit_new_outside_interval(long si, double rmin, double rmin_new, long k, double rk, double rk1, int ismin, double E, double J)

?

Parameters

si – ?
rmin – ?
rmin_new – ?
k – ?
rk – ?
rk1 – ?
ismin – ?
E – ?
J – ?

orbit_rs_t calc_orbit_new(long index, double E, double J)

finds out new orbit of star, primarily peri and apo- center distances

Parameters

index – star index
E – energy
J – angular momentum

Returns

orbit structure containing properties of new orbit

orbit_rs_t calc_orbit_new_J(long index, double J, struct star_coords pos_old, orbit_rs_t orbit_old)

?

Parameters

index – ?
J – ?
pos_old – ?
orbit_old – ?

Returns

?

void set_a_b(long index, long k, double *a, double *b)

?

Parameters

index – ?
k – ?
a – ?
b – ?

long find_zero_Q_slope(long index, long k, double E, double J, int positive)

?

Parameters

index – ?
k – ?
E – ?
J – ?
positive – ?

Returns

?

cmc_relaxation.c

Functions

struct star_coords get_star_coords_from_star(long index)

These routines are for calculating the orbital perturbations of a pair of stars due to relaxation. They still depend on the global ‘star’ array and are thus not fully modular.

Parameters: index – index of star
Returns: ?

void set_star_coords_for_star(struct star_coords pos)

?

Parameters: pos – ?

struct encounter get_encounter_dyns(struct star_coords star1, struct star_coords star2, gsl_rng *rng)

?

Parameters

star1 – ?
star2 – ?
rng – ?

Returns

?

struct relaxation_params get_relaxation_params(double dt, struct encounter enc)

?

Parameters

dt – ?
enc – ?

Returns

?

double get_relaxation_time(struct encounter enc)

?

Parameters: enc – ?
Returns: ?

double scattering_angle(double dt, double Trel)

?

Parameters

dt – ?
Trel – ?

Returns

?

struct perturbation scatter_relax(struct encounter enc, struct relaxation_params rparams)

?

Parameters

enc – ?
rparams – ?

Returns

?

struct perturbation scatter_relax_old(struct star_coords pos[2], double dt, gsl_rng *rng)

?

Parameters

pos[2] – ?
dt – ?
rng – ?

Returns

?

cmc_remove_star.c

Functions

int remove_old_star(double time, long k)

Sourav: checks whether a star should die or not before performing dynamics.

Parameters

time – ?
k – index of star

Returns

?

cmc_search_grid.c

Functions

struct Search_Grid *search_grid_initialize(double power_law_exponent, double fraction, long starsPerBin, double part_frac)

?

Parameters

power_law_exponent – ?
fraction – ?
starsPerBin – ?
part_frac – ?

Returns

?

double search_grid_estimate_prop_const(struct Search_Grid *grid)

?

Parameters: grid – ?
Returns: ?

void search_grid_allocate(struct Search_Grid *grid)

?

Parameters: grid – ?

void search_grid_update(struct Search_Grid *grid)

?

Parameters: grid – ?

double search_grid_get_r(struct Search_Grid *grid, long index)

This function gives only the approximate radius of the search grid at index “index” and cannot be used to infer the right search interval, due to round-off errors.

Parameters

grid – ?
index – ?

Returns

?

struct Interval search_grid_get_interval(struct Search_Grid *grid, double r)

Returns the interval of particle indices such that star[min].r>r>star[max].r.

Parameters

grid – ?
r – ?

Returns

?

long search_grid_get_grid_index(struct Search_Grid *grid, double r)

?

Parameters

grid – ?
r – ?

Returns

?

double search_grid_get_grid_indexf(struct Search_Grid *grid, double r)

?

Parameters

grid – ?
r – ?

Returns

?

void search_grid_free(struct Search_Grid *grid)

?

Parameters: grid – ?

void search_grid_print_binsizes(struct Search_Grid *grid)

?

Parameters: grid – ?

cmc_sort.c

Defines

GENSORT_NAME

GENSORT_TYPE

GENSORT_KEYTYPE

GENSORT_GETKEY(a)

GENSORT_COMPAREKEYS(k1, k2)

GENSORT_USEPOINTERS

Functions

keyType getKey(type *buf)

returns the key of a given data element

Parameters: buf – input data element
Returns: key of the data element

int compare_keyType(const void *a, const void *b)

comparison function for sort.

Parameters

a – first key
b – second key

Returns

returns 1 if a > b, -1 if a < b, 0 if equal.

int compare_type(const void *a, const void *b)

comparison function for sort

Parameters

a – first data element
b – second data element

Returns

returns 1 if a > b, -1 if a < b, 0 if equal.

void find_expected_count(int *expected_count, int N, int numprocs)

given the total number of data elements, number or processors, computes the expected count based on the data partitioning scheme

Parameters

expected_count – array where the expected count will be store after computation
N – total number of data elements
numprocs – number of processors

keyType *sample(type *buf, keyType *sample_array, int N, int n_samples)

returns n_samples samples from the given N data elements in array buf

Parameters

buf – array that needs to be sampled
sample_array – array in which samples will be stored
N – number of data elements in buf
n_samples – number of samples required

Returns

int binary_search(type *buf, keyType r, int kmin, int kmax)

binary search on buf for r

Parameters

buf – input array on which search is done
r – target value
kmin – left start index
kmax – right start index

Returns

index k such that buf[k] < r < buf[k+1]

void remove_stripped_stars(type *buf, int *local_N)

Stripped stars have their r values set to SF_INFINITY. The first step of the parallel sort is a local sort by each processor. During this, these stars are pushed to the end of the local array. This routine removes all of these stars starting from the end of the sorted array.

Parameters

buf – array sorted by r
local_N – number of stars in the array

int sample_sort(type *buf, int *local_N, MPI_Datatype dataType, binary_t *b_buf, MPI_Datatype b_dataType, MPI_Comm commgroup, int n_samples)

Parallel sample sort. Following are the steps:

Local sort - each processor sorts the local data using sequential sort
Sampling - each processor samples s keys from the local dataset and sends them to the root node. The root node sorts these s*p keys, and picks p regularly sampled keys/points among these, called splitter points which are broadcasted to all processors. These points are supposed to divide the entire dataset into p buckets with approximately equaly number of points in them. This extent to which this division will be equal is a factor of the number of samples s contributed by each processor. If s=p, and regular sampling is used, the load balance will be ideal i.e. the maximum number of points ending up in a processor would be at most ~ 2*N/p. If s < p, it’ll be worse, and is s > p, it’ll be better.
Each processor places these splitter points into their sorted local array using binary search.
All processors have an all-to-all communication and exchange data
Each processor sorts the received chunks of data which completes the sort
Since the number of points ending up on each processor is non-deterministic, an optional phase is to exchange data between processors so that the number of data points on each processor is in accordance with out data partitioning scheme.

Parameters

buf – the local data set (star) which is a part of the entire data set which is divided among many processors which is to be sorted in parallel
local_N – number of local data points
dataType – MPI datatype for the star data structure
b_buf – binary local data set
b_dataType – MPI datatype for the binary data structure
commgroup – MPI communication group
n_samples – number of samples to be contributed by each processor

Returns

total number of stars that were sorted

void load_balance(type *inbuf, type *outbuf, binary_t *b_inbuf, binary_t *b_outbuf, int *expected_count, int *actual_count, int myid, int procs, MPI_Datatype dataType, MPI_Datatype b_dataType, MPI_Comm commgroup): after sample sort, exchanges data between processors to make sure the number of stars in each processor is in accordance with the data partitioning scheme

cmc_sscollision.c

Functions

void sscollision_do(long k, long kp, double rperimax, double w[4], double W, double rcm, double vcm[4], gsl_rng *rng)

Does single single collision.

Parameters

k – index of first star
kp – index of second star
rperimax – ?
w[4] – ?
W – ?
rcm – ?
vcm[4] – ?
rng – gsl rng

void merge_two_stars(star_t *star1, star_t *star2, star_t *merged_star, double *vs, struct rng_t113_state *s)

merge two stars using stellar evolution if it’s enabled

Parameters

star1 – pointer to star 1
star2 – pointer to star 2
merged_star – pointer to merged star
vs – ?
s – rng state

double coll_CE(double Mrg, double Mint, double Mwd, double Rrg, double vinf)

calculate resulting semi-major axis from collisional common envelope event

Parameters

Mrg – mass of red giant
Mint – mass of intruder (NS, BH, MS, etc.)
Mwd – mass of RG core that will become WD
Rrg – radius of RG
vinf – relative velocity at infinity between RG and intruder

Returns

resulting semi-major axis from collisional common envelope event

double coll_CE_twogiant(double M1, double M2, double Mc1, double Mc2, double R1, double R2, double vinf)

calculate resulting semi-major axis from collisional common envelope event for two giant stars

Parameters

M1 – mass of giant
M2 – mass of giant
Mc1 – mass of giant core that will become WD
Mc2 – mass of giant core that will become WD
R1 – radius of giant
R2 – radius of giant
vinf – relative velocity at infinity between the two giant stars

Returns

resulting semi-major axis from collisional common envelope event

double r_of_m(double M)

the mass—radius relationship, from Freitag, et al. (2004) M is expected to be input in the same units as star[].m radius is output in code (N-body) units

This is ONLY used if stellar evolution is off (for a toy rejuvination problem)

Parameters: M – ?
Returns: radius in code (N-body) units

cmc_stellar_evolution.c

Functions

void restart_stellar_evolution(void): This is the same as the stellar_evolution_init, except we’re only setting the global BSE properties (the individual star and binary properties were set in the restart binary files)

void stellar_evolution_init(void): Initializes stellar evolution, both the global BSE properties and the individual properties for each star. Runs BSE for a dt->0 time to set a bunch of the parameters for individual stars.

void do_stellar_evolution(gsl_rng *rng)

does stellar evolution using sse and bse packages.

Parameters: rng – gsl rng

void write_stellar_data(void): ?

void handle_bse_outcome(long k, long kb, double *vs, double tphysf, int kprev0, int kprev1)

?

Parameters

k – index of star 1
kb – index of star 2
vs – ?
tphysf – ?

void pulsar_write(long k, double kick)

Output info of pulsars.

Parameters

k – star index
kick – ?

void write_morepulsar(long i)

Output info of morepulsars.

Parameters: k – star index

void getCMCvalues(long s_id, double kick, int remtype1, int remtype2, int progtype1, int progtype2, double formation, int boomtype, int boomstar)

?

Parameters

s_id – ?
kick – ?
remtype1 – ?
remtype2 – ?
progtype1 – ?
progtype2 – ?
formation – ?
boomtype – ?
boomstar – ?

void bh_count(long k)

Meagan: count different types of bh-objects; at end of timestep, we’ll print these totals.

Parameters: k – index of star

void cp_binmemb_to_star(long k, int kbi, long knew)

copies corresponding binary member variables of given star to the star variables of new star

Parameters

k – star index
kbi – 0 or 1 based on which binary member
knew – index of new star

void cp_SEvars_to_newstar(long oldk, int kbi, long knew)

copies variables of old star to stellar evolution variables of new star

Parameters

oldk – old star index
kbi – 0 or 1 for binary, -1 for non-binary
knew – index of new star

void cp_m_to_newstar(long oldk, int kbi, long knew)

copies mass from old to new star

Parameters

oldk – old star index
kbi – 0 or 1 for binar, -1 for non-binary
knew – index of new star

void cp_SEvars_to_star(long oldk, int kbi, star_t *target_star)

same as cpSEvars_to_newstar, only difference in function signature

Parameters

oldk – old star index
kbi – 0 or 1 for binary, -1 for non-binary
target_star – pointer to target star

void cp_m_to_star(long oldk, int kbi, star_t *target_star)

same as cp_m_to_newstar, only difference in function signature

Parameters

oldk – old star index
kbi – 0 or 1 for binary, -1 for non-binary
target_star – pointer to target star

void cp_SEvars_to_newbinary(long oldk, int oldkbi, long knew, int kbinew)

copies stellar evoltion variables to new binary (set everything except tb)

Parameters

oldk – old star index
oldkbi – 0 or 1 for binary, -1 for non-binary
knew – index of new star
kbinew – ?

void cp_starSEvars_to_binmember(star_t instar, long binindex, int bid)

?

Parameters

instar – ?
binindex – ?
bid – ?

void cp_starmass_to_binmember(star_t instar, long binindex, int bid)

?

Parameters

instar – ?
binindex – ?
bid – ?

double peters_t_insp_integral(double e, void *params)

double peters_t_insp_eccentric(double a, double e, double beta)

double peters_t_insp(double mG3c5, double a, double e)

int rhs_peters(double t, const double y[], double f[], void *params)

void integrate_a_e_peters_eqn(long kb)

cmc_trace.c

Warning

doxygenfile: Cannot find file “cmc_trace.c

cmc_utils.c

Defines

GENSEARCH_NAME

GENSEARCH_TYPE

GENSEARCH_KEYTYPE

GENSEARCH_GETKEY(a)

GENSEARCH_COMPAREKEYS(k1, k2)

NR_END

FREE_ARG

Functions

double sqr(double x)

a fast square function

Parameters: x – number to be squared
Returns: square of x

double cub(double x)

a fast cube function

Parameters: x – number to be cubed
Returns: cube of x

double potential_serial(double r)

The potential computed using the star[].phi computed at the star locations in star[].r sorted by increasing r. Note: this is the pure serial version.

Parameters: r – position at which potential is required
Returns: potential at r

double potential(double r)

The potential computed using the phi computed at the star locations in r, sorted by increasing r.

Parameters: r – position at which potential is required
Returns: potential at position r

void toggle_debugging(int signal)

toggle debugging

Parameters: signal – signal

void exit_cleanly_old(int signal)

old exit function

Parameters: signal – signal

void exit_cleanly(int signal, const char *fn)

close buffers, then exit

Parameters

signal – signal
fn – function which called exit

void free_arrays(void): free’s dynamically allocated arrays

void sf_gsl_errhandler(const char *reason, const char *file, int line, int gsl_errno)

GSL error handler.

Parameters

reason – string containing reason for error
file – filename
line – line number
gsl_errno – gsl returned error number

void set_velocities3(void): set velocities a la Stodolkiewicz to be able to conserve energy (parallel version of set_velocities3)

void energy_conservation3(void): Wrapper function for energy conservation scheme.

void ComputeIntermediateEnergy(void): computes intermediate energies, and transfers “new” dynamical params to the standard variables

void update_tspent(struct tms tmsbufref)

long DynamicalFrictionTimescale()

compute the dynamical friciton timescale

Param

long CheckStop()

makes a few checks at the beginning of each timestep to make sure the simulation is proceeding normally, and expected on some abnormal activity, or if simulation reaches some of the user-specified termination conditions.

Parameters: tmsbufref –
Returns: 0 if error and simulation needs to be terminated, 1 otherwise

int CheckCheckpoint()

void ComputeEnergy(void): Calculates E,J for every star. Also, calculates, global energy variabies (parallel version of ComputeEnergy)

long potential_calculate(void)

Parallel version of potential_calculate(). Currently the entire calculation is done by all nodes since it uses only the duplicated arrays r and m, but might consider parallelization in the long term to improve scalability. Computing the potential at each star sorted by increasing radius. Units: G = 1 and Mass is in units of total INITIAL mass. Total mass is computed by SUMMING over all stars that have NOT ESCAPED i.e., over all stars upto N_MAX <= N_STAR. N_MAX is computed in this routine by counting all stars with radius < SF_INFINITY and Radius of the (N_MAX+1)th star is set to infinity i.e., star[N_MAX+1].r = SF_INFINITY. Also setting star[N_MAX+1].phi = 0. Assuming star[0].r = 0. star[].phi is also indexed i.e. it is the value of the potential at radius star[k].r NOTE: Assming here that NO two stars are at the SAME RADIUS upto double precision. Returns N_MAX. Potential given in star[].phi.

Returns: total number of stars

int find_stars_mass_bin(double smass)

find the star[i]’s mass bin return -1 on failure

Parameters: smass – ?
Returns: ?

void comp_multi_mass_percent(): computes the Lagrange radii for various mass bins stored in the array mass_bins[NO_MASS_BINS]

void comp_mass_percent(): Computes radii containing mass_pc[] % of the mass.

double fastpotential(double r, long kmin, long kmax)

computes potential at position r

Parameters

r – position at which potential is required
kmin – kmin for bisection
kmax – kmax for bisection

Returns

potential at r

long check_if_r_around_last_index(long last_index, double r)

?

Parameters

last_index – ?
r – ?

Returns

?

void nrerror(char error_text[])

Numerical Recipes standard error handler.

Parameters: error_text[] – error text

double *vector(long nl, long nh)

allocate a double vector with subscript range v[nl..nh]

Parameters

nl – lower bound for subscript range
nh – upper bound for subscript range

Returns

?

int *ivector(long nl, long nh)

allocate an int vector with subscript range v[nl..nh]

Parameters

nl – lower bound for subscript range
nh – upper bound for subscript range

Returns

?

void free_vector(double *v, long nl, long nh)

free a double vector allocated with vector()

Parameters

v – double vector pointer
nl – lower bound for subscript range
nh – upper bound for subscript range

void free_ivector(int *v, long nl, long nh)

free an int vector allocated with ivector()

Parameters

v – int vector pointer
nl – lower bound for subscript range
nh – upper bound for subscript range

void update_vars(void): update some important global variables - total number, mass, and binding energy of binaries in cluster

void units_set(void): set the units

double compute_tidal_boundary(void): compute the tidal boundary given the current time of the cluster and the dominant component of the tidal tensor

void central_calculate(void): calculate central quantities: see description in Fregeau & Rasio (2007) based upon Casertano & Hut (1985)

double local_kT(long si, int p)

?

Parameters

si – ?
p – ?

Returns

?

double core_kT(int Ncore, int p)

?

Parameters

Ncore – ?
p – ?

Returns

?

central_t central_hard_binary(double ktmin, central_t old_cent)

?

Parameters

ktmin – ?
old_cent – ?

Returns

?

void clusdyn_calculate(void): Parallel version of clusdyn_calculate. Done by all nodes since accessing only duplicate arrays.

double function_q(long j, long double r, long double pot, long double E, long double J)

calculates function Q (or vr^2)

This is probably the nth incarnation of vr^2 (aka Q). It differs in as much as it is not a macro but a real function with input parameters all being long double. This has the effect that, on IA-32 systems, vr^2 has always a predictable precision, even when the 387 fpu is used. The reason is simply that vr^2 is now forced to be calculated with extended double precision matching the precision that is used internally by the 387 fpu. This way, possibly mixed type (double, extended) expressions are avoided when 387 instructions are generated by the compiler (This problem does not exist with -msse2 -mfpmath=sse). Previously, the results depended on the order of evaluation of the subexpressions, which changes with the optimization level and surrounding code, sometimes varying by much more than DBL_EPSILON. It is clear that there is little one can do about this problem as long as vr^2 is a macro, so I decided to write vr^2 as a function with all variables declared long double. I also declared it as inlined (see cmc.h), so there should be no performance impact. — Stefan, 10/01/07

Parameters

j – star index
r – position
pot – potential
E – energy
J – angular momentum

Returns

vr^2 (aka Q)

double timeStartSimple()

starts timer - returns double precision time when this function is called

Returns: double precision time in microseconds when this function was called

void timeEndSimple(double timeStart, double *timeAccum)

ends timer - computes difference between current time and input time

Parameters

timeStart – start time which needs to be subtracted from current time
timeAccum – variable to store difference between current time and start time

void timeStart2(double *st)

void timeEnd2(char *fileName, char *funcName, double *st, double *end, double *tot)

void create_timing_files(): creates timing files

void set_global_vars1()

sets_global_variable_for_tracking_events

sets some important global variables

void set_global_vars2()

void calc_sigma_new(): wrapper function for calc_sigma_r

void bin_vars_calculate(): Calculates some global binary variables - total binary mass,and E.

void calc_potential_new(): Wrapper function for potential calculation.

void compute_energy_new(): wrapper function for ComputeEnergy

void set_bh_counters(): set some of the BH counters to zero

void set_energy_vars(): initializing some energy variables

void reset_interaction_flags(): reset interaction flags

void calc_clusdyn_new(): wrapper function for clusdyn_calculate()

void calc_timestep(gsl_rng *rng)

calculates timestep

Parameters: rng – gsl rng

void energy_conservation1(): computes some numbers necessary to implement Stodolkiewicz’s energy conservation scheme

void toy_rejuvenation(): Sourav: checking all stars for their possible extinction from old age Sourav: toy rejuvenation: DMrejuv storing amount of mass loss per time step.

void new_orbits_calculate(): wrapper function for get_positions (first part of older tidally_strip_stars)

void energy_conservation2(): more numbers necessary to implement Stodolkiewicz’s energy conservation scheme

void qsorts_new(void): Sort stars based on their radial positions. In the serial version, simple quick sort is used. In the parallel version, parallel sample sort is used.

void post_sort_comm(): The duplicated arrays have to be collected and synchronized after the sorting. This is where this is done.

void findIndices(long N, int blkSize, int i, int *begin, int *end)

Given a number of blocks N each of size blkSize, and processor id i, returns the begin and end indices of an array of size N*blkSize that the ith processor will get if the data was partitioned in the following way.

Parameters

N – number of blocks of size blkSize that needs to be partitioned among processors.
blkSize – the factor blkSize which data in each processor is required to be a multiple of (except the last one).
i – processor id
begin – the index of the first data element that will go to proc i
end – the index of the last data element that will go to proc i

void findLimits(long N, int blkSize)

Populates the Start and End arrays given the initial data size N, and the factor blkSize which data in each processor is required to be a multiple of (except the last processor). Followin is the data partitioning scheme - Each processor gets data that is a multiple of blkSize, and the remainder after division goes to the last processor. For more details, refer to: http://adsabs.harvard.edu/abs/2013ApJS..204…15P.

Parameters

N – data size that needs to be partitioned
blkSize – the factor blkSize which data in each processor is required to be a multiple of (except the last one).

int findProcForIndex(int j)

For the serial version to mimic the parallel, random numbers have to be generated from the right stream. For this purpose, we need a mapping between the star index in the array and the processor they will belong to in a corresponding parallel run. This function provides this mapping - for a given index j, it returns its processor id.

Parameters: j – star id for which we want to know the processor id
Returns: processor id that the star would belong to in a corresponding parallel run

void set_rng_states(): Initializes the rng related variables, and also sets the starting states based on the initial seed for both serial and parallel versions.

int get_global_idx(int i)

Function which performs the index transformation from the given local index of a particular processor to the global index based on the data partitioning scheme. If the given index is less than the allotted number of stars a processor is supposed to hold, the corresponding global value is returned. If it is beyond, that means it is a newly created star. In this case, an appropriate value beyond N_MAX is returned since the global values of the newly created stars are stored beyond N_MAX appropriately.

Parameters: i – local index for which we need the global index
Returns: global index corresponding to the local index

int get_local_idx(int i)

This does the opposite of get_global_idx, and performs the index transformation from the given global index to the local index of a particular processor based on the data partitioning scheme.

Parameters: i – global index for which we need the local index
Returns: local index corresponding to the global index

FEWBODY C Code

binbin.c

Functions

void print_usage(FILE *stream)

int calc_units(fb_obj_t *obj[2], fb_units_t *units)

int main(int argc, char *argv[])

binsingle.c

Functions

void print_usage(FILE *stream)

int calc_units(fb_obj_t *obj[2], fb_units_t *units)

int main(int argc, char *argv[])

cluster.c

Functions

void print_usage(FILE *stream)

void calc_units(fb_hier_t hier, fb_units_t *units)

double fv(double v, void *params)

double vf(double f)

int main(int argc, char *argv[])

fewbody.c

Functions

fb_ret_t fewbody(fb_input_t input, fb_units_t units, fb_hier_t *hier, double *t, gsl_rng *rng, struct rng_t113_state *curr_st)

Variables

int fb_debug = 0

fewbody_classify.c

Functions

int fb_classify(fb_hier_t *hier, double t, double tidaltol, double speedtol, fb_units_t units, fb_input_t input)

int fb_is_stable(fb_obj_t *obj, double speedtol, fb_units_t units, fb_input_t input)

int fb_is_stable_binary(fb_obj_t *obj, double speedtol, fb_units_t units, fb_input_t input)

int fb_is_stable_triple(fb_obj_t *obj)

int fb_is_stable_quad(fb_obj_t *obj)

int fb_mardling(fb_obj_t *obj, int ib, int is)

fewbody_coll.c

Functions

int fb_is_collision(double r, double R1, double R2, double M1, double M2, double k1, double k2, double mass_units, double length_units, double bhns_tde_flag)

int fb_collide(fb_hier_t *hier, double f_exp, fb_units_t units, gsl_rng *rng, struct rng_t113_state *curr_st, double bh_reff, fb_input_t input)

void fb_merge(fb_obj_t *obj1, fb_obj_t *obj2, int nstarinit, double f_exp, fb_units_t units, gsl_rng *rng, struct rng_t113_state *curr_st, double bh_reff)

void fb_bh_merger(double m1, double m2, double a1, double a2, double *mass_frac, double *afinal, double *v_para, double *v_perp, struct rng_t113_state *curr_st)

fewbody_hier.c

Functions

void fb_malloc_hier(fb_hier_t *hier)

void fb_init_hier(fb_hier_t *hier)

void fb_free_hier(fb_hier_t hier)

void fb_trickle(fb_hier_t *hier, double t)

void fb_elkcirt(fb_hier_t *hier, double t)

int fb_create_indices(int *hi, int nstar)

int fb_n_hier(fb_obj_t *obj)

char *fb_sprint_hier(fb_hier_t hier, char string[FB_MAX_STRING_LENGTH])

char *fb_sprint_hier_hr(fb_hier_t hier, char string[FB_MAX_STRING_LENGTH])

void fb_upsync(fb_obj_t *obj, double t)

void fb_randorient(fb_obj_t *obj, gsl_rng *rng, struct rng_t113_state *curr_st)

void fb_downsync(fb_obj_t *obj, double t)

void fb_objcpy(fb_obj_t *obj1, fb_obj_t *obj2)

fewbody_int.c

Functions

void fb_malloc_ks_params(fb_ks_params_t *ks_params)

void fb_init_ks_params(fb_ks_params_t *ks_params, fb_hier_t hier)

void fb_free_ks_params(fb_ks_params_t ks_params)

void fb_malloc_nonks_params(fb_nonks_params_t *nonks_params)

void fb_init_nonks_params(fb_nonks_params_t *nonks_params, fb_hier_t hier)

void fb_free_nonks_params(fb_nonks_params_t nonks_params)

fewbody_io.c

Functions

void fb_print_version(FILE *stream)

void fb_print_story(fb_obj_t *star, int nstar, double t, char *logentry)

fewbody_isolate.c

Functions

int fb_collapse(fb_hier_t *hier, double t, double tidaltol, double speedtol, fb_units_t units, fb_nonks_params_t nonks_params, fb_input_t input)

int fb_expand(fb_hier_t *hier, double t, double tidaltol)

fewbody_ks.c

Functions

double fb_ks_dot(double x[4], double y[4])

double fb_ks_mod(double x[4])

void fb_calc_Q(double q[4], double Q[4])

void fb_calc_ksmat(double Q[4], double Qmat[4][4])

void fb_calc_amat(double **a, int nstar, int kstar)

void fb_calc_Tmat(double **a, double *m, double **T, int nstar, int kstar)

int fb_ks_func(double s, const double *y, double *f, void *params)

double fb_ks_Einit(const double *y, fb_ks_params_t params)

void fb_euclidean_to_ks(fb_obj_t **star, double *y, int nstar, int kstar)

void fb_ks_to_euclidean(double *y, fb_obj_t **star, int nstar, int kstar)

fewbody_nonks.c

Defines

FB_FM(i, j, k)

FB_REL(i, j, k)

pi2

Functions

int fb_nonks_func(double t, const double *y, double *f, void *params)

int fb_nonks_jac(double t, const double *y, double *dfdy, double *dfdt, void *params)

void fb_euclidean_to_nonks(fb_obj_t **star, double *y, int nstar)

void fb_nonks_to_euclidean(double *y, fb_obj_t **star, int nstar)

fewbody_scat.c

Functions

void fb_init_scattering(fb_obj_t *obj[2], double vinf, double b, double rtid)

void fb_normalize(fb_hier_t *hier, fb_units_t units)

fewbody_utils.c

Functions

double *fb_malloc_vector(int n)

double **fb_malloc_matrix(int nr, int nc)

void fb_free_vector(double *v)

void fb_free_matrix(double **m)

double fb_sqr(double x)

double fb_cub(double x)

double fb_dot(double x[3], double y[3])

double fb_mod(double x[3])

int fb_cross(double x[3], double y[3], double z[3])

int fb_angmom(fb_obj_t *star, int nstar, double L[3])

void fb_angmomint(fb_obj_t *star, int nstar, double L[3])

double fb_einttot(fb_obj_t *star, int nstar)

double fb_einttot_rel(fb_obj_t *star, int nstar)

double fb_petot(fb_obj_t *star, int nstar)

double fb_ketot(fb_obj_t *star, int nstar)

double fb_outerpetot(fb_obj_t **obj, int nobj)

double fb_dedt_gw(fb_obj_t *star, int nstar, fb_units_t units, fb_nonks_params_t nonks_params)

void fb_n_ecc(fb_obj_t *obj1, fb_obj_t *obj2, double *a, double *e, fb_units_t units)

void fb_check_ecc_for_inspiral(fb_obj_t *obj1, fb_obj_t *obj2, double sep_M, fb_units_t units)

double fb_pn_ecc_t(fb_obj_t *star1, fb_obj_t *star2, fb_units_t units, fb_nonks_params_t nonks_params)

double fb_E_rel(fb_obj_t *star1, fb_obj_t *star2, fb_units_t units, fb_nonks_params_t nonks_params)

double E_rel(fb_obj_t *star1, fb_obj_t *star2, fb_units_t units, fb_nonks_params_t nonks_params)

double J_rel(fb_obj_t *star1, fb_obj_t *star2, fb_units_t units, fb_nonks_params_t nonks_params)

double fb_compute_distance_in_M(fb_obj_t *star1, fb_obj_t *star2, fb_units_t units)

double fb_Etot_rel(fb_obj_t *star, int nstar, fb_units_t units, fb_nonks_params_t nonks_params)

double fb_outerketot(fb_obj_t **obj, int nobj)

double fb_kepler(double e, double mean_anom)

double fb_keplerfunc(double ecc_anom, void *params)

double fb_reltide(fb_obj_t *bin, fb_obj_t *single, double r)

sigma_binsingle.c

Functions

void print_usage(FILE *stream)

int calc_units(fb_obj_t *obj[2], fb_units_t *units)

int main(int argc, char *argv[])

triple.c

Functions

void print_usage(FILE *stream)

int calc_units(fb_obj_t *obj[1], fb_units_t *units)

int main(int argc, char *argv[])

triplebin.c

Functions

void print_usage(FILE *stream)

int calc_units(fb_obj_t *obj[2], fb_units_t *units)

int main(int argc, char *argv[])

BSE C CODE

Functions

void bse_zcnsts(double *z, double *zpars)

calculate metallicity constants

Parameters

z – ?
zpars – ?

void bse_evolve_single(int *kw, double *mass, double *mt, double *r, double *lum, double *mc, double *rc, double *menv, double *renv, double *ospin, double *epoch, double *tms, double *tphys, double *tphysf, double *dtp, double *z, double *zpars, double *vs, double *bhspin)

void bse_evolv2(int *kstar, double *mass0, double *mass, double *rad, double *lum, double *massc, double *radc, double *menv, double *renv, double *ospin, double *B_0, double *bacc, double *tacc, double *epoch, double *tms, double *tphys, double *tphysf, double *dtp, double *z, double *zpars, double *tb, double *ecc, double *vs, double *bhspin)

evolve a binary

Parameters

kstar – ?
mass0 – ?
mass – ?
rad – ?
lum – ?
massc – ?
radc – ?
menv – ?
renv – ?
ospin – ?
B_0 – ?
bacc – ?
tacc – ?
epoch – ?
tms – ?
tphys – ?
tphysf – ?
dtp – ?
z – ?
zpars – ?
tb – ?
ecc – ?
vs – ?

void bse_evolv2_safely(int *kstar, double *mass0, double *mass, double *rad, double *lum, double *massc, double *radc, double *menv, double *renv, double *ospin, double *B_0, double *bacc, double *tacc, double *epoch, double *tms, double *tphys, double *tphysf, double *dtp, double *z, double *zpars, double *tb, double *ecc, double *vs, double *bhspin)

evolve a binary star safely: in some cases, a merger can have non self-consistent properties, leading to crazy things like NaN radii—this is the easiest way to get around that problem

Parameters

kstar – ?
mass0 – ?
mass – ?
rad – ?
lum – ?
massc – ?
radc – ?
menv – ?
renv – ?
ospin – ?
B_0 – ?
bacc – ?
tacc – ?
epoch – ?
tms – ?
tphys – ?
tphysf – ?
dtp – ?
z – ?
zpars – ?
tb – ?
ecc – ?
vs – ?

void bse_instar(void): set collision matrix

void bse_star(int *kw, double *mass, double *mt, double *tm, double *tn, double *tscls, double *lums, double *GB, double *zpars)

star routine; shouldn’t need to be used often outside of BSE

Parameters

kw – stellar type
mass – initial mass
mt – current mass
tm – main sequence timescale
tn – nuclear timescale
tscls – timescales array (20)
lums – luminosities array (10)
GB – giant branch parameters (10)
zpars – evolution parameters based on metallicity (20)

void bse_hrdiag(double *mass, double *aj, double *mt, double *tm, double *tn, double *tscls, double *lums, double *GB, double *zpars, double *r, double *lum, int *kw, double *mc, double *rc, double *menv, double *renv, double *k2, double *bhspin)

void bse_kick(int *kw, double *m1, double *m1n, double *m2, double *ecc, double *sep, double *jorb, double *vk, int *snstar, double *r2, double *fallback, double *vs)

kick routine; shouldn’t need to be used often outside of BSE

Parameters

kw – stellar type
m1 – mass of star 1
m1n – new mass of star 1
m2 – mass of star 2
ecc – orbit eccentrity
sep – orbit semimajor axis
jorb – orbital angular momentum
vk – kick magnitude, can be used to help set initial pulsar properties if desired
snstar – which star (primary 1 or secondary 2) that goes SN
r2 – radius of companion, stars may collide if lucky kick
fallback – fraction of pre-SN mass that may fall back onto remnant. Can be 0. Kick strength is limited by such fall back of material
vs – old -> vs = velocity (3) of center of mass if bound, relative velocity at infinity if unbound new -> vs = three possible sets of velocities (3). Is an array of 12, correctly accounts for

void bse_mix(double *mass, double *mt, double *aj, int *kw, double *zpars, double *bhspin)

Mix routine, call from cmc_sscollision.c within the merge_two_stars routine. Helps merge troublesome systems that also wont pass through a common envelope.

Parameters

mass – ?
mt – ?
aj – ?
kw – ?
zpars – ?
ecsnp – ?

void bse_comenv(bse_binary *tempbinary, double *zpars, double *vs, int *fb)

?

Parameters

tempbinary – ?
zpars – ?
vs – ?
fb – ?
ecsnp – ?
ecsn_mlow – ?

void bse_set_using_cmc(void)

void bse_set_idum(int idum)

void bse_set_eddlimflag(int eddlimflag)

void bse_set_sigmadiv(double sigmadiv)

void bse_set_grflag(int grflag)

void bse_set_kickflag(int kickflag)

void bse_set_zsun(double zsun)

void bse_set_rembar_massloss(double rembar_massloss)

void bse_set_natal_kick_array(double natal_kick_array[10], long len_kick)

void bse_set_fprimc_array(double fprimc_array[16], long len_fprimc)

void bse_set_qcrit_array(double qcrit_array[16], long len_qcrit)

void bse_set_aic(int aic)

void bse_set_bdecayfac(int bdecayfac)

void bse_set_st_cr(int st_cr)

void bse_set_st_tide(int st_tide)

void bse_set_htpmb(int htpmb)

void bse_set_rejuvflag(int rejuvflag)

void bse_set_qcflag(int qcflag)

void bse_set_don_lim(double don_lim)

void bse_set_acc_lim(double acc_lim)

void bse_set_ussn(int ussn)

void bse_set_neta(double neta)

void bse_set_bwind(double bwind)

void bse_set_hewind(double hewind)

void bse_set_windflag(int windflag)

void bse_set_pisn(double pisn)

void bse_set_ecsn(double ecsn)

void bse_set_ecsn_mlow(double ecsn_mlow)

void bse_set_sigma(double sigma)

void bse_set_bhsigmafrac(double bhsigmafrac)

void bse_set_polar_kick_angle(int polar_kick_angle)

void bse_set_ifflag(int ifflag)

void bse_set_wdflag(int wdflag)

void bse_set_bhflag(int bhflag)

void bse_set_bhspinflag(int bhspinflag)

void bse_set_bhms_coll_flag(int bhms_coll_flag)

void bse_set_bhspinmag(double bhspinmag)

void bse_set_remnantflag(int remnantflag)

void bse_set_mxns(double mxns)

void bse_set_bconst(double bconst)

void bse_set_CK(double CK)

void bse_set_rejuv_fac(double rejuv_fac)

void bse_set_pts1(double pts1)

void bse_set_pts2(double pts2)

void bse_set_pts3(double pts3)

void bse_set_alpha1(double alpha1)

void bse_set_lambdaf(double lambdaf)

void bse_set_ceflag(int ceflag)

void bse_set_cemergeflag(int cemergeflag)

void bse_set_cekickflag(int cekickflag)

void bse_set_cehestarflag(int cehestarflag)

void bse_set_tflag(int tflag)

void bse_set_beta(double beta)

void bse_set_xi(double xi)

void bse_set_acc2(double acc2)

void bse_set_epsnov(double epsnov)

void bse_set_eddfac(double eddfac)

void bse_set_gamma(double gamma)

void bse_set_merger(double merger)

void bse_set_id1_pass(long int id1_pass)

void bse_set_id2_pass(long int id2_pass)

void bse_set_taus113state(struct rng_t113_state state, int first)

copies the C tausworthe rng state variables to the Fortran states

Parameters

state – C rng state
first – ?

double bse_get_alpha1(void)

double bse_get_lambdaf(void)

double bse_get_spp(int i, int j)

double bse_get_scm(int i, int j)

double bse_get_bpp(int i, int j)

double bse_get_bcm(int i, int j)

struct rng_t113_state bse_get_taus113state(void)

copies back the Fortran tausworthe rng state variables to the C state

Returns: state with variables copied from the Fortran rng

int icase_get(int i, int j)

?

Parameters

i – ?
j – ?

Returns

?

char *bse_get_sselabel(int kw)

?

Parameters: kw – ?
Returns: ?

char *bse_get_bselabel(int kw)

?

Parameters: kw – ?
Returns: ?

double bse_kick_speed(int *startype)

kick speed from distribution, taken directly from BSE code.

bse_evolv1() and bse_evolv2() return the kick speed, but the kick function is included here just in case you want to use it. may change startype in certain cases

Parameters: startype – star type
Returns: ?

double bse_rl(double q)

Roche lobe formula, taken directly from BSE code.

Parameters: q – ?
Returns: Returns R_L1/a, where q=m_1/m_2.