From 9fe43751bda44a6457834f20d9aff865bd7b4087 Mon Sep 17 00:00:00 2001 From: Christopher O'Connor Date: Fri, 27 Mar 2026 13:01:37 -0500 Subject: [PATCH 1/2] generalized/patched 3bbf + various bug fixes --- include/cmc/cmc.h | 20 +++-- include/cmc/cmc_vars.h | 2 + src/cmc/cmc_dynamics.c | 26 ++++--- src/cmc/cmc_dynamics_helper.c | 136 +++++++++++++++++++++------------- src/cmc/cmc_io.c | 9 ++- 5 files changed, 125 insertions(+), 68 deletions(-) diff --git a/include/cmc/cmc.h b/include/cmc/cmc.h index e06a007..b0a001d 100644 --- a/include/cmc/cmc.h +++ b/include/cmc/cmc.h @@ -680,11 +680,16 @@ typedef struct{ #define PARAMDOC_STREAMS "to run the serial version with the given number of random streams - primarily used to mimic the parallel version running with the same no.of processors" int STREAMS; /* Meagan - 3bb */ -#define PARAMDOC_THREEBODYBINARIES "toggles three-body binary formation (0=off, 1=on)" +#define PARAMDOC_THREEBODYBINARIES "toggles three-body binary formation (0=off, 1=on w/ pairing by mass, 2=on w/ pairing at random)" /** -* @brief toggles three-body binary formation (0=off, 1=on) +* @brief toggles three-body binary formation (0=off, 1=on w/ pairing by mass, 2=on w/ pairing at random) */ int THREEBODYBINARIES; +#define PARAMDOC_BINARY_HARDNESS_POWER "power law index for hardness of newly formed three-body binaries" +/** + * @brief Power-law index for hardness of three-body binaries (default is -3.0). + */ + double BINARY_HARDNESS_POWER; #define PARAMDOC_MIN_BINARY_HARDNESS "minimum hardness for newly formed three-body binaries" /** * @brief minimum hardness for newly formed three-body binaries @@ -1988,7 +1993,7 @@ void binint_do(long k, long kp, double rperi, double w[4], double W, double rcm, double simul_relax(gsl_rng *rng); double simul_relax_new(void); -int destroy_bbh(double m1, double m2,double a,double e,double nlocal,double sigma,struct rng_t113_state* rng_st); +int destroy_bbh(double m1,double m2,double a,double e,double nlocal,double sigma,struct rng_t113_state* rng_st); double simul_relax_new(void); void calc_sigma_r(long p, long N_LIMIT, double *sig_r, double *sig_sigma, long* sig_n, int r_0_mave_1); void break_wide_binaries(struct rng_t113_state* rng_st); @@ -1998,12 +2003,17 @@ double sigma_r(double r); // Meagan /* three-body binary formation */ void sort_three_masses(long sq, long *k1, long *k2, long *k3); -double get_eta(double eta_min, long k1, long k2, long k3, double vrel12[4], double vrel3[4]); +//double get_eta(double eta_min, long k1, long k2, long k3, double vrel12[4], double vrel3[4]); 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); -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_running, gsl_rng *rng); +//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_running, gsl_rng *rng); void calc_sigma_local(long k1, long p, long N_LIMIT, double *ave_local_mass, double *sigma_local); int remove_old_star(double time, long k); +// Chris O'C -- generalized 3BBF +void make_threebodybinary(double P_3bb, long k1, long k2, long k3, long form_binary, double eta_min, double eta_power, 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_running, gsl_rng *rng); +void shuffle_three_masses(long sq, long *k1, long *k2, long *k3, gsl_rng *rng); +double get_eta(double eta_min, double eta_power, long k1, long k2, long k3, double vrel12[4], double vrel3[4]); + // Meagan /* extra output for bhs */ void bh_count(long k); diff --git a/include/cmc/cmc_vars.h b/include/cmc/cmc_vars.h index 123a6e2..042337d 100644 --- a/include/cmc/cmc_vars.h +++ b/include/cmc/cmc_vars.h @@ -75,6 +75,8 @@ _EXTERN_ int BININITKT, STOPATCORECOLLAPSE; /* Meagan - 3bb */ _EXTERN_ int THREEBODYBINARIES, ONLY_FORM_BH_THREEBODYBINARIES;; _EXTERN_ double MIN_BINARY_HARDNESS; +// Chris: generalized 3bb +_EXTERN_ double BINARY_HARDNESS_POWER; _EXTERN_ double BINARY_DISTANCE_BREAKING; _EXTERN_ int BINARY_BREAKING_MIN; _EXTERN_ int USE_TT_FILE; diff --git a/src/cmc/cmc_dynamics.c b/src/cmc/cmc_dynamics.c index fadf6f2..aeb3b0d 100644 --- a/src/cmc/cmc_dynamics.c +++ b/src/cmc/cmc_dynamics.c @@ -30,7 +30,7 @@ void dynamics_apply(double dt, gsl_rng *rng) long sq, k1, k2, k3, form_binary; double n_threshold, triplet_count, num_triplets_averaged=200; double ave_local_mass, sigma_local, vrel_ave, v1[4], v2[4], v3[4], vrel12[4], vrel3[4]; - double eta_min=MIN_BINARY_HARDNESS, Y1, rate_3bb, rate_ave=0.0, P_3bb, P_ave=0.0; + double eta_min=MIN_BINARY_HARDNESS, eta_power=BINARY_HARDNESS_POWER, Y1, rate_3bb, rate_ave=0.0, P_3bb, P_ave=0.0; double clight10o7; double collisions_multiple; @@ -66,8 +66,6 @@ void dynamics_apply(double dt, gsl_rng *rng) pararootfprintf(logfile, "%s(): performing interactions:", __FUNCTION__); - - /* Added by Meagan 3/10/11 */ //================================ // Three-body binary formation | @@ -96,8 +94,13 @@ void dynamics_apply(double dt, gsl_rng *rng) { dt = SaveDt; form_binary = 0; // reset this to zero; later we decide whether to form a binary, and if so, set form_binary=1 - // Sort stars by mass (k1 is most massive) - sort_three_masses(sq, &k1, &k2, &k3); + // Sort stars by mass (k1 is most massive); >Feb 2026 only if using Morscher prescription for 3BBF + if (THREEBODYBINARIES == 1) { + sort_three_masses(sq, &k1, &k2, &k3); + } else { + shuffle_three_masses(sq, &k1, &k2, &k3, rng); + } + //sort_three_masses(sq, &k1, &k2, &k3); n_local = calc_n_local(get_global_idx(k1), BH_AVEKERNEL, N_LIMIT); // If density above threshold, check for 3bb formation if (n_local > n_threshold) { @@ -122,17 +125,18 @@ void dynamics_apply(double dt, gsl_rng *rng) // Calculate RATE of binary formation // Below is rate_3bb with all the velocity terms vrel_3 and vrel_12 replaced with the averaged local relative velocity, vrel_ave. We did this because we were finding that when we used the actual relative velocities, vrel12, if too large, the 3bb rate would be extremely low and binaries would not form (since it depends strongly on v: v^-9). When we replaced vrel12 with the average relative velocity (over 20 stars), the 3bb formation rate was high enough that binaries would form. We decided to use the average relative velocity for all relative velocity terms, vrel_3 and vrel_12. - - // *Note* Factor of 0.5 in front of rate_3bb ensures that our sampling method produces the correct overall analytic 3bb rate - rate_3bb = 0.5 * sqrt(2) * sqr(PI) * sqr(n_local) * pow(vrel_ave, -9) * pow(((star_m[get_global_idx(k1)] + star_m[get_global_idx(k2)]) * madhoc), 5.0) * pow(eta_min, -5.5) * (1.0 + 2.0*eta_min) * (1.0 + 2.0 * ((star_m[get_global_idx(k1)] + star_m[get_global_idx(k2)] + star_m[get_global_idx(k3)]) / (star_m[get_global_idx(k1)] + star_m[get_global_idx(k2)])) * eta_min); + // old rate from Morscher et al. 2013 + //rate_3bb = 0.5 * sqrt(2) * sqr(PI) * sqr(n_local) * pow(vrel_ave, -9) * pow(((star_m[get_global_idx(k1)] + star_m[get_global_idx(k2)]) * madhoc), 5.0) * pow(eta_min, -eta_power) * (1.0 + 2.0*eta_min) * (1.0 + 2.0 * ((star_m[get_global_idx(k1)] + star_m[get_global_idx(k2)] + star_m[get_global_idx(k3)]) / (star_m[get_global_idx(k1)] + star_m[get_global_idx(k2)])) * eta_min); + // updated rate from Atallah et al. 2024, plus user-specified power law + rate_3bb = 0.1328 * sqrt(10.) * sqr(n_local) * pow(vrel_ave,-9) * pow(((star_m[get_global_idx(k1)] + star_m[get_global_idx(k2)]) * madhoc), 5.0) * pow(eta_min, eta_power) * (1.0 + 2.0*eta_min) * (1.0 + 2.0 * ((star_m[get_global_idx(k1)] + star_m[get_global_idx(k2)] + star_m[get_global_idx(k3)]) / (star_m[get_global_idx(k1)] + star_m[get_global_idx(k2)])) * eta_min); // Calculate PROBABILITY of binary formation P_3bb = rate_3bb * (dt * ((double) clus.N_STAR)/log(GAMMA*((double) clus.N_STAR))); // print info on probability calculation, for innermost 200 triplets, in each timestep if (triplet_count <= num_triplets_averaged) { - P_ave += P_3bb; - rate_ave += rate_3bb; + P_ave += P_3bb; + rate_ave += rate_3bb; } /* //MPI: Commenting out since this will cause a deadlock, and is also this needs a reduction for each star. // For each timestep, output average rate and probability of 3bb formation, @@ -171,7 +175,7 @@ void dynamics_apply(double dt, gsl_rng *rng) /* Here is where new binary properties are calculated. if form_binary=1, the new binary will be created, but if form_binary=0, the interaction will be logged in lightcollision log */ - make_threebodybinary(P_3bb, k1, k2, k3, form_binary, eta_min, ave_local_mass, n_local, sigma_local, v1, v2, v3, vrel12, vrel3, delta_E_3bb, rng); + make_threebodybinary(P_3bb, k1, k2, k3, form_binary, eta_min, eta_power, ave_local_mass, n_local, sigma_local, v1, v2, v3, vrel12, vrel3, delta_E_3bb, rng); } else { // Probability of 3bb formation too low ==> No binary formed form_binary = 0; /* do nothing; check 3bb formation for next 3 stars */ diff --git a/src/cmc/cmc_dynamics_helper.c b/src/cmc/cmc_dynamics_helper.c index 8d16bac..451f348 100644 --- a/src/cmc/cmc_dynamics_helper.c +++ b/src/cmc/cmc_dynamics_helper.c @@ -257,8 +257,6 @@ long create_binary(int idx, int dyn_0_se_1) return(j); } - - /** * @brief Sort by mass: k1, k2, k3 in order of most massive to least massive. (parallel version of sort_three_masses) * @@ -312,6 +310,52 @@ void sort_three_masses(long sq, long *k1, long *k2, long *k3) { //parafprintf(threebbfile, "end of sort_three_masses function" ); } +/** +* @brief randomly permute three objects -- alternative to sort_three_masses +* +* @param sq index of the first of the three consecutive stars +* @param k1 variable to store the first star's index +* @param k2 variable to store the second star's index +* @param k3 variable to store the third star's index +*/ +void shuffle_three_masses(long sq, long *k1, long *k2, long *k3, gsl_rng *rng) { +/** Returns a random permutation of k1, k2, k3 */ + //long g_sq0 = get_global_idx(sq); + //long g_sq1 = get_global_idx(sq+1); + //long g_sq2 = get_global_idx(sq+2); + double y = rng_t113_dbl_new(curr_st); + if (y < 1.0/6.0) { + *k1 = sq; + *k2 = sq + 1; + *k3 = sq + 2; + } + else if (y < 1.0/3.0) { + *k1 = sq; + *k2 = sq + 2; + *k3 = sq + 1; + } + else if (y < 0.5) { + *k1 = sq + 1; + *k2 = sq; + *k3 = sq + 2; + } + else if (y < 2.0/3.0) { + *k1 = sq + 1; + *k2 = sq + 2; + *k3 = sq; + } + else if (y < 5.0/6.0) { + *k1 = sq + 2; + *k2 = sq; + *k3 = sq + 1; + } + else { + *k1 = sq + 2; + *k2 = sq + 1; + *k3 = sq; + } +} + /** * @brief Function for calculating quantities relevant for three-body binary formation @@ -381,11 +425,11 @@ void calc_3bb_encounter_dyns(long k1, long k2, long k3, double v1[4], double v2[ } - /** * @brief ? * * @param eta_min ? +* @param eta_power ? * @param k1 index of star 1 * @param k2 index of star 2 * @param k3 index of star 3 @@ -394,7 +438,7 @@ void calc_3bb_encounter_dyns(long k1, long k2, long k3, double v1[4], double v2[ * * @return ? */ -double get_eta(double eta_min, long k1, long k2, long k3, double vrel12[4], double vrel3[4]) +double get_eta(double eta_min, double eta_power, long k1, long k2, long k3, double vrel12[4], double vrel3[4]) { double kk, eta, norm, area_max, area, eta_test, comp_ymax, comp_y, true_y; double eta_max=50.0; @@ -414,31 +458,28 @@ double get_eta(double eta_min, long k1, long k2, long k3, double vrel12[4], doub comp_y=1e6; found_eta = 0; - /* normalize distribution, d(Rate)/d(eta) (differential rate of 3bb formation); take abs. value*/ - norm = fabs(1.0/(-1*(pow(eta_max, -5.5)) + (pow(eta_min, -5.5)) + (kk + 2) * ( -1*(pow(eta_max, -4.5)) + (pow(eta_min, -4.5)) ) + 2 * kk * ( -1*(pow(eta_max, -3.5)) + (pow(eta_min, -3.5))))); + if (eta_power==-2.0){ + eta_power=-2.01; // avoid divergence of distribution at eta^-2 + //TODO implement logarithmic distribution for this case + } + + /* normalize distribution, d(Rate)/d(eta) ~ eta^-power (differential rate of 3bb formation); take abs. value*/ + norm = fabs(1.0/(-1*(pow(eta_max, eta_power)) + (pow(eta_min, eta_power)) + (kk + 2) * ( -1*(pow(eta_max, eta_power+1.)) + (pow(eta_min, eta_power+1.)) ) + 2 * kk * ( -1*(pow(eta_max, eta_power+2.)) + (pow(eta_min, eta_power+2.))))); /* now start a while loop that will keep restarting as long as the last computed comp_value was rejected */ while (found_eta == 0) { - /* total area under comparison curve in range of eta from eta_min to eta_max */ - // area_max = fabs(3.25*norm*((pow(eta_max,-4.0)) - (pow(eta_min, -4.0)))); - area_max = fabs((2.0/7.0) * norm * kk * ((pow(eta_max,-3.5)) - (pow(eta_min, -3.5)))); + area_max = fabs(norm * kk * ((pow(eta_max,2.+eta_power)) - (pow(eta_min, 2.+eta_power))) / (2.+eta_power)); /* choose area in range (0, area_max); this area corresponds to a particular value of eta */ area = rng_t113_dbl_new(curr_st)*area_max; /* find the eta that corresponds to chosen area */ - /* test coordinate 1 (x) */ - // eta_test = pow(((-2.0*area/(13.0*norm)) + pow(eta_min, -4.0)), -0.25); - - // TODO: figure out this negative sign in front of 2 below - where did it come from?? - eta_test = pow(((-7.0*area/(2.0*norm*kk)) + pow(eta_min, -3.5)), (-2.0/7.0)); + eta_test = pow((pow(eta_min,2.+eta_power)+(pow(eta_max,2.+eta_power) - pow(eta_min,2.+eta_power))*(area/area_max)),1./(2.+eta_power)); //pow(((-3.0*area/(norm*kk)) + pow(eta_min, -3)), (-1.0/3.0)); /* now choose second test coordinate (y) by choosing a number between zero and the value of the comparison function at eta_test, comp_ymax */ - comp_ymax = kk * norm * pow(eta_test, -4.5); + comp_ymax = kk * norm * pow(eta_test, 1.+eta_power); comp_y = rng_t113_dbl_new(curr_st)*comp_ymax; - /* now accept or reject the point if it lies below/above true distribution: - d(rate)/d(eta) ~ norm*(5*eta^-6 + 8*eta^-5) */ - // true_y = norm*(5.0*pow(eta_test, -6.0) + 8.0*pow(eta_test, -5.0)); - true_y = norm*(5.5*pow(eta_test, -6.5) + (4.5 * kk +9.0) * (pow(eta_test, -5.5)) +7.0 * kk * (pow(eta_test, -4.5))); - if (comp_y < true_y) { /* fall within true distribution */ + /* now accept or reject the point if it lies below/above true distribution: */ + true_y = norm*(pow(eta_min,eta_power) - pow(eta_test,1.+eta_power)) + (2.+kk)*(pow(eta_min,1.+eta_power) - pow(eta_test,1.+eta_power)) + 2.*kk*(pow(eta_min,2.+eta_power) - pow(eta_test,2.+eta_power)); + if (comp_y < true_y) { /* fall within true distribution */ found_eta = 1; eta = eta_test; } @@ -467,16 +508,19 @@ double get_eta(double eta_min, long k1, long k2, long k3, double vrel12[4], doub * @param delta_E_3bb ? * @param rng gsl rng */ -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) +void make_threebodybinary(double P_3bb, long k1, long k2, long k3, long form_binary, double eta_min, double eta_power, 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) { - double PE_i, PE_f, KE_i, KE_f, KE_cmf_i, KE_cmf_f, delta_PE, delta_KE, delta_E, binary_cm; + double PE_i, PE_f, KE_i, KE_f, KE_cmf_i, KE_cmf_f, delta_PE, delta_KE, delta_E, binary_cm; double m1, m2, m3; // double ms, mb; double r1, r2, r3; + //double starrad1, starrad2; double cm_vel[4], v1_cmf[4], v2_cmf[4], v3_cmf[4]; - double eta, Eb, ecc_max, ecc, r_p, semi_major; + double eta, Eb, ecc_max, ecc, r_p, semi_major; //, roche1, roche2, qb; double vs_cmf[4], vb_cmf[4], vs[4], vb[4], angle3, angle4; long j, knew; + double Yp; //, mba, mbb, ms, rba, rbb, rs; + // long kba, kbb, ks; // Form new binary, set new binary/stellar properties, destroy old stars // Calculate Total potential energy of three stars (enforce cons. of potential energy at end to get new positions) // Find energy of new binary @@ -484,7 +528,6 @@ void make_threebodybinary(double P_3bb, long k1, long k2, long k3, long form_bin // Choose random direction for kick to single in COM frame - given direction of vs, know that vb points oppositely // Solve cons. of momentum and cons. of energy to obtain scalar velocities of binary and single, vb and vs - // INITIAL VALUES m1=star_m[get_global_idx(k1)]; m2=star_m[get_global_idx(k2)]; @@ -492,12 +535,11 @@ void make_threebodybinary(double P_3bb, long k1, long k2, long k3, long form_bin r1=star_r[get_global_idx(k1)]; r2=star_r[get_global_idx(k2)]; r3=star_r[get_global_idx(k3)]; + //starrad1=star[k1].rad; + //starrad2=star[k2].rad; + //rad3=star[k3].rad; - // Initial energy - PE_i = madhoc*(m1*star_phi[get_global_idx(k1)] + m2*star_phi[get_global_idx(k2)] + m3*star_phi[get_global_idx(k3)]); - KE_i = 0.5*madhoc*(m1 * sqr(v1[0]) + m2 * sqr(v2[0]) + m3 * sqr(v3[0])); - - // COM of candidate binary pair + // COM of candidate binary pair binary_cm = (m1 * r1 + m2 * r2)/(m1 + m2); for (j=1; j<=3; j++) { @@ -522,37 +564,31 @@ void make_threebodybinary(double P_3bb, long k1, long k2, long k3, long form_bin // Initial kinetic energy of three stars in COM frame (relative to COM motion) KE_cmf_i = 0.5*madhoc*(m1 * sqr(v1_cmf[0]) + m2 * sqr(v2_cmf[0]) + m3 * sqr(v3_cmf[0])); - - // COMPUTE NEW QUANTITIES - // Binary orbital properties: - // binding energy: choose eta, then Eb = eta * / sigma^2 - // eccentricity: chosen from thermal distribution - eta = get_eta(eta_min, k1, k2, k3, vrel12, vrel3); - Eb = -0.5*eta * ave_local_mass * sqr(sigma_local); //Note: ave_local_mass is already multiplied by madhoc + eta = get_eta(eta_min, eta_power, k1, k2, k3, vrel12, vrel3); + Eb = -0.5*eta * ave_local_mass * sqr(sigma_local); //Note: ave_local_mass is already multiplied by madhoc // **Note binding energy is NEGATIVE - ecc_max = 1.0; + //ecc_max = 1.0; + semi_major = m1 * m2 * sqr(madhoc) / (-2*Eb); + // calculate Roche lobe radii (Eggleton 1983) + //roche1 = ((0.6 * pow(m2/m1, -2./3.) + log(1.+pow(m2/m1, -1./3.))) / (0.49 * pow(m2/m1, -2./3.))) * star[k1].rad; + //roche2 = ((0.6 * pow(m2/m1, 2./3.) + log(1.+pow(m2/m1, 1./3.))) / (0.49 * pow(m2/m1, 2./3.))) * star[k2].rad; + //ecc_max -= 2.*(roche1+roche2)/semi_major; // prevent interacting binaries from 3BBF; + ecc_max = 0.7; //1. - 2.*(((0.6 * pow(m2/m1, -2./3.) + log(1.+pow(m2/m1, -1./3.))) / (0.49 * pow(m2/m1, -2./3.))) * star[k1].rad + ((0.6 * pow(m2/m1, 2./3.) + log(1.+pow(m2/m1, 1./3.))) / (0.49 * pow(m2/m1, 2./3.))) * star[k2].rad)/semi_major; // prevent interacting binaries from 3BBF; ecc = ecc_max * sqrt(rng_t113_dbl_new(curr_st)); r_p = m1 * m2 * madhoc / (eta * sqr(sigma_local)); //note, for units, have madhoc^2 in numerator, and madhoc in the denominator ====> single power madhoc in numerator - semi_major = m1 * m2 * sqr(madhoc) / (-2*Eb); // Using cons. of momentum and energy, can calculate scalar velocities of binary and single vs_cmf[0] = sqrt((2.0*(KE_cmf_i-Eb)*((m1 + m2)/m3))/((madhoc)*(m1 + m2 + m3))); vb_cmf[0] = 1.0*vs_cmf[0]*m3/(m1 + m2); - // Choose random direction for motion of single star - angle3 = rng_t113_dbl_new(curr_st) * PI; // polar angle - [0, PI) + // Choose random direction for motion of single star + angle3 = acos(2.0*rng_t113_dbl_new(curr_st) - 1.); // polar angle - [0, PI) angle4 = rng_t113_dbl_new(curr_st) * 2.0 * PI; // azimuthal angle - [0, 2*PI) // Compute vector velocities of single and binary vs_cmf[1] = vs_cmf[0] * sin(angle3) * cos(angle4); vs_cmf[2] = vs_cmf[0] * sin(angle3) * sin(angle4); vs_cmf[3] = vs_cmf[0] * cos(angle3); - /* Fix direction of velocity of binary to be opposite to that of single - add Pi to each of the angles - (polar angle and azimuthal) that were used to orient the velocity of the single */ - vb_cmf[1] = vb_cmf[0] * sin(PI - angle3) * cos(angle4 + PI); - vb_cmf[2] = vb_cmf[0] * sin(PI - angle3) * sin(angle4 + PI); - vb_cmf[3] = vb_cmf[0] * cos(PI - angle3); - - vs_cmf[0] = sqrt(sqr(vs_cmf[1]) + sqr(vs_cmf[2]) + sqr(vs_cmf[3])); + vs_cmf[0] = sqrt(sqr(vs_cmf[1]) + sqr(vs_cmf[2]) + sqr(vs_cmf[3])); vb_cmf[0] = sqrt(sqr(vb_cmf[1]) + sqr(vb_cmf[2]) + sqr(vb_cmf[3])); /* Final COM energy - Binding energy of new binary plus kinetic energy of single and binary, @@ -678,10 +714,8 @@ void make_threebodybinary(double P_3bb, long k1, long k2, long k3, long form_bin } - } - /** * @brief ? * @@ -2795,8 +2829,8 @@ void break_wide_binaries(struct rng_t113_state* rng_st) star[i].vr *= exc_ratio; star[i].vt *= exc_ratio; Eexcess -= E_dump; - set_star_EJ(k); - star[k].interacted = 1; + set_star_EJ(i); // formerly set_star_EJ(k); Chris O'C fixed this bug in August 2024 + star[i].interacted = 1; // formerly star[k].interacted = 1; see above } } } diff --git a/src/cmc/cmc_io.c b/src/cmc/cmc_io.c index d25ecdf..74207ab 100644 --- a/src/cmc/cmc_io.c +++ b/src/cmc/cmc_io.c @@ -989,10 +989,15 @@ if(myid==0) { sscanf(values, "%d", &procs); parsed.STREAMS = 1; /* Meagan: new parameters for three-body binaries: THREEBODYBINARIES, MIN_BINARY_HARDNESS, ONLY_FORM_BH_THREEBODYBINARIES */ + // Chris O'C: added BINARY_HARDNESS_POWER March 2026 } else if (strcmp(parameter_name, "THREEBODYBINARIES") == 0) { PRINT_PARSED(PARAMDOC_THREEBODYBINARIES); sscanf(values, "%d", &THREEBODYBINARIES); parsed.THREEBODYBINARIES = 1; + } else if (strcmp(parameter_name, "BINARY_HARDNESS_POWER") == 0) { + PRINT_PARSED(PARAMDOC_BINARY_HARDNESS_POWER); + sscanf(values, "%lf", &BINARY_HARDNESS_POWER); + parsed.BINARY_HARDNESS_POWER = 1; } else if (strcmp(parameter_name, "MIN_BINARY_HARDNESS") == 0) { PRINT_PARSED(PARAMDOC_MIN_BINARY_HARDNESS); sscanf(values, "%lf", &MIN_BINARY_HARDNESS); @@ -1700,6 +1705,8 @@ if(myid==0) { /*Meagan: new parameters for 3-body binary formation*/ CHECK_PARSED(THREEBODYBINARIES, 0, PARAMDOC_THREEBODYBINARIES); CHECK_PARSED(MIN_BINARY_HARDNESS, 5.0, PARAMDOC_MIN_BINARY_HARDNESS); + // Chris O'C - new parameter to adjust the hardness distribution for 3BBF + CHECK_PARSED(BINARY_HARDNESS_POWER, -5.5, PARAMDOC_BINARY_HARDNESS_POWER); CHECK_PARSED(ONLY_FORM_BH_THREEBODYBINARIES, 1, PARAMDOC_ONLY_FORM_BH_THREEBODYBINARIES); // default - 1: three-body binary formation only allowed for black holes CHECK_PARSED(BH_LOSS_CONE, 0, PARAMDOC_BH_LOSS_CONE); @@ -2496,7 +2503,7 @@ MPI: In the parallel version, IO is done in the following way. Some files requir //MPI: Headers are written out only by the root node. // print header if(RESTART_TCOUNT <= 0){ - pararootfprintf(escfile, "#1:tcount #2:t #3:m[MSUN] #4:r #5:vr #6:vt #7:r_peri #8:r_apo #9:Rtidal #10:phi_rtidal #11:phi_zero #12:E #13:J #14:id #15:binflag #16:m0[MSUN] #17:m1[MSUN] #18:id0 #19:id1 #20:a #21:e #22:startype #23:bin_startype0 #24:bin_startype1 #25:rad0 #26:rad1 #27:tb #28:lum0 #29:lum1 #30:massc0 #31:massc1 #32:radc0 #33:radc1 #34:menv0 #35:menv1 #36:renv0 #37:renv1 #38:tms0 #39:tms1 #40:dmdt0 #41:dmdt1 #42:radrol0 #43:radrol1 #44:ospin0 #45:ospin1 #46:B0 #47:B1 #48:formation0 #49:formation1 #50:bacc0 #51:bacc1 #52:tacc0 $53:tacc1 #54:mass0_0 #55:mass0_1 #56:epoch0 #57:epoch1 #58:bhspin #59:bhspin1 #60:bhspin2 #61:ospin #62:B #63:formation\n"); + pararootfprintf(escfile, "#1:tcount #2:t #3:m[MSUN] #4:r #5:vr #6:vt #7:r_peri #8:r_apo #9:Rtidal #10:phi_rtidal #11:phi_zero #12:E #13:J #14:id #15:binflag #16:m0[MSUN] #17:m1[MSUN] #18:id0 #19:id1 #20:a #21:e #22:startype #23:bin_startype0 #24:bin_startype1 #25:rad0 #26:rad1 #27:tb #28:lum0 #29:lum1 #30:massc0 #31:massc1 #32:radc0 #33:radc1 #34:menv0 #35:menv1 #36:renv0 #37:renv1 #38:tms0 #39:tms1 #40:dmdt0 #41:dmdt1 #42:radrol0 #43:radrol1 #44:ospin0 #45:ospin1 #46:B0 #47:B1 #48:formation0 #49:formation1 #50:bacc0 #51:bacc1 #52:tacc0 #53:tacc1 #54:mass0_0 #55:mass0_1 #56:epoch0 #57:epoch1 #58:bhspin #59:bhspin1 #60:bhspin2 #61:ospin #62:B #63:formation\n"); // there was a $ in there and it was goofing up automatic parsing of outputs -- fixed by Chris O'C, March 2026 // print header pararootfprintf(triplefile, "#1:time #2:min0 #3:min1 #4:mout #5:Rin0 #6:Rin1 #7:Rout #8:ain #9:aout #10:ein #11:eout #12:ktypein0 #13:ktypein1 #14:ktypeout #15:Tlk_quad #16:Tlk_oct#17:eps_oct #18:T_GR #19:eps_GR\n"); // print header From f1fdeb0f6816327fffffa69b5e6c75bbd8d9f0c0 Mon Sep 17 00:00:00 2001 From: Christopher O'Connor Date: Fri, 27 Mar 2026 14:51:23 -0500 Subject: [PATCH 2/2] restored mistakenly deleted lines in cmc_dynamics_helper.c --- src/cmc/cmc_dynamics_helper.c | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/cmc/cmc_dynamics_helper.c b/src/cmc/cmc_dynamics_helper.c index 451f348..fe92fe9 100644 --- a/src/cmc/cmc_dynamics_helper.c +++ b/src/cmc/cmc_dynamics_helper.c @@ -588,6 +588,12 @@ void make_threebodybinary(double P_3bb, long k1, long k2, long k3, long form_bin vs_cmf[2] = vs_cmf[0] * sin(angle3) * sin(angle4); vs_cmf[3] = vs_cmf[0] * cos(angle3); + /* Fix direction of velocity of binary to be opposite to that of single - add Pi to each of the angles + (polar angle and azimuthal) that were used to orient the velocity of the single */ + vb_cmf[1] = vb_cmf[0] * sin(PI - angle3) * cos(angle4 + PI); + vb_cmf[2] = vb_cmf[0] * sin(PI - angle3) * sin(angle4 + PI); + vb_cmf[3] = vb_cmf[0] * cos(PI - angle3); + vs_cmf[0] = sqrt(sqr(vs_cmf[1]) + sqr(vs_cmf[2]) + sqr(vs_cmf[3])); vb_cmf[0] = sqrt(sqr(vb_cmf[1]) + sqr(vb_cmf[2]) + sqr(vb_cmf[3]));