From b8bc7d78846372749e1dce2ca485a0158fe81fcf Mon Sep 17 00:00:00 2001 From: elenagonzalez870 Date: Thu, 26 Mar 2026 12:53:22 -0500 Subject: [PATCH] * Definining pred_id1 and pred_id2 in do_stellar_evolution to prevent id=0 in nsformation and bhformation files. The reason for this is that within handle_bse_outcome (called before writing out to these files), an object can be destroyed if a merger happens, and hence the IDs in the binary[kb] structure are zeroed. I would recommned to move the writing to these files inside handle_bse_outcome such that there's more consistency. * Passing VKO as a pointer into handle_bse_outcome such that it gets updated when the net kicks are calculated. This bug caused neutron stars (NSs) and black holes (BHs) born in binaries in the nsformation and bhformation files to appear as though they did not receive kicks. The compact objects were receiving kicks; they were simply not written to the files correctly. --- include/cmc/cmc.h | 2 +- src/cmc/cmc_stellar_evolution.c | 78 ++++++++++++++++++--------------- 2 files changed, 43 insertions(+), 37 deletions(-) diff --git a/include/cmc/cmc.h b/include/cmc/cmc.h index e06a007..f153b6b 100644 --- a/include/cmc/cmc.h +++ b/include/cmc/cmc.h @@ -1898,7 +1898,7 @@ void stellar_evolution_init(void); void restart_stellar_evolution(void); void do_stellar_evolution(gsl_rng *rng); void write_stellar_data(void); -void handle_bse_outcome(long k, long kb, double *vs, double tphysf, int kprev0, int kprev1); +void handle_bse_outcome(long k, long kb, double *vs, double tphysf, int kprev0, int kprev1, double *VKO); void cp_binmemb_to_star(long k, int kbi, long knew); void cp_SEvars_to_newstar(long oldk, int kbi, long knew); void cp_m_to_newstar(long oldk, int kbi, long knew); diff --git a/src/cmc/cmc_stellar_evolution.c b/src/cmc/cmc_stellar_evolution.c index 3e8c369..6469250 100644 --- a/src/cmc/cmc_stellar_evolution.c +++ b/src/cmc/cmc_stellar_evolution.c @@ -106,6 +106,7 @@ void stellar_evolution_init(void){ int kprev0=-100; int kprev1=-100; binary_t tempbinary; + double VKO; /* SSE */ /* bse_set_hewind(0.5); */ @@ -343,7 +344,7 @@ void stellar_evolution_init(void){ &(binary[kb].bse_tb), &(binary[kb].e), vs,&(binary[kb].bse_bhspin[0])); *curr_st=bse_get_taus113state(); - handle_bse_outcome(k, kb, vs, tphysf, kprev0, kprev1); + handle_bse_outcome(k, kb, vs, tphysf, kprev0, kprev1, &VKO); } else { eprintf("totally confused!\n"); exit_cleanly(-1, __FUNCTION__); @@ -568,12 +569,12 @@ void do_stellar_evolution(gsl_rng *rng) } } - if (WRITE_MOREPULSAR_INFO) { - if (kprev!=13 && star[k].se_k==13) { // newly formed NS - parafprintf(newnsfile, "%.18g %g 0 %ld %g %g %g %g %g %d", TotalTime, star_r[g_k], star[k].id,star[k].se_zams_mass,star[k].se_mass, star[k].se_mt, star[k].se_scm_formation, VKO, kprev); - parafprintf (newnsfile, "\n"); - } - } + if (WRITE_MOREPULSAR_INFO) { + if (kprev!=13 && star[k].se_k==13) { // newly formed NS + parafprintf(newnsfile, "%.18g %g 0 %ld %g %g %g %g %g %d", TotalTime, star_r[g_k], star[k].id,star[k].se_zams_mass,star[k].se_mass, star[k].se_mt, star[k].se_scm_formation, VKO, kprev); + parafprintf (newnsfile, "\n"); + } + } if (WRITE_BH_INFO) { if (kprev!=14 && star[k].se_k==14) { // newly formed BH @@ -657,20 +658,25 @@ void do_stellar_evolution(gsl_rng *rng) fprintf(stderr, "tphys=%g tphysf=%g kstar1=%d kstar2=%d m1=%g m2=%g r1=%g r2=%g l1=%g l2=%g \n", binary[kb].bse_tphys, tphysf, binary[kb].bse_kw[0], binary[kb].bse_kw[1], binary[kb].bse_mass[0], binary[kb].bse_mass[1], binary[kb].bse_radius[0], binary[kb].bse_radius[1], binary[kb].bse_lum[0], binary[kb].bse_lum[1]); fprintf(stderr, "k= %ld kb=%ld star_id=%ld bin_id1=%ld bin_id2=%ld \n", k, kb, star[k].id, binary[kb].id1, binary[kb].id2); exit(1); - } + } + + /*EGP: Adding these two new variables such that IDs can be used to write out to the files below even if the objects are deleted inside handle_bse_outcome. */ + /*EGP: Note that when we delete objects, we zero out some quantities like the IDs, but a lot of bse params are not, hence why only the IDs were affected. */ + long prev_id1 = binary[kb].id1; + long prev_id2 = binary[kb].id2; - handle_bse_outcome(k, kb, vs, tphysf, kprev0, kprev1); + handle_bse_outcome(k, kb, vs, tphysf, kprev0, kprev1, &VKO); if (WRITE_BH_INFO) { if (kprev0!=14 && binary[kb].bse_kw[0]==14) { // newly formed BH - parafprintf(newbhfile, "%.18g %g 1 %ld %g %g %g %g %g", TotalTime, star_r[g_k], binary[kb].id1, binary[kb].bse_zams_mass[0], binary[kb].bse_mass0[0], binary[kb].bse_mass[0], binary[kb].bse_bhspin[0], VKO); + parafprintf(newbhfile, "%.18g %g 1 %ld %g %g %g %g %g", TotalTime, star_r[g_k], prev_id1, binary[kb].bse_zams_mass[0], binary[kb].bse_mass0[0], binary[kb].bse_mass[0], binary[kb].bse_bhspin[0], VKO); for (ii=0; ii<16; ii++){ parafprintf (newbhfile, " %g", vs[ii]); } parafprintf (newbhfile, "\n"); } if (kprev1!=14 && binary[kb].bse_kw[1]==14 && binary[kb].id2 != 0) { // newly formed BH - parafprintf(newbhfile, "%.18g %g 1 %ld %g %g %g %g %g", TotalTime, star_r[g_k], binary[kb].id2, binary[kb].bse_zams_mass[1],binary[kb].bse_mass0[1], binary[kb].bse_mass[1], binary[kb].bse_bhspin[1],VKO); + parafprintf(newbhfile, "%.18g %g 1 %ld %g %g %g %g %g", TotalTime, star_r[g_k], prev_id2, binary[kb].bse_zams_mass[1],binary[kb].bse_mass0[1], binary[kb].bse_mass[1], binary[kb].bse_bhspin[1],VKO); for (ii=0; ii<16; ii++){ parafprintf (newbhfile, " %g", vs[ii]); } @@ -680,11 +686,11 @@ void do_stellar_evolution(gsl_rng *rng) if (WRITE_MOREPULSAR_INFO) { if (kprev0!=13 && binary[kb].bse_kw[0]==13) { // newly formed NS - parafprintf(newnsfile, "%.18g %g 1 %ld %g %g %g %g %g %d", TotalTime, star_r[g_k], binary[kb].id1, binary[kb].bse_zams_mass[0], binary[kb].bse_mass0[0], binary[kb].bse_mass[0], binary[kb].bse_bcm_formation[0], VKO, kprev0); + parafprintf(newnsfile, "%.18g %g 1 %ld %g %g %g %g %g %d", TotalTime, star_r[g_k], prev_id1, binary[kb].bse_zams_mass[0], binary[kb].bse_mass0[0], binary[kb].bse_mass[0], binary[kb].bse_bcm_formation[0], VKO, kprev0); parafprintf(newnsfile, "\n"); } if (kprev1!=13 && binary[kb].bse_kw[1]==13 && binary[kb].id2 != 0) { // newly formed NS - parafprintf(newnsfile, "%.18g %g 1 %ld %g %g %g %g %g %d", TotalTime, star_r[g_k], binary[kb].id2, binary[kb].bse_zams_mass[1],binary[kb].bse_mass0[1], binary[kb].bse_mass[1], binary[kb].bse_bcm_formation[1],VKO,kprev1); + parafprintf(newnsfile, "%.18g %g 1 %ld %g %g %g %g %g %d", TotalTime, star_r[g_k], prev_id2, binary[kb].bse_zams_mass[1],binary[kb].bse_mass0[1], binary[kb].bse_mass[1], binary[kb].bse_bcm_formation[1],VKO,kprev1); parafprintf(newnsfile, "\n"); } } @@ -804,14 +810,14 @@ void write_stellar_data(void){ * @param vs ? * @param tphysf ? */ -void handle_bse_outcome(long k, long kb, double *vs, double tphysf, int kprev0, int kprev1) +void handle_bse_outcome(long k, long kb, double *vs, double tphysf, int kprev0, int kprev1, double *VKO) { int j, jj; long knew=0, knewp=0, convert; - double dtp, VKO; + double dtp; knew = 0; - VKO = 0.0; + *VKO = 0.0; /* PK: extract some stellar/binary info from BSE's bcm array */ /* but don't loop forever, just until maximum array length. */ @@ -908,14 +914,14 @@ void handle_bse_outcome(long k, long kb, double *vs, double tphysf, int kprev0, star[k].vr += 0.0; star[k].vt += 0.0; set_star_EJ(k); - VKO = 0.0; + *VKO = 0.0; } else if (vs[4] <= 0.0) { /* If one kick */ star[k].vr += vs[3] * 1.0e5 / (units.l/units.t); vt_add_kick(&(star[k].vt),vs[1],vs[2], curr_st); //star[k].vt += sqrt(vs[1]*vs[1]+vs[2]*vs[2]) * 1.0e5 / (units.l/units.t); set_star_EJ(k); - VKO = sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]); + *VKO = sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]); } else { /* Two kicks */ star[k].vr += vs[3] * 1.0e5 / (units.l/units.t) + vs[7] * 1.0e5 / (units.l/units.t); @@ -923,7 +929,7 @@ void handle_bse_outcome(long k, long kb, double *vs, double tphysf, int kprev0, vt_add_kick(&(star[k].vt),vs[5],vs[6], curr_st);//will mean a different angle for each kick, should be ok as this is a randomised process... //star[k].vt += sqrt(vs[1]*vs[1]+vs[2]*vs[2]) * 1.0e5 / (units.l/units.t) + sqrt(vs[5]*vs[5]+vs[6]*vs[6]) * 1.0e5 / (units.l/units.t); set_star_EJ(k); - VKO = sqrt(vs[1]*vs[1] + vs[2]*vs[2] + vs[3]*vs[3]) + sqrt(vs[5]*vs[5]+vs[6]*vs[6]+vs[7]*vs[7]); + *VKO = sqrt(vs[1]*vs[1] + vs[2]*vs[2] + vs[3]*vs[3]) + sqrt(vs[5]*vs[5]+vs[6]*vs[6]+vs[7]*vs[7]); } if (sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]) != 0.0) { //dprintf("birth kick of %f km/s\n", sqrt(vs[0]*vs[0]+vs[1]*vs[1]+vs[2]*vs[2])); @@ -982,7 +988,7 @@ void handle_bse_outcome(long k, long kb, double *vs, double tphysf, int kprev0, vt_add_kick(&(star[knew].vt),vs[1],vs[2], curr_st); set_star_EJ(knew); - VKO = sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]); + *VKO = sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]); star[knewp].vr += vs[7] * 1.0e5 / (units.l/units.t); vt_add_kick(&(star[knewp].vt),vs[5],vs[6], curr_st); //star[knewp].vt += sqrt(vs[5]*vs[5]+vs[6]*vs[6]) * 1.0e5 / (units.l/units.t); @@ -993,7 +999,7 @@ void handle_bse_outcome(long k, long kb, double *vs, double tphysf, int kprev0, vt_add_kick(&(star[knewp].vt),vs[1],vs[2], curr_st); //star[knewp].vt += sqrt(vs[1]*vs[1]+vs[2]*vs[2]) * 1.0e5 / (units.l/units.t); set_star_EJ(knewp); - VKO = sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]); + *VKO = sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]); star[knew].vr += vs[7] * 1.0e5 / (units.l/units.t); vt_add_kick(&(star[knew].vt),vs[5],vs[6], curr_st); //star[knew].vt += sqrt(vs[5]*vs[5]+vs[6]*vs[6]) * 1.0e5 / (units.l/units.t); //minus at front? @@ -1009,7 +1015,7 @@ void handle_bse_outcome(long k, long kb, double *vs, double tphysf, int kprev0, vt_add_kick(&(star[knew].vt),vs[9],vs[10], curr_st); //star[knew].vt += sqrt(vs[1]*vs[1]+vs[2]*vs[2]) * 1.0e5 / (units.l/units.t) + sqrt(vs[9]*vs[9]+vs[10]*vs[10]) * 1.0e5 / (units.l/units.t); set_star_EJ(knew); - VKO = sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]+vs[9]*vs[9]+vs[10]*vs[10]+vs[11]*vs[11]); + *VKO = sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]+vs[9]*vs[9]+vs[10]*vs[10]+vs[11]*vs[11]); star[knewp].vr += vs[3] * 1.0e5 / (units.l/units.t) + vs[7] * 1.0e5 / (units.l/units.t); vt_add_kick(&(star[knewp].vt),vs[5],vs[6], curr_st);//finalise companion vt... @@ -1022,7 +1028,7 @@ void handle_bse_outcome(long k, long kb, double *vs, double tphysf, int kprev0, vt_add_kick(&(star[knewp].vt),vs[9],vs[10], curr_st); //star[knewp].vt += sqrt(vs[1]*vs[1]+vs[2]*vs[2]) * 1.0e5 / (units.l/units.t) + sqrt(vs[9]*vs[9]+vs[10]*vs[10]) * 1.0e5 / (units.l/units.t); set_star_EJ(knewp); - VKO = sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]+vs[9]*vs[9]+vs[10]*vs[10]+vs[11]*vs[11]); + *VKO = sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]+vs[9]*vs[9]+vs[10]*vs[10]+vs[11]*vs[11]); star[knew].vr += vs[3] * 1.0e5 / (units.l/units.t) + vs[7] * 1.0e5 / (units.l/units.t); vt_add_kick(&(star[knew].vt),vs[5],vs[6], curr_st); //star[knew].vt += sqrt(vs[1]*vs[1]+vs[2]*vs[2]) * 1.0e5 / (units.l/units.t) + sqrt(vs[5]*vs[5]+vs[6]*vs[6]) * 1.0e5 / (units.l/units.t); @@ -1036,7 +1042,7 @@ void handle_bse_outcome(long k, long kb, double *vs, double tphysf, int kprev0, //star[knew].vt += sqrt(vs[1]*vs[1]+vs[2]*vs[2]) * 1.0e5 / (units.l/units.t); vt_add_kick(&(star[knew].vt),vs[1],vs[2], curr_st); set_star_EJ(knew); - VKO = sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]); + *VKO = sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]); star[knewp].vr += vs[7] * 1.0e5 / (units.l/units.t) + vs[11] * 1.0e5 / (units.l/units.t); vt_add_kick(&(star[knewp].vt),vs[5],vs[6], curr_st); vt_add_kick(&(star[knewp].vt),vs[9],vs[10], curr_st); @@ -1047,7 +1053,7 @@ void handle_bse_outcome(long k, long kb, double *vs, double tphysf, int kprev0, vt_add_kick(&(star[knewp].vt),vs[1],vs[2], curr_st); //star[knewp].vt += sqrt(vs[1]*vs[1]+vs[2]*vs[2]) * 1.0e5 / (units.l/units.t); set_star_EJ(knewp); - VKO = sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]); + *VKO = sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]); star[knew].vr += vs[7] * 1.0e5 / (units.l/units.t) + vs[11] * 1.0e5 / (units.l/units.t); vt_add_kick(&(star[knew].vt),vs[5],vs[6], curr_st); vt_add_kick(&(star[knew].vt),vs[9],vs[10], curr_st); @@ -1059,7 +1065,7 @@ void handle_bse_outcome(long k, long kb, double *vs, double tphysf, int kprev0, star[knew].vr += 0.0; star[knew].vt += 0.0; set_star_EJ(knew); - VKO = 0.0; + *VKO = 0.0; star[knewp].vr += 0.0; star[knewp].vt += 0.0; set_star_EJ(knewp); @@ -1100,14 +1106,14 @@ void handle_bse_outcome(long k, long kb, double *vs, double tphysf, int kprev0, vt_add_kick(&(star[knew].vt),vs[1],vs[2], curr_st); //star[knew].vt += sqrt(vs[1]*vs[1]+vs[2]*vs[2]) * 1.0e5 / (units.l/units.t); set_star_EJ(knew); - VKO = sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]); + *VKO = sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]); } else if (vs[0]>0.0 && vs[4]>0.0 && vs[8]<=0.0 && (vs[0]==vs[4])) { /* 1 SN occured disrupted system then one star killed itself */ star[knew].vr += vs[3] * 1.0e5 / (units.l/units.t); //star[knew].vt += sqrt(vs[1]*vs[1]+vs[2]*vs[2]) * 1.0e5 / (units.l/units.t); vt_add_kick(&(star[knew].vt),vs[1],vs[2], curr_st); set_star_EJ(knew); - VKO = sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]); + *VKO = sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]); } else if (vs[0]>0.0 && vs[4]>0.0 && vs[8]<=0.0 && (vs[0]!=vs[4])) { /* 2 SNe then system mergers, star feels both kicks (i.e. it is COM) */ /* NOTE: merger remnants of double compact (NS or BH) binaries do not receive kicks in BSE as of yet! @@ -1117,13 +1123,13 @@ void handle_bse_outcome(long k, long kb, double *vs, double tphysf, int kprev0, vt_add_kick(&(star[knew].vt),vs[5],vs[6], curr_st); //star[knew].vt += sqrt(vs[1]*vs[1]+vs[2]*vs[2]) * 1.0e5 / (units.l/units.t) + sqrt(vs[5]*vs[5]+vs[6]*vs[6]) * 1.0e5 / (units.l/units.t); set_star_EJ(knew); - VKO = sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]+vs[5]*vs[5]+vs[6]*vs[6]+vs[7]*vs[7]); + *VKO = sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]+vs[5]*vs[5]+vs[6]*vs[6]+vs[7]*vs[7]); } else { /* No kick - going to set_star_EJ seems overkill (here and elsewhere). */ star[knew].vr += 0.0; star[knew].vt += 0.0; set_star_EJ(knew); - VKO = 0.0; + *VKO = 0.0; } if (sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]) != 0.0) { //dprintf("birth kick of %f km/s\n", sqrt(vs[0]*vs[0]+vs[1]*vs[1]+vs[2]*vs[2])); @@ -1193,14 +1199,14 @@ void handle_bse_outcome(long k, long kb, double *vs, double tphysf, int kprev0, vt_add_kick(&(star[knew].vt),vs[1],vs[2], curr_st); //star[knew].vt += sqrt(vs[1]*vs[1]+vs[2]*vs[2]) * 1.0e5 / (units.l/units.t); set_star_EJ(knew); - VKO = sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]); + *VKO = sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]); } else if (vs[0]>0.0 && vs[4]>0.0 && vs[8]<=0.0 && (vs[0]==vs[4])) { /* 1 SN occured disrupted system then one star killed itself */ star[knew].vr += vs[3] * 1.0e5 / (units.l/units.t); vt_add_kick(&(star[knew].vt),vs[1],vs[2], curr_st); //star[knew].vt += sqrt(vs[1]*vs[1]+vs[2]*vs[2]) * 1.0e5 / (units.l/units.t); set_star_EJ(knew); - VKO = sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]); + *VKO = sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]); } else if (vs[0]>0.0 && vs[4]>0.0 && vs[8]<=0.0 && (vs[0]!=vs[4])) { /* 2 SNe then system mergers, star feels both kicks (i.e. it is COM) */ /* NOTE: merger remnants of double compact (NS or BH) binaries do not receive kicks in BSE as of yet! @@ -1210,13 +1216,13 @@ void handle_bse_outcome(long k, long kb, double *vs, double tphysf, int kprev0, vt_add_kick(&(star[knew].vt),vs[5],vs[6], curr_st); //star[knew].vt += sqrt(vs[1]*vs[1]+vs[2]*vs[2]) * 1.0e5 / (units.l/units.t) + sqrt(vs[5]*vs[5]+vs[6]*vs[6]) * 1.0e5 / (units.l/units.t); set_star_EJ(knew); - VKO = sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]+vs[5]*vs[5]+vs[6]*vs[6]+vs[7]*vs[7]); + *VKO = sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]+vs[5]*vs[5]+vs[6]*vs[6]+vs[7]*vs[7]); } else { /* No kick - going to set_star_EJ seems overkill (here and elsewhere). */ star[knew].vr += 0.0; star[knew].vt += 0.0; set_star_EJ(knew); - VKO = 0.0; + *VKO = 0.0; } if (sqrt(vs[1]*vs[1]+vs[2]*vs[2]+vs[3]*vs[3]) != 0.0) { //dprintf("birth kick of %f km/s\n", sqrt(vs[0]*vs[0]+vs[1]*vs[1]+vs[2]*vs[2])); @@ -1267,10 +1273,10 @@ void handle_bse_outcome(long k, long kb, double *vs, double tphysf, int kprev0, { if (knew) { //bcm_boom_search(knew, vs, getCMCvalues); - pulsar_write(knew, VKO); + pulsar_write(knew, *VKO); } else { //bcm_boom_search(k, vs, getCMCvalues); - pulsar_write(k, VKO); + pulsar_write(k, *VKO); } }