From 7bcc8bb00ad2dbb834ec9d55767620b349970e14 Mon Sep 17 00:00:00 2001 From: Vera Delfavero Date: Sun, 20 Apr 2025 14:16:36 -0400 Subject: [PATCH 01/19] bugfix (issue 698) --- src/cosmic/src/comenv.f | 55 ++++++++++++++++++++++++++--------------- 1 file changed, 35 insertions(+), 20 deletions(-) diff --git a/src/cosmic/src/comenv.f b/src/cosmic/src/comenv.f index 6062b7365..80cdab0db 100644 --- a/src/cosmic/src/comenv.f +++ b/src/cosmic/src/comenv.f @@ -21,6 +21,21 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, * Update : P. D. Kiel (for ECSN, fallback and bugs) * Date : cmc version mid 2010 * +* Update : V. E. Delfavero (bugs in bpp reporting) +* Date : 20th, April 2025 +* +* Note on indexing binary components: +* +* M01, M1, MC1, AJ1, JSPIN1, KW1, formation1, bhspin1 are selected with star1 +* M02, M2, MC2, AJ2, JSPIN2, KW2, formation2, bhspin1 are selected with star2 +* +* deltam_1 and deltam_2 are calculated before being passed into comenv +* +* teff, radc, and ospin are calculated inside comenv (or a function) +* +* Other quantities are passed as vectors. +* These are lumin, menv_bpp, tms, rad, renv, bacc, tacc, epoch, and B_0. +* tms and rad are immediately separated into tms1/2_bpp and rad1/2_bpp * INTEGER KW1,KW2,KW,KW1i,KW2i,snp INTEGER star1,star2 @@ -305,11 +320,11 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, & rrl1_bpp,rrl2_bpp, & aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp, & massc1_bpp,massc2_bpp,rad1_bpp,rad2_bpp, - & M02,M01,lumin(2),lumin(1),teff2,teff1, - & RC2,RC1,menv_bpp(2),menv_bpp(1),renv_bpp(2), - & renv_bpp(1),OSPIN2,OSPIN1,B_0(2),B_0(1), - & bacc(2),bacc(1),tacc(2),tacc(1),epoch(2), - & epoch(1),bhspin2,bhspin1, + & M02,M01,lumin(1),lumin(2),teff1,teff2, + & RC2,RC1,menv_bpp(1),menv_bpp(2),renv_bpp(1), + & renv_bpp(2),OSPIN2,OSPIN1,B_0(1),B_0(2), + & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), + & epoch(2),bhspin2,bhspin1, & deltam_2,deltam_1,formation2,formation1, & binstate,mergertype,'bpp') else @@ -610,11 +625,11 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, & rrl1_bpp,rrl2_bpp, & aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp, & massc1_bpp,massc2_bpp,rad1_bpp,rad2_bpp, - & M02,M01,lumin(2),lumin(1),teff2,teff1, - & RC2,RC1,menv_bpp(2),menv_bpp(1),renv_bpp(2), - & renv_bpp(1),OSPIN2,OSPIN1,B_0(2),B_0(1), - & bacc(2),bacc(1),tacc(2),tacc(1),epoch(2), - & epoch(1),bhspin2,bhspin1, + & M02,M01,lumin(1),lumin(2),teff1,teff2, + & RC2,RC1,menv_bpp(1),menv_bpp(2),renv_bpp(1), + & renv_bpp(2),OSPIN2,OSPIN1,B_0(1),B_0(2), + & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), + & epoch(2),bhspin2,bhspin1, & deltam_2,deltam_1,formation2,formation1, & binstate,mergertype,'bpp') else @@ -779,11 +794,11 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, & rrl1_bpp,rrl2_bpp, & aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp, & massc1_bpp,massc2_bpp,rad1_bpp,rad2_bpp, - & M02,M01,lumin(2),lumin(1),teff2,teff1, - & RC2,RC1,menv_bpp(2),menv_bpp(1),renv_bpp(2), - & renv_bpp(1),OSPIN2,OSPIN1,B_0(2),B_0(1), - & bacc(2),bacc(1),tacc(2),tacc(1),epoch(2), - & epoch(1),bhspin2,bhspin1, + & M02,M01,lumin(1),lumin(2),teff1,teff2, + & RC2,RC1,menv_bpp(1),menv_bpp(2),renv_bpp(1), + & renv_bpp(2),OSPIN2,OSPIN1,B_0(1),B_0(2), + & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), + & epoch(2),bhspin2,bhspin1, & deltam_2,deltam_1,formation2,formation1, & binstate,mergertype,'bpp') else @@ -1011,11 +1026,11 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, & rrl1_bpp,rrl2_bpp, & aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp, & massc1_bpp,massc2_bpp,rad1_bpp,rad2_bpp, - & M02,M01,lumin(2),lumin(1),teff2,teff1, - & RC2,RC1,menv_bpp(2),menv_bpp(1),renv_bpp(2), - & renv_bpp(1),OSPIN2,OSPIN1,B_0(2),B_0(1), - & bacc(2),bacc(1),tacc(2),tacc(1),epoch(2), - & epoch(1),bhspin2,bhspin1, + & M02,M01,lumin(1),lumin(2),teff1,teff2, + & RC2,RC1,menv_bpp(1),menv_bpp(2),renv_bpp(1), + & renv_bpp(2),OSPIN2,OSPIN1,B_0(1),B_0(2), + & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), + & epoch(2),bhspin2,bhspin1, & deltam_2,deltam_1,formation2,formation1, & binstate,mergertype,'bpp') else From 4318963fe15eba6f7640c856fb592ea6f8035bcf Mon Sep 17 00:00:00 2001 From: Vera Delfavero Date: Mon, 21 Apr 2025 10:21:29 -0400 Subject: [PATCH 02/19] Added preSN mass tracking --- src/cosmic/src/assign_remnant.f | 6 ++++ src/cosmic/src/comenv.f | 58 +++++++++++++++++++-------------- src/cosmic/src/evolv2.f | 14 ++++++++ 3 files changed, 54 insertions(+), 24 deletions(-) diff --git a/src/cosmic/src/assign_remnant.f b/src/cosmic/src/assign_remnant.f index 858cde738..68a3464d9 100644 --- a/src/cosmic/src/assign_remnant.f +++ b/src/cosmic/src/assign_remnant.f @@ -4,6 +4,8 @@ SUBROUTINE assign_remnant(zpars,mc,mcbagb,mass,mt,kw,bhspin,kidx) INCLUDE 'const_bse.h' common /fall/fallback + REAL*8 mass_preSN, mHe_preSN, massc_preSN + COMMON mass_preSN, mHe_preSN, massc_preSN REAL*8 fallback REAL ran3 EXTERNAL ran3 @@ -70,6 +72,10 @@ SUBROUTINE assign_remnant(zpars,mc,mcbagb,mass,mt,kw,bhspin,kidx) mass = mt * else +* Store values in common block + mass_preSN = mt + mHe_preSN = mcbagb - mc + massc_preSN = mc if(ecsn.gt.0.d0.and.mcbagb.lt.ecsn_mlow)then * * Star is not massive enough to ignite C burning. diff --git a/src/cosmic/src/comenv.f b/src/cosmic/src/comenv.f index 80cdab0db..25b17eb59 100644 --- a/src/cosmic/src/comenv.f +++ b/src/cosmic/src/comenv.f @@ -56,6 +56,8 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, REAL*8 bhspin1,bhspin2 REAL*8 deltam_1,deltam_2 common /fall/fallback + REAL*8 mass_preSN, mHe_preSN, massc_preSN + COMMON mass_preSN, mHe_preSN, massc_preSN INTEGER formation1,formation2 REAL*8 sigmahold REAL*8 AURSUN,K3 @@ -319,9 +321,10 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, & kstar2_bpp,SEP_postCE,TB,ECC, & rrl1_bpp,rrl2_bpp, & aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp, - & massc1_bpp,massc2_bpp,rad1_bpp,rad2_bpp, - & M02,M01,lumin(1),lumin(2),teff1,teff2, - & RC2,RC1,menv_bpp(1),menv_bpp(2),renv_bpp(1), + & massc1_bpp,massc_preSN,rad1_bpp,rad2_bpp, + & M02,mass_preSN,lumin(1),lumin(2), + & teff1,teff2, + & RC2,RC1,menv_bpp(1),mHe_preSN,renv_bpp(1), & renv_bpp(2),OSPIN2,OSPIN1,B_0(1),B_0(2), & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin2,bhspin1, @@ -337,9 +340,10 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, & kstar2_bpp,SEP_postCE,TB,ECC, & rrl1_bpp,rrl2_bpp, & aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp, - & massc1_bpp,massc2_bpp,rad1_bpp,rad2_bpp, - & M01,M02,lumin(1),lumin(2),teff1,teff2, - & RC1,RC2,menv_bpp(1),menv_bpp(2),renv_bpp(1), + & massc_preSN,massc2_bpp,rad1_bpp,rad2_bpp, + & mass_preSN,M02,lumin(1),lumin(2), + & teff1,teff2, + & RC1,RC2,mHe_preSN,menv_bpp(2),renv_bpp(1), & renv_bpp(2),OSPIN1,OSPIN2,B_0(1),B_0(2), & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin1,bhspin2, @@ -624,9 +628,10 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, & kstar2_bpp,SEP_postCE,TB,ECC, & rrl1_bpp,rrl2_bpp, & aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp, - & massc1_bpp,massc2_bpp,rad1_bpp,rad2_bpp, - & M02,M01,lumin(1),lumin(2),teff1,teff2, - & RC2,RC1,menv_bpp(1),menv_bpp(2),renv_bpp(1), + & massc1_bpp,massc_preSN,rad1_bpp,rad2_bpp, + & M02,mass_preSN,lumin(1),lumin(2), + & teff1,teff2, + & RC2,RC1,menv_bpp(1),mHe_preSN,renv_bpp(1), & renv_bpp(2),OSPIN2,OSPIN1,B_0(1),B_0(2), & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin2,bhspin1, @@ -642,9 +647,10 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, & kstar2_bpp,SEP_postCE,TB,ECC, & rrl1_bpp,rrl2_bpp, & aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp, - & massc1_bpp,massc2_bpp,rad1_bpp,rad2_bpp, - & M01,M02,lumin(1),lumin(2),teff1,teff2, - & RC1,RC2,menv_bpp(1),menv_bpp(2),renv_bpp(1), + & massc_preSN,massc2_bpp,rad1_bpp,rad2_bpp, + & mass_preSN,M02,lumin(1),lumin(2), + & teff1,teff2, + & RC1,RC2,mHe_preSN,menv_bpp(2),renv_bpp(1), & renv_bpp(2),OSPIN1,OSPIN2,B_0(1),B_0(2), & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin1,bhspin2, @@ -793,9 +799,10 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, & kstar2_bpp,SEP_postCE,TB,ECC, & rrl1_bpp,rrl2_bpp, & aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp, - & massc1_bpp,massc2_bpp,rad1_bpp,rad2_bpp, - & M02,M01,lumin(1),lumin(2),teff1,teff2, - & RC2,RC1,menv_bpp(1),menv_bpp(2),renv_bpp(1), + & massc_preSN,massc2_bpp,rad1_bpp,rad2_bpp, + & mass_preSN,M01,lumin(1),lumin(2), + & teff1,teff2, + & RC2,RC1,mHe_preSN,menv_bpp(2),renv_bpp(1), & renv_bpp(2),OSPIN2,OSPIN1,B_0(1),B_0(2), & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin2,bhspin1, @@ -811,9 +818,10 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, & kstar2_bpp,SEP_postCE,TB,ECC, & rrl1_bpp,rrl2_bpp, & aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp, - & massc1_bpp,massc2_bpp,rad1_bpp,rad2_bpp, - & M01,M02,lumin(1),lumin(2),teff1,teff2, - & RC1,RC2,menv_bpp(1),menv_bpp(2),renv_bpp(1), + & massc1_bpp,massc_preSN,rad1_bpp,rad2_bpp, + & M01,mass_preSN,lumin(1),lumin(2), + & teff1,teff2, + & RC1,RC2,menv_bpp(1),mHe_preSN,renv_bpp(1), & renv_bpp(2),OSPIN1,OSPIN2,B_0(1),B_0(2), & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin1,bhspin2, @@ -1025,9 +1033,10 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, & kstar2_bpp,-1.d0,TB,0.d0, & rrl1_bpp,rrl2_bpp, & aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp, - & massc1_bpp,massc2_bpp,rad1_bpp,rad2_bpp, - & M02,M01,lumin(1),lumin(2),teff1,teff2, - & RC2,RC1,menv_bpp(1),menv_bpp(2),renv_bpp(1), + & massc1_bpp,massc_preSN,rad1_bpp,rad2_bpp, + & M02,mass_preSN,lumin(1),lumin(2), + & teff1,teff2, + & RC2,RC1,menv_bpp(1),mHe_preSN,renv_bpp(1), & renv_bpp(2),OSPIN2,OSPIN1,B_0(1),B_0(2), & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin2,bhspin1, @@ -1043,9 +1052,10 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, & kstar2_bpp,-1.d0,TB,0.d0, & rrl1_bpp,rrl2_bpp, & aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp, - & massc1_bpp,massc2_bpp,rad1_bpp,rad2_bpp, - & M01,M02,lumin(1),lumin(2),teff1,teff2, - & RC1,RC2,menv_bpp(1),menv_bpp(2),renv_bpp(1), + & massc_preSN,massc2_bpp,rad1_bpp,rad2_bpp, + & mass_preSN,M02,lumin(1),lumin(2), + & teff1,teff2, + & RC1,RC2,mHe_preSN,menv_bpp(2),renv_bpp(1), & renv_bpp(2),OSPIN1,OSPIN2,B_0(1),B_0(2), & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin1,bhspin2, diff --git a/src/cosmic/src/evolv2.f b/src/cosmic/src/evolv2.f index 083b6c7b8..d022a36da 100644 --- a/src/cosmic/src/evolv2.f +++ b/src/cosmic/src/evolv2.f @@ -195,6 +195,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, REAL*8 deltam1_bcm,deltam2_bcm,b01_bcm,b02_bcm REAL*8 B(2),Bbot,omdot,b_mdot,b_mdot_lim,evolve_type COMMON /fall/fallback + REAL*8 mass_preSN, mHe_preSN, massc_preSN + COMMON mass_preSN, mHe_preSN, massc_preSN REAL ran3 EXTERNAL ran3 * @@ -1328,6 +1330,10 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, else b02_bcm = B(2) endif +* Load preSN values for the SN writetab + mass0(k) = mass_preSN + menv(k) = mHe_preSN + massc(k) = massc_preSN CALL writetab(jp,tphys,evolve_type, & mass(1),mass(2),kstar(1),kstar(2), & sep,tb,ecc,rrl1,rrl2, @@ -1368,6 +1374,10 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, b02_bcm = B(2) endif +* Load preSN values for the SN writetab + mass0(k) = mass_preSN + menv(k) = mHe_preSN + massc(k) = massc_preSN CALL writetab(jp,tphys,evolve_type, & mass(1),mass(2),kstar(1),kstar(2), & sep,tb,ecc,rrl1,rrl2, @@ -3629,6 +3639,10 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, else b02_bcm = B(2) endif +* Load preSN values for the SN writetab + mass0(k) = mass_preSN + menv(k) = mHe_preSN + massc(k) = massc_preSN CALL writetab(jp,tphys,evolve_type, & mass(1),mass(2),kstar(1),kstar(2), & sep,tb,ecc,rrl1,rrl2, From 1aa12b9e95a9eb5629243772ca19d24395f7a075 Mon Sep 17 00:00:00 2001 From: Vera Delfavero Date: Tue, 20 May 2025 11:47:05 -0400 Subject: [PATCH 03/19] fixed M0 write --- src/cosmic/src/evolv2.f | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/cosmic/src/evolv2.f b/src/cosmic/src/evolv2.f index d022a36da..f7038eb3a 100644 --- a/src/cosmic/src/evolv2.f +++ b/src/cosmic/src/evolv2.f @@ -1376,6 +1376,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, * Load preSN values for the SN writetab mass0(k) = mass_preSN + m0 = mass_preSN menv(k) = mHe_preSN massc(k) = massc_preSN CALL writetab(jp,tphys,evolve_type, @@ -3641,6 +3642,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, endif * Load preSN values for the SN writetab mass0(k) = mass_preSN + m0 = mass_preSN menv(k) = mHe_preSN massc(k) = massc_preSN CALL writetab(jp,tphys,evolve_type, From f6f178fc79a4138e612558e66f1598068c65e11b Mon Sep 17 00:00:00 2001 From: Vera Delfavero Date: Wed, 21 May 2025 11:04:33 -0400 Subject: [PATCH 04/19] Missed one --- src/cosmic/src/evolv2.f | 1 + 1 file changed, 1 insertion(+) diff --git a/src/cosmic/src/evolv2.f b/src/cosmic/src/evolv2.f index f7038eb3a..af7f1dd6a 100644 --- a/src/cosmic/src/evolv2.f +++ b/src/cosmic/src/evolv2.f @@ -1332,6 +1332,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, endif * Load preSN values for the SN writetab mass0(k) = mass_preSN + m0 = mass_preSN menv(k) = mHe_preSN massc(k) = massc_preSN CALL writetab(jp,tphys,evolve_type, From f24bd0f27c36c7fd03c5467ec652e0ed36a34477 Mon Sep 17 00:00:00 2001 From: Katie Breivik Date: Mon, 9 Jun 2025 16:18:02 -0400 Subject: [PATCH 05/19] this fixes the test data --- src/cosmic/tests/data/unit_tests_results.hdf5 | Bin 6470124 -> 7518844 bytes 1 file changed, 0 insertions(+), 0 deletions(-) diff --git a/src/cosmic/tests/data/unit_tests_results.hdf5 b/src/cosmic/tests/data/unit_tests_results.hdf5 index f08ffb33f1c56abed87be610741dbc3f78787195..eea66b61f3be22bfd2db8fee9c6671d3933eacf0 100644 GIT binary patch delta 3150 zcmc(hdsI_L8o)CNgj^5^8c9S#0!9tUBcY1;h<6f2a4m``hs76^l_Q=aAV$~M9@#x8 z91FYP)K2RQkRXbv#JEjRx<1yfRoGqephd7pi|uYb!q(bosYTr{j~ndc&pl`FIlntI zcfR=^GjnIYfp7N+hV9n{)oQ_xb0sxXcA+h1vkzZe!E-$+5?6_M!8`<2nftq#CwA5< zVqG7w!~ZA4z(I60rRGB#2n#Hf^96HZ23$Tvc%Uw*V}Sf|ZU zV}458jVOe^+=WDzC{X|kQoD_ea+hkdfLXxj1yW4|U*=TW(8{3Kl~6?^Kg4p&_Yp&K zBw}(=B0}mKvx49x(V4ONQgUOyS50Lo!}D@!Of;LnLyk^zpDz-m@OX4*0g9ojLZr4l zj8^1JG_HX^j*DnHBBPV`xLS2dIS5ug*CI#`WSH&^ zE_uHn$^RNzs+Sam(#m=yQVg=`9cg7T*n4}5eD;nm{9@UPe$@{FxcZ&TcgC6SlR$Ym zd9zSSrq0aq$vo9r`Tg=*!n}|J6Oy2IP^rNAaYgw*5OccY-*xFYJNNnms2= zyK->OgF0jP`Dy;H{>fMp^V&;ZE>Q4w^pGhng`L|T$K#?9v$Y({# z-p1nmtsCi*IwU14iz|t}xQCV%A~7x4ghVwLWsJXDg+KjWA(GO{-AF{*bt1p@c_nM} zO0tbxi}Nh^Ht3%_X9Tutg|(#sJ)TO6LuG+i=A6BhFG*&(P-;V$1aYA^-$V;@IzGYG z1j%dZ-6Kdwy&8}@{MC7V*r)k1c)L|;QLl)>v1DNIHnhKbv3%))Ra#t6GaArN=7YNf ztN*=wFRrDZH6ZvkXqXAbswX8MUTrVz!GDl<#Z?%d;+&pmy*YpFV_#+*T=x)fWtEAu zHW$Qau(Z;Iicvj%VM1|YUCZ6&M7PsQ{c*4d?}|`i zd&DzrW*4C~m-^T$`vmj9hOe^wPJ6fAPRwHI7mY~6w6|8|b}p;NHYz`gCNsAOpXCkC zJ%;~AQ=vJocZ=qxrKt*^|LD^@_jdH@ez>9^oxbV!?0dB8D9TRcR(E5`_%RcwMd2%A zo&ULrUTm|B+FooN!oOMjq^s`8FPNU)gp~O!RiC_M+x#)N+V=1z)~fS$xLwYzy-C`MH zY3OkzW3p?bwzX{QF*H%#aWtKV9z*ItZdbZ3ke#Y=1Wypt|F?UYQv0f=*;_NTp8sa? zP*$-|d=z*X3A<+SI8vY<>fekM;oKUyCGg7`o7<)=-^ZGUMz?-&<^i5hQ=5@gl)zh- zdZYO+OO$dY`kGcXqZG8C#7R1;zuS{NVk_wJ*oNO~UA8nd>4ei(9uz;bb^T*` zRzW{u<3CQnkWv2-`_aY|h(&h#<%tpXPxkiu{$25gCzPJh>VpB>BAuaV=%E&;`j1+W zR2#SWm4!CG<5_qvq?WG4_W#UU$C}oxQ^z(u#OF3_hvK5zI~6O! z^4a(I{`oKR4~6=yuwAxppJ2U>t~v#iQ`g@5y?@Fxmit#!y`Oh&>pO(Ip> zT4jb)Q9MN?<^p_OH#M~)CG4GF)XbgE!`HW-#E?Jd4-^Cehoq zDA>|D3NEc>Cceyh>XOIPMYL=;55x@YxfFCdl@CO2bKtK!xS%s$2&s_-}IIBkG(M_X)F@OXZ3j_mF zKn9EhLV)o=C?E&I00p1~!mZIw5sBhaE8EhY6I%*zvQjX$++>XaCIBj6A}|S<45)!f zAPSfQL<2EEED#5104)#?ya70LO#>2u=|Cbd1DFZS0%ij`U=F|n7|>fyR)ZQ5rBEgq zuo=iK;i~d^Z3bX2Fb_xq<^#z<3a|k96|fL^6G#OX0gHh&U2m<3b1ce$zcyHnmIBSoTjX4m7p*_YHgj(uax6^8? zgJXXj;M!eRe-I&vC_h4D8&RxXhhtXeprw`y=M=H5WrYs))iK9p4q||Ty7J68!g)(|Y8U{vsTk6Axo(8#CBdu<~ z&pOM&k3Bp9a%-gA9XzA#5+XN8jG~!brdXq3+Iq3Ci@SPYYrTa7PUK=A*D!U1*x&l` zqOUmG>m>(Wb$c8oY4W#JMxPoNo!Amn-5$W)N>IMgr5NS}Xs_!z>tw9p;-W%1YZmbW zI-Stz8RAGOcEk{RB{5S>#P-y4bG46Hp!z)pdwmy3>G_l}t6KMwjC0-?LVvn$DKpiI z0~RHe*+neyjk~kRUSJPj1}U3<8RT%j~X*!58S>RZ9Bgx^PaU zvrka6g$Y?C52!NjpnbY2tq`za*;Yxq2gfNL^B4ci{GWie>hWTKTIepN?Vv0WGf@-Q>A)P(U#J5`g!ts>j%YB zuzvpZ+6Qx&W+PV#K14#Fdx z-_Bgn2Yh0TAurlA*6>X+B$kd(T*vFNf8_n}$mV!rVxjf$BKPLd_?9nzSwoJomUox;BZUL^ZxsG^ke|spx?_wy&uetayNg%m@XqAy z`(`MU)yaQ#tJ`gumIg#{sm()qC%uZuW;S7lQipf1`2;M#0epR_5ZwN!>2Cf|n`_5) zUG9S1{sz!;?^!B~dVW+zZ0uqKgmE8D+<#(X-6_(|Ch&6E$aC5ct&Ojd3l~56^s`?K zrrhq*E3=>7#~)@fC!rwgX7|CwadiZFX4mqip^1@>#0v(x6)8*3ued=jTHStLRz5|x zZokuKsvG3L+fpoHT_-_f49?!_^9rYOsEI^VL66&ByI>+C+T(EeeGPGL8)w$-$F~Pw zzbe-&o~h_c3;K+>W$ho$?lzOtw*4XHm|f9i`N3DaA2=5-YJXSJ<({mK z+|#^gz|h1>8lk}5xfku?)H~Sl|25Te3O0J*+uXSJ`JDLqgJ9qBbZ3YGs#sbR1VSHs z0z;DX(KwC-=}+3-5ta7_Z=O=o zAOnuF*wf(cT~+wzp@D*hte%sq`k%$q--Q~0umh*(vNm7{`uD8n%j(?X5z4Q;oP8p2K=jzy6>R@VFu8#%V!`k(^-n??T7k83X}_G zH_TXB!z=Uqz`0n~r|bMWc?WOawj(^YevREp~-CNJIm}WG)g*`SKx>!jwo=N)} zeN=aHj(n(HaB;yc=Rw6{+91%w`H&Z6d2t9~ z``bV>yQ#iQ*ZDI&z0sU;Rxy$?!baP$yf3j(3(Bili3O7Bkb5e-W&v;8crCmN)=tGC ziM4lE_?NVo*Kp3h8@15G%i{0`$pjuAuX&)2KH)^0a^+}RbR-i4e zf%ZO23=i0~I(lyb&M77g>h*A zfrthPvPCzp&?y%8oHNWlhrJ>Z<4xk#>;>D3@QQ1;6-Y1=f`lU9L&A{Ru?Qp*i9({0 z7$g>nL$pXdvJzQ^Bp`{%Bgkq*ha@3ukYr>nl7gfnJVFpXyJj;yoc*W6wsmBiDg<;R zpc?_*2y8v)%2=te*{0=f}wTAN^t|8j;+mw=*1E2qO7 Date: Tue, 10 Jun 2025 08:17:36 -0400 Subject: [PATCH 06/19] hopefully this matches the python versions on the github action runners --- meson.build | 2 +- src/cosmic/tests/data/unit_tests_results.hdf5 | Bin 7518844 -> 8567500 bytes 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/meson.build b/meson.build index 7e0a669c9..f21768d98 100644 --- a/meson.build +++ b/meson.build @@ -1,7 +1,7 @@ project('cosmic', 'c', 'fortran', - version : '3.5.0', + version : '3.6.0', default_options: ['warning_level=0', 'optimization=3'], ) diff --git a/src/cosmic/tests/data/unit_tests_results.hdf5 b/src/cosmic/tests/data/unit_tests_results.hdf5 index eea66b61f3be22bfd2db8fee9c6671d3933eacf0..138c564f5ad91af20ccfcac4f9c2c76a64bc345c 100644 GIT binary patch delta 1052 zcmb``O-~a+90u^tmX_HTu`VqjSOipFtI+ZyDoQ~sRnQiEMNpveNKLHaU_zqNm{v~3 zNgQH~yf!3;1jruRV*CJT6MFE-%{fYLN(!sXq*UUZYnmj|eBjK{} zy=1`8rXEuE)J#>kIV!U;tA(Yll$0{3r=-A7v$gm$&q*&yRyOZ^#Tr`L*>_s0(&8Pv z{t*IwZW~eeTAPX1YHGKea^1Vxp01wB?1|4kOqk-C)-KO^#tEh%r0RQ}a&Ev7c^(~z zM+Zo>EYpM6c^O`vAo0}+skKEb-sgTyShGy|{RL@4d*!a4zr@0~7bTmVHf8G_*0Hid%t>d$Zv;UsjxDL4&h zpcBr*Ip_j6bVCoEhYR3=i_i;RxCED>5BlK>Tm>KaVSuL+gR}CtWI$V;NH~}rgllje zZoo~r1-D@c?!YjNz+Jcpqc8^JJeUkTVQvq;iR7Td6x)NaK%zjRK%zjRK%zjRK%zjR zK%zjRK%zjRK%zjRK%zjRNCqg2x=g%gNh)UEZVP^)%KRAWo!SiW(Me)c)1fcF|L<-& F^bbb#w=4hv delta 941 zcmb`_-Aj{E90zcJ&$G?*OsmJatuo(xQFEo6FIic8vX=JdroEb$IU<4-YNQx=A%w2H z+BndOh_qCUFmgN>+lL zjF7^Yt&B>;kpixHn0P*EQH5dFf0u`xx5fQoD%6?%Yo5G@BVNr$Vl3ve>S|3cB|{ap zdFG9^@2QBpp*J=0YVE`=o{{A8sHcliuV0rNoyx|6HhLIpiH2Gzl#_@c^~b^KZi-HK zOVw4$9KG5D;&Y6Nn=hsAWO?p}KeSmaypnA4O586#-k~DF#w9BsGn&QlY``w+$ED2F zR^@%rG$DZu1~7sN(t&{jX2^g{*aVwl3#gD4HBDHakmQJo7n2tMMoMKAbWB+w8?0ah zJ8Xq*upM$>2joH?;<=Q8 zm@D~ReEgox%@@t8Sc$rtq&dcOs!j&yV=kzIYN&x)H~@8U5Dvj%sD~qP6x^V}F=&9} zZ~{D#xU&(O;1ryOGjJBpK{I&42Q6Yg*7{Ve#{HIHU%VC0!v(kqm*6s7fi}1b?a%?& x;5u}|4d@bmaetp!?;({hn|*vaL* Date: Tue, 10 Jun 2025 15:42:01 -0400 Subject: [PATCH 07/19] I think we are running into a numpy change specifically --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index b7db21f3c..578a9fffb 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,5 @@ scipy >= 0.12.1 -numpy == 1.23.5 +numpy >= 1.26.0 astropy >= 1.1.1 configparser tqdm >= 4.0 From ff761e1ba01ec074992a4c0152ad285fd8024d26 Mon Sep 17 00:00:00 2001 From: Katie Breivik Date: Thu, 26 Jun 2025 14:02:59 -0400 Subject: [PATCH 08/19] intermediate updates --- src/cosmic/src/bpp_array.f | 89 ++++++++++++++++++++------------------ src/cosmic/src/const_bse.h | 4 +- src/cosmic/src/hrdiag.f | 2 +- 3 files changed, 51 insertions(+), 44 deletions(-) diff --git a/src/cosmic/src/bpp_array.f b/src/cosmic/src/bpp_array.f index 5ea53257e..b71ddd88a 100644 --- a/src/cosmic/src/bpp_array.f +++ b/src/cosmic/src/bpp_array.f @@ -3,26 +3,30 @@ SUBROUTINE WRITETAB(jp,tphys,evolve_type, & mass1,mass2,kstar1,kstar2,sep, & tb,ecc,rrl1,rrl2, & aj1,aj2,tms1,tms2, - & massc1,massc2,rad1,rad2, + & massc_he_1,massc_he_2, + & massc_co_1,massc_co_2, rad1,rad2, & mass0_1,mass0_2,lumin1,lumin2, & teff1,teff2,radc1,radc2,menv1, & menv2,renv1,renv2,ospin1,ospin2, & b_0_1,b_0_2,bacc1,bacc2,tacc1,tacc2, & epoch1,epoch2,bhspin1,bhspin2, & deltam_1,deltam_2,SN_1,SN_2, - & bin_state,merger_type,tabname) + & bin_state,merger_type,metallicity, + & tabname) IMPLICIT NONE INCLUDE 'const_bse.h' * * Write results to bpp or bcm array. * -* Author : Scott Coughlin, Tom Wagg -* Date : 12th March 2019, September 2024 +* Author : Scott Coughlin, Tom Wagg, Katie Breivik +* Date : 12th March 2019, September 2024, June 2025 * - REAL*8 mass1,mass2 + REAL*8 mass1,mass2,metallicity REAL*8 evolve_type,sep,tb,ecc,tphys,rrl1,rrl2 - REAL*8 aj1,aj2,tms1,tms2,massc1,massc2,rad1,rad2 + REAL*8 aj1,aj2,tms1,tms2,rad1,rad2 + REAL*8 massc_he_1,massc_he_2 + REAL*8 massc_co_1,massc_co_2 REAL*8 mass0_1,mass0_2,lumin1,lumin2,radc1,radc2 REAL*8 menv1,menv2,renv1,renv2,ospin1,ospin2 REAL*8 b_0_1,b_0_2,bacc1,bacc2,tacc1,tacc2,epoch1,epoch2 @@ -33,7 +37,7 @@ SUBROUTINE WRITETAB(jp,tphys,evolve_type, INTEGER jp, col_ind INTEGER kstar1,kstar2 REAL*8 yeardy,aursun,rsunau - REAL*8 all_cols(49) + REAL*8 all_cols(52) CHARACTER*3 tabname PARAMETER(yeardy=365.24d0,aursun=214.95d0) @@ -60,40 +64,43 @@ SUBROUTINE WRITETAB(jp,tphys,evolve_type, all_cols(13) = aj2 all_cols(14) = tms1 all_cols(15) = tms2 - all_cols(16) = massc1 - all_cols(17) = massc2 - all_cols(18) = rad1 - all_cols(19) = rad2 - all_cols(20) = mass0_1 - all_cols(21) = mass0_2 - all_cols(22) = lumin1 - all_cols(23) = lumin2 - all_cols(24) = teff1 - all_cols(25) = teff2 - all_cols(26) = radc1 - all_cols(27) = radc2 - all_cols(28) = menv1 - all_cols(29) = menv2 - all_cols(30) = renv1 - all_cols(31) = renv2 - all_cols(32) = ospin1 - all_cols(33) = ospin2 - all_cols(34) = b_0_1 - all_cols(35) = b_0_2 - all_cols(36) = bacc1 - all_cols(37) = bacc2 - all_cols(38) = tacc1 - all_cols(39) = tacc2 - all_cols(40) = epoch1 - all_cols(41) = epoch2 - all_cols(42) = bhspin1 - all_cols(43) = bhspin2 - all_cols(44) = deltam_1 - all_cols(45) = deltam_2 - all_cols(46) = float(SN_1) - all_cols(47) = float(SN_2) - all_cols(48) = bin_state - all_cols(49) = merger_type + all_cols(16) = massc_he_1 + all_cols(17) = massc_he_2 + all_cols(18) = massc_co_1 + all_cols(19) = massc_co_2 + all_cols(20) = rad1 + all_cols(21) = rad2 + all_cols(22) = mass0_1 + all_cols(23) = mass0_2 + all_cols(24) = lumin1 + all_cols(25) = lumin2 + all_cols(26) = teff1 + all_cols(27) = teff2 + all_cols(28) = radc1 + all_cols(29) = radc2 + all_cols(30) = menv1 + all_cols(31) = menv2 + all_cols(32) = renv1 + all_cols(33) = renv2 + all_cols(34) = ospin1 + all_cols(35) = ospin2 + all_cols(36) = b_0_1 + all_cols(37) = b_0_2 + all_cols(38) = bacc1 + all_cols(39) = bacc2 + all_cols(40) = tacc1 + all_cols(41) = tacc2 + all_cols(42) = epoch1 + all_cols(43) = epoch2 + all_cols(44) = bhspin1 + all_cols(45) = bhspin2 + all_cols(46) = deltam_1 + all_cols(47) = deltam_2 + all_cols(48) = float(SN_1) + all_cols(49) = float(SN_2) + all_cols(50) = bin_state + all_cols(51) = merger_type + all_cols(52) = metallicity * check which table we are writing to and write the appropriate columns if (tabname .eq. 'bpp') then diff --git a/src/cosmic/src/const_bse.h b/src/cosmic/src/const_bse.h index 4595c1809..e760b891b 100644 --- a/src/cosmic/src/const_bse.h +++ b/src/cosmic/src/const_bse.h @@ -55,9 +55,9 @@ COMMON /TSTEPC/ dmmax,drmax REAL*8 scm(50000,14),spp(20,3) COMMON /SINGLE/ scm,spp - REAL*8 bcm(50000,49),bpp(1000,49) + REAL*8 bcm(50000,52),bpp(1000,52) COMMON /BINARY/ bcm,bpp INTEGER n_col_bpp, n_col_bcm - INTEGER col_inds_bpp(49), col_inds_bcm(49) + INTEGER col_inds_bpp(52), col_inds_bcm(52) COMMON /COL/ n_col_bpp,col_inds_bpp,n_col_bcm,col_inds_bcm * diff --git a/src/cosmic/src/hrdiag.f b/src/cosmic/src/hrdiag.f index e554ef695..dc19af197 100644 --- a/src/cosmic/src/hrdiag.f +++ b/src/cosmic/src/hrdiag.f @@ -27,7 +27,7 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, * real*8 mass,aj,mt,tm,tn,tscls(20),lums(10),GB(10),zpars(20),met real*8 bhspin - real*8 r,lum,mc,rc,menv,renv,k2 + real*8 r,lum,mc_he,mc_co,rc,menv,renv,k2 real*8 mch,mlp,tiny * parameter(mch=1.44d0,mlp=12.d0,tiny=1.0d-14) parameter(mlp=12.d0,tiny=1.0d-14) From 4868f4436772a33d6dddb17a09eae09b7b94eb5e Mon Sep 17 00:00:00 2001 From: Katie Breivik Date: Mon, 11 Aug 2025 11:41:40 -0400 Subject: [PATCH 09/19] lots of debugging still to do but this compiles --- src/cosmic/evolve.py | 19 ++++-- src/cosmic/src/const_bse.h | 3 +- src/cosmic/src/evolv2.f | 134 +++++++++++++++++++++++-------------- src/cosmic/src/hrdiag.f | 40 ++++++++++- 4 files changed, 135 insertions(+), 61 deletions(-) diff --git a/src/cosmic/evolve.py b/src/cosmic/evolve.py index 0d46e0254..6c9d105e9 100644 --- a/src/cosmic/evolve.py +++ b/src/cosmic/evolve.py @@ -49,30 +49,31 @@ ALL_COLUMNS = ['tphys', 'mass_1', 'mass_2', 'kstar_1', 'kstar_2', 'sep', 'porb', 'ecc', 'RRLO_1', 'RRLO_2', 'evol_type', 'aj_1', 'aj_2', 'tms_1', - 'tms_2', 'massc_1', 'massc_2', 'rad_1', 'rad_2', 'mass0_1', + 'tms_2', 'massc_he_1', 'massc_he_2', 'massc_co_1', 'massc_co_2', + 'rad_1', 'rad_2', 'mass0_1', 'mass0_2', 'lum_1', 'lum_2', 'teff_1', 'teff_2', 'radc_1', 'radc_2', 'menv_1', 'menv_2', 'renv_1', 'renv_2', 'omega_spin_1', 'omega_spin_2', 'B_1', 'B_2', 'bacc_1', 'bacc_2', 'tacc_1', 'tacc_2', 'epoch_1', 'epoch_2', 'bhspin_1', 'bhspin_2', - 'deltam_1', 'deltam_2', 'SN_1', 'SN_2', 'bin_state', 'merger_type'] + 'deltam_1', 'deltam_2', 'SN_1', 'SN_2', 'bin_state', 'merger_type', 'metallicity'] INTEGER_COLUMNS = ["bin_state", "bin_num", "kstar_1", "kstar_2", "SN_1", "SN_2", "evol_type"] -BPP_COLUMNS = ['tphys', 'mass_1', 'mass_2', 'kstar_1', 'kstar_2', +BPP_COLUMNS = ['tphys', 'metallicity', 'mass_1', 'mass_2', 'kstar_1', 'kstar_2', 'sep', 'porb', 'ecc', 'RRLO_1', 'RRLO_2', 'evol_type', 'aj_1', 'aj_2', 'tms_1', 'tms_2', - 'massc_1', 'massc_2', 'rad_1', 'rad_2', + 'massc_he_1', 'massc_he_2', 'massc_co_1', 'massc_co_2', 'rad_1', 'rad_2', 'mass0_1', 'mass0_2', 'lum_1', 'lum_2', 'teff_1', 'teff_2', 'radc_1', 'radc_2', 'menv_1', 'menv_2', 'renv_1', 'renv_2', 'omega_spin_1', 'omega_spin_2', 'B_1', 'B_2', 'bacc_1', 'bacc_2', 'tacc_1', 'tacc_2', 'epoch_1', 'epoch_2', 'bhspin_1', 'bhspin_2'] -BCM_COLUMNS = ['tphys', 'kstar_1', 'mass0_1', 'mass_1', 'lum_1', 'rad_1', - 'teff_1', 'massc_1', 'radc_1', 'menv_1', 'renv_1', 'epoch_1', +BCM_COLUMNS = ['tphys', 'metallicity', 'kstar_1', 'mass0_1', 'mass_1', 'lum_1', 'rad_1', + 'teff_1', 'massc_he_1', 'massc_co_1', 'radc_1', 'menv_1', 'renv_1', 'epoch_1', 'omega_spin_1', 'deltam_1', 'RRLO_1', 'kstar_2', 'mass0_2', 'mass_2', - 'lum_2', 'rad_2', 'teff_2', 'massc_2', 'radc_2', 'menv_2', + 'lum_2', 'rad_2', 'teff_2', 'massc_he_2', 'massc_co_2', 'radc_2', 'menv_2', 'renv_2', 'epoch_2', 'omega_spin_2', 'deltam_2', 'RRLO_2', 'porb', 'sep', 'ecc', 'B_1', 'B_2', 'SN_1', 'SN_2', 'bin_state', 'merger_type'] @@ -368,6 +369,8 @@ def evolve(self, initialbinarytable, pool=None, bpp_columns=None, bcm_columns=No for i in range(len(initial_conditions)): initial_conditions[i]["n_col_bpp"] = len(bpp_columns) initial_conditions[i]["col_inds_bpp"] = col_inds_bpp + print("n_col_bpp: ", initial_conditions[0]["n_col_bpp"]) + print("col_inds_bpp: ", initial_conditions[0]["col_inds_bpp"], len(initial_conditions[0]["col_inds_bpp"])) # same for bcm col_inds_bcm = np.zeros(len(ALL_COLUMNS), dtype=int) @@ -375,6 +378,8 @@ def evolve(self, initialbinarytable, pool=None, bpp_columns=None, bcm_columns=No for i in range(len(initial_conditions)): initial_conditions[i]["n_col_bcm"] = len(bcm_columns) initial_conditions[i]["col_inds_bcm"] = col_inds_bcm + print("n_col_bcm: ", initial_conditions[0]["n_col_bcm"]) + print("col_inds_bcm: ", initial_conditions[0]["col_inds_bcm"], len(initial_conditions[0]["col_inds_bcm"])) # check if a pool was passed if pool is None: diff --git a/src/cosmic/src/const_bse.h b/src/cosmic/src/const_bse.h index e760b891b..2327a95d6 100644 --- a/src/cosmic/src/const_bse.h +++ b/src/cosmic/src/const_bse.h @@ -20,7 +20,8 @@ INTEGER ceflag,cekickflag,cemergeflag,cehestarflag,ussn COMMON /CEFLAGS/ ceflag,cekickflag,cemergeflag,cehestarflag,ussn INTEGER pisn_track(2) - COMMON /TRACKERS/ pisn_track + REAL*8 mc_he(2),mc_co(2) + COMMON /TRACKERS/ pisn_track,mc_he,mc_co * REAL*8 zsun COMMON /METVARS/ zsun diff --git a/src/cosmic/src/evolv2.f b/src/cosmic/src/evolv2.f index af7f1dd6a..c48e19e37 100644 --- a/src/cosmic/src/evolv2.f +++ b/src/cosmic/src/evolv2.f @@ -414,6 +414,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, CALL star(kstar(k),mass0(k),mass(k),tm,tn,tscls,lums,GB,zpars) CALL hrdiag(mass0(k),age,mass(k),tm,tn,tscls,lums,GB,zpars, & rm,lum,kstar(k),mc,rc,me,re,k2,bhspin(k),k) + WRITE(*,*)'k=',k,' kstar=',kstar(k),'HRDIAG ran' aj(k) = age epoch(k) = tphys - age rad(k) = rm @@ -631,6 +632,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, djorb = djorb + djgr delet = delet + delet1 endif + * do 502 , k = 1,2 @@ -1215,6 +1217,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, * goto 140 endif * +* WRITE(*,*)'hrdiag call 1220' CALL star(kw,m0,mt,tm,tn,tscls,lums,GB,zpars) CALL hrdiag(m0,age,mt,tm,tn,tscls,lums,GB,zpars, & rm,lum,kw,mc,rc,me,re,k2,bhspin(k),k) @@ -1222,6 +1225,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, if(kw.ne.15)then ospin(k) = jspin(k)/(k2*(mt-mc)*rm*rm+k3*mc*rc*rc) endif +* WRITE(*,*)'hrdiag call finished' * * At this point there may have been a supernova. * @@ -1339,7 +1343,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & mass(1),mass(2),kstar(1),kstar(2), & sep,tb,ecc,rrl1,rrl2, & aj(1),aj(2),tms(1),tms(2), - & massc(1),massc(2),rad(1),rad(2), + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad(1),rad(2), & mass0(1),mass0(2),lumin(1),lumin(2), & teff1,teff2,radc(1),radc(2), & menv(1),menv(2),renv(1),renv(2), @@ -1347,7 +1352,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin(1),bhspin(2), & deltam1_bcm,deltam2_bcm,formation(1), - & formation(2),binstate,mergertype,'bpp') + & formation(2),binstate,mergertype,z,'bpp') CALL kick(kw,mass(k),mt,0.d0,0.d0,-1.d0,0.d0,vk,k, & 0.d0,fallback,sigmahold,kick_info,disrupt,bkick) @@ -1384,7 +1389,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & mass(1),mass(2),kstar(1),kstar(2), & sep,tb,ecc,rrl1,rrl2, & aj(1),aj(2),tms(1),tms(2), - & massc(1),massc(2),rad(1),rad(2), + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad(1),rad(2), & mass0(1),mass0(2),lumin(1),lumin(2), & teff1,teff2,radc(1),radc(2), & menv(1),menv(2),renv(1),renv(2), @@ -1392,7 +1398,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin(1),bhspin(2), & deltam1_bcm,deltam2_bcm,formation(1), - & formation(2),binstate,mergertype,'bpp') + & formation(2),binstate,mergertype,z,'bpp') CALL kick(kw,mass(k),mt,mass(3-k),ecc,sep,jorb,vk,k, & rad(3-k),fallback,sigmahold,kick_info,disrupt,bkick) @@ -1600,7 +1606,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & mass(1),mass(2),kstar(1),kstar(2),sep, & tb,ecc,rrl1,rrl2, & aj(1),aj(2),tms(1),tms(2), - & massc(1),massc(2),rad(1),rad(2), + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad(1),rad(2), & mass0(1),mass0(2),lumin(1),lumin(2), & teff1,teff2,radc(1),radc(2), & menv(1),menv(2),renv(1),renv(2), @@ -1608,7 +1615,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin(1),bhspin(2), & deltam1_bcm,deltam2_bcm,formation(1), - & formation(2),binstate,mergertype,'bpp') + & formation(2),binstate,mergertype,z,'bpp') if(snova)then bpp(jp,11) = 2.0 dtm = 0.d0 @@ -1664,7 +1671,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & mass(1),mass(2),kstar(1),kstar(2), & sep,tb,ecc,rrl1,rrl2, & aj(1),aj(2),tms(1),tms(2), - & massc(1),massc(2),rad(1),rad(2), + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad(1),rad(2), & mass0(1),mass0(2),lumin(1),lumin(2), & teff1,teff2,radc(1),radc(2), & menv(1),menv(2),renv(1),renv(2), @@ -1672,7 +1680,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin(1),bhspin(2), & deltam1_bcm,deltam2_bcm,formation(1), - & formation(2),binstate,mergertype,'bcm') + & formation(2),binstate,mergertype,z,'bcm') if(isave) tsave = tsave + dtp if(output) write(*,*)'bcm1',kstar(1),kstar(2),mass(1), & mass(2),rad(1),rad(2),ospin(1),ospin(2),jspin(1) @@ -1812,7 +1820,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & mass(1),mass(2),kstar(1),kstar(2),sep, & tb,ecc,rrl1,rrl2, & aj(1),aj(2),tms(1),tms(2), - & massc(1),massc(2),rad(1),rad(2), + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad(1),rad(2), & mass0(1),mass0(2),lumin(1),lumin(2), & teff1,teff2,radc(1),radc(2), & menv(1),menv(2),renv(1),renv(2), @@ -1820,7 +1829,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin(1),bhspin(2), & deltam1_bcm,deltam2_bcm,formation(1), - & formation(2),binstate,mergertype,'bpp') + & formation(2),binstate,mergertype,z,'bpp') endif * iter = iter + 1 @@ -1904,7 +1913,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & mass(1),mass(2),kstar(1),kstar(2),sep, & tb,ecc,rrl1,rrl2, & aj(1),aj(2),tms(1),tms(2), - & massc(1),massc(2),rad(1),rad(2), + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad(1),rad(2), & mass0(1),mass0(2),lumin(1),lumin(2), & teff1,teff2,radc(1),radc(2), & menv(1),menv(2),renv(1),renv(2), @@ -1912,7 +1922,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin(1),bhspin(2), & deltam1_bcm,deltam2_bcm,formation(1), - & formation(2),binstate,mergertype,'bpp') + & formation(2),binstate,mergertype,z,'bpp') * if(check_dtp.eq.1)then CALL checkstate(dtp,dtp_original,tsave,tphys,tphysf, @@ -1960,7 +1970,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & mass(1),mass(2),kstar(1),kstar(2), & sep,tb,ecc,rrl1,rrl2, & aj(1),aj(2),tms(1),tms(2), - & massc(1),massc(2),rad(1),rad(2), + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad(1),rad(2), & mass0(1),mass0(2),lumin(1),lumin(2), & teff1,teff2,radc(1),radc(2), & menv(1),menv(2),renv(1),renv(2), @@ -1968,7 +1979,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin(1),bhspin(2), & deltam1_bcm,deltam2_bcm,formation(1), - & formation(2),binstate,mergertype,'bcm') + & formation(2),binstate,mergertype,z,'bcm') if(output) write(*,*)'bcm2:',kstar(1),kstar(2),mass(1), & mass(2),rad(1),rad(2),ospin(1),ospin(2),jspin(1) * & mass(2),rad(1),rad(2),ospin(1),ospin(2),b01_bcm,b02_bcm,jspin(1) @@ -2401,7 +2412,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & kstar(1),kstar(2),sep, & tb,ecc,rrl1,rrl2, & aj(1),aj(2),tms(1),tms(2), - & massc(1),massc(2),rad(1),rad(2), + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad(1),rad(2), & mass0(1),mass0(2),lumin(1),lumin(2), & teff1,teff2,radc(1),radc(2), & menv(1),menv(2),renv(1),renv(2), @@ -2409,7 +2421,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin(1),bhspin(2), & deltam1_bcm,deltam2_bcm,formation(1), - & formation(2),binstate,mergertype,'bpp') + & formation(2),binstate,mergertype,z,'bpp') CALL comenv(mass0(j1),mass(j1),massc(j1),aj(j1),jspin(j1), & kstar(j1),mass0(j2),mass(j2),massc(j2),aj(j2), @@ -2524,7 +2536,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & kstar(1),kstar(2),sep, & tb,ecc,rrl1,rrl2, & aj(1),aj(2),tms(1),tms(2), - & massc(1),massc(2),rad(1),rad(2), + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad(1),rad(2), & mass0(1),mass0(2),lumin(1),lumin(2), & teff1,teff2,radc(1),radc(2), & menv(1),menv(2),renv(1),renv(2), @@ -2532,7 +2545,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin(1),bhspin(2), & deltam1_bcm,deltam2_bcm,formation(1), - & formation(2),binstate,mergertype,'bpp') + & formation(2),binstate,mergertype,z,'bpp') * epoch(j1) = tphys - aj(j1) com = .false. @@ -3650,7 +3663,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & mass(1),mass(2),kstar(1),kstar(2), & sep,tb,ecc,rrl1,rrl2, & aj(1),aj(2),tms(1),tms(2), - & massc(1),massc(2),rad(1),rad(2), + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad(1),rad(2), & mass0(1),mass0(2),lumin(1),lumin(2), & teff1,teff2,radc(1),radc(2), & menv(1),menv(2),renv(1),renv(2), @@ -3658,7 +3672,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin(1),bhspin(2), & deltam1_bcm,deltam2_bcm,formation(1), - & formation(2),binstate,mergertype,'bpp') + & formation(2),binstate,mergertype,z,'bpp') CALL kick(kw,mass(k),mt,mass(3-k),ecc,sep,jorb,vk,k, & rad(3-k),fallback,sigmahold,kick_info,disrupt,bkick) sigma = sigmahold !reset sigma after possible ECSN kick dist. Remove this if u want some kick link to the intial pulsar values... @@ -3791,7 +3805,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & mass(1),mass(2),kstar(1),kstar(2), & sep,tb,ecc,rrl1,rrl2, & aj(1),aj(2),tms(1),tms(2), - & massc(1),massc(2),rad(1),rad(2), + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad(1),rad(2), & mass0(1),mass0(2),lumin(1),lumin(2), & teff1,teff2,radc(1),radc(2), & menv(1),menv(2),renv(1),renv(2), @@ -3799,7 +3814,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin(1),bhspin(2), & deltam1_bcm,deltam2_bcm,formation(1), - & formation(2),binstate,mergertype,'bcm') + & formation(2),binstate,mergertype,z,'bcm') if(isave) tsave = tsave + dtp if(output) write(*,*)'bcm3:',kstar(1),kstar(2),mass(1), & mass(2),rad(1),rad(2),ospin(1),ospin(2),jspin(1) @@ -3836,7 +3851,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & mass(1),mass(2),kstar(1),kstar(2), & sep,tb,ecc,rrl1,rrl2, & aj(1),aj(2),tms(1),tms(2), - & massc(1),massc(2),rad(1),rad(2), + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad(1),rad(2), & mass0(1),mass0(2),lumin(1),lumin(2), & teff1,teff2,radc(1),radc(2), & menv(1),menv(2),renv(1),renv(2), @@ -3844,7 +3860,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin(1),bhspin(2), & deltam1_bcm,deltam2_bcm,formation(1), - & formation(2),binstate,mergertype,'bpp') + & formation(2),binstate,mergertype,z,'bpp') endif * @@ -3884,7 +3900,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & mass(1),mass(2),kstar(1),kstar(2),sep, & tb,ecc,rrl1,rrl2, & aj(1),aj(2),tms(1),tms(2), - & massc(1),massc(2),rad(1),rad(2), + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad(1),rad(2), & mass0(1),mass0(2),lumin(1),lumin(2), & teff1,teff2,radc(1),radc(2), & menv(1),menv(2),renv(1),renv(2), @@ -3892,7 +3909,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin(1),bhspin(2), & deltam1_bcm,deltam2_bcm,formation(1), - & formation(2),binstate,mergertype,'bpp') + & formation(2),binstate,mergertype,z,'bpp') dtm = 0.d0 goto 4 endif @@ -3935,7 +3952,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & mass(1),mass(2),kstar(1),kstar(2),sep, & tb,ecc,rrl1,rrl2, & aj(1),aj(2),tms(1),tms(2), - & massc(1),massc(2),rad(1),rad(2), + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad(1),rad(2), & mass0(1),mass0(2),lumin(1),lumin(2), & teff1,teff2,radc(1),radc(2), & menv(1),menv(2),renv(1),renv(2), @@ -3943,7 +3961,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin(1),bhspin(2), & deltam1_bcm,deltam2_bcm,formation(1), - & formation(2),binstate,mergertype,'bpp') + & formation(2),binstate,mergertype,z,'bpp') * kcomp1 = kstar(j1) kcomp2 = kstar(j2) @@ -3982,7 +4000,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & kstar(1),kstar(2),sep, & tb,ecc,rrl1,rrl2, & aj(1),aj(2),tms(1),tms(2), - & massc(1),massc(2),rad(1),rad(2), + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad(1),rad(2), & mass0(1),mass0(2),lumin(1),lumin(2), & teff1,teff2,radc(1),radc(2), & menv(1),menv(2),renv(1),renv(2), @@ -3990,7 +4009,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin(1),bhspin(2), & deltam1_bcm,deltam2_bcm,formation(1), - & formation(2),binstate,mergertype,'bpp') + & formation(2),binstate,mergertype,z,'bpp') CALL comenv(mass0(j1),mass(j1),massc(j1),aj(j1),jspin(j1), & kstar(j1),mass0(j2),mass(j2),massc(j2),aj(j2), & jspin(j2),kstar(j2),zpars,ecc,sep,jorb,coel,j1,j2, @@ -4066,7 +4085,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & kstar(1),kstar(2),sep, & tb,ecc,rrl1,rrl2, & aj(1),aj(2),tms(1),tms(2), - & massc(1),massc(2),rad(1),rad(2), + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad(1),rad(2), & mass0(1),mass0(2),lumin(1),lumin(2), & teff1,teff2,radc(1),radc(2), & menv(1),menv(2),renv(1),renv(2), @@ -4074,7 +4094,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin(1),bhspin(2), & deltam1_bcm,deltam2_bcm,formation(1), - & formation(2),binstate,mergertype,'bpp') + & formation(2),binstate,mergertype,z,'bpp') CALL comenv(mass0(j2),mass(j2),massc(j2),aj(j2),jspin(j2), & kstar(j2),mass0(j1),mass(j1),massc(j1),aj(j1), & jspin(j1),kstar(j1),zpars,ecc,sep,jorb,coel,j2,j1, @@ -4155,7 +4175,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & kstar(1),kstar(2),sep, & tb,ecc,rrl1,rrl2, & aj(1),aj(2),tms(1),tms(2), - & massc(1),massc(2),rad(1),rad(2), + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad(1),rad(2), & mass0(1),mass0(2),lumin(1),lumin(2), & teff1,teff2,radc(1),radc(2), & menv(1),menv(2),renv(1),renv(2), @@ -4163,7 +4184,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin(1),bhspin(2), & deltam1_bcm,deltam2_bcm,formation(1), - & formation(2),binstate,mergertype,'bpp') + & formation(2),binstate,mergertype,z,'bpp') endif epoch(1) = tphys - aj(1) epoch(2) = tphys - aj(2) @@ -4212,7 +4233,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & mass(1),mass(2),kstar(1),kstar(2),sep, & tb,ecc,rrl1,rrl2, & aj(1),aj(2),tms(1),tms(2), - & massc(1),massc(2),rad(1),rad(2), + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad(1),rad(2), & mass0(1),mass0(2),lumin(1),lumin(2), & teff1,teff2,radc(1),radc(2), & menv(1),menv(2),renv(1),renv(2), @@ -4220,7 +4242,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin(1),bhspin(2), & deltam1_bcm,deltam2_bcm,formation(1), - & formation(2),binstate,mergertype,'bpp') + & formation(2),binstate,mergertype,z,'bpp') dtm = 0.d0 * * Reset orbital parameters as separation may have changed. @@ -4285,7 +4307,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & kstar(1),kstar(2),0.d0, & 0.d0,-1.d0,0.d0,ngtv, & aj(1),aj(2),tms(1),tms(2), - & massc(1),massc(2),rad(1),rad(2), + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad(1),rad(2), & mass0(1),mass0(2),lumin(1),lumin(2), & teff1,teff2,radc(1),radc(2), & menv(1),menv(2),renv(1),renv(2), @@ -4293,7 +4316,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin(1),bhspin(2), & deltam1_bcm,deltam2_bcm,formation(1), - & formation(2),binstate,mergertype,'bpp') + & formation(2),binstate,mergertype,z,'bpp') elseif(ecc.gt.1.d0)then * * Binary dissolved by a supernova or tides. @@ -4327,7 +4350,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & kstar(1),kstar(2),sep, & tb,ecc,0.d0,ngtv2, & aj(1),aj(2),tms(1),tms(2), - & massc(1),massc(2),rad(1),rad(2), + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad(1),rad(2), & mass0(1),mass0(2),lumin(1),lumin(2), & teff1,teff2,radc(1),radc(2), & menv(1),menv(2),renv(1),renv(2), @@ -4335,7 +4359,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin(1),bhspin(2), & deltam1_bcm,deltam2_bcm,formation(1), - & formation(2),binstate,mergertype,'bpp') + & formation(2),binstate,mergertype,z,'bpp') else evolve_type = 9.0 teff1 = 1000.d0*((1130.d0*lumin(1)/ @@ -4361,7 +4385,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & kstar(1),kstar(2),0.d0, & 0.d0,0.d0,0.d0,ngtv, & aj(1),aj(2),tms(1),tms(2), - & massc(1),massc(2),rad(1),rad(2), + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad(1),rad(2), & mass0(1),mass0(2),lumin(1),lumin(2), & teff1,teff2,radc(1),radc(2), & menv(1),menv(2),renv(1),renv(2), @@ -4369,7 +4394,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin(1),bhspin(2), & deltam1_bcm,deltam2_bcm,formation(1), - & formation(2),binstate,mergertype,'bpp') + & formation(2),binstate,mergertype,z,'bpp') endif endif if(kstar(2).eq.15)then @@ -4443,7 +4468,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & mass(1),mass(2),kstar(1),kstar(2),sep, & tb,ecc,rrl1,rrl2, & aj(1),aj(2),tms(1),tms(2), - & massc(1),massc(2),rad(1),rad(2), + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad(1),rad(2), & mass0(1),mass0(2),lumin(1),lumin(2), & teff1,teff2,radc(1),radc(2), & menv(1),menv(2),renv(1),renv(2), @@ -4451,7 +4477,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin(1),bhspin(2), & deltam1_bcm,deltam2_bcm,formation(1), - & formation(2),binstate,mergertype,'bpp') + & formation(2),binstate,mergertype,z,'bpp') endif * if(kstar(1).eq.15.and.bpp(jp,4).lt.15.0)then @@ -4487,7 +4513,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & kstar(1),kstar(2),0.d0, & 0.d0,-1.d0,0.d0,ngtv, & aj(1),aj(2),tms(1),tms(2), - & massc(1),massc(2),rad(1),rad(2), + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad(1),rad(2), & mass0(1),mass0(2),lumin(1),lumin(2), & teff1,teff2,radc(1),radc(2), & menv(1),menv(2),renv(1),renv(2), @@ -4495,7 +4522,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin(1),bhspin(2), & deltam1_bcm,deltam2_bcm,formation(1), - & formation(2),binstate,mergertype,'bpp') + & formation(2),binstate,mergertype,z,'bpp') elseif(kstar(1).eq.15.and.kstar(2).eq.15)then * * Cases of accretion induced supernova or single star supernova. @@ -4525,7 +4552,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & kstar(1),kstar(2),0.d0, & 0.d0,0.d0,0.d0,ngtv2, & aj(1),aj(2),tms(1),tms(2), - & massc(1),massc(2),rad(1),rad(2), + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad(1),rad(2), & mass0(1),mass0(2),lumin(1),lumin(2), & teff1,teff2,radc(1),radc(2), & menv(1),menv(2),renv(1),renv(2), @@ -4533,7 +4561,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin(1),bhspin(2), & deltam1_bcm,deltam2_bcm,formation(1), - & formation(2),binstate,mergertype,'bpp') + & formation(2),binstate,mergertype,z,'bpp') else evolve_type = 10.0 !added by PA for systems that stop evolving halfway @@ -4563,7 +4591,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & kstar(1),kstar(2),sep, & tb,ecc,rrl1,rrl2, & aj(1),aj(2),tms(1),tms(2), - & massc(1),massc(2),rad(1),rad(2), + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad(1),rad(2), & mass0(1),mass0(2),lumin(1),lumin(2), & teff1,teff2,radc(1),radc(2), & menv(1),menv(2),renv(1),renv(2), @@ -4571,7 +4600,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin(1),bhspin(2), & deltam1_bcm,deltam2_bcm,formation(1), - & formation(2),binstate,mergertype,'bpp') + & formation(2),binstate,mergertype,z,'bpp') endif endif * @@ -4624,7 +4653,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & mass(1),mass(2),kstar(1),kstar(2), & sep,tb,ecc,rrl1,rrl2, & aj(1),aj(2),tms(1),tms(2), - & massc(1),massc(2),rad(1),rad(2), + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad(1),rad(2), & mass0(1),mass0(2),lumin(1),lumin(2), & teff1,teff2,radc(1),radc(2), & menv(1),menv(2),renv(1),renv(2), @@ -4632,7 +4662,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), & epoch(2),bhspin(1),bhspin(2), & deltam1_bcm,deltam2_bcm,formation(1), - & formation(2),binstate,mergertype,'bcm') + & formation(2),binstate,mergertype,z,'bcm') if(output) write(*,*)'bcm4:',kstar(1),kstar(2),mass(1), & mass(2),rad(1),rad(2),ospin(1),ospin(2),jspin(1), & tphys,tphysf diff --git a/src/cosmic/src/hrdiag.f b/src/cosmic/src/hrdiag.f index dc19af197..080a9bf41 100644 --- a/src/cosmic/src/hrdiag.f +++ b/src/cosmic/src/hrdiag.f @@ -27,7 +27,7 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, * real*8 mass,aj,mt,tm,tn,tscls(20),lums(10),GB(10),zpars(20),met real*8 bhspin - real*8 r,lum,mc_he,mc_co,rc,menv,renv,k2 + real*8 r,lum,rc,menv,renv,k2,mc real*8 mch,mlp,tiny * parameter(mch=1.44d0,mlp=12.d0,tiny=1.0d-14) parameter(mlp=12.d0,tiny=1.0d-14) @@ -178,8 +178,14 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, endif eta = mctmsf(mass) tau = (aj - tm)/thg +* WRITE(*,*)'hrdiag: 184, k=',kidx,' mass=',mass +* WRITE(*,*)'mc=',mc,' mcx=',mcx + mc = ((1.d0 - tau)*eta + tau)*mc mc = MAX(mc,mcx) + + mc_he(kidx) = mc + mc_co(kidx) = 0.0 * * Test whether core mass has reached total mass. * @@ -190,6 +196,8 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, * Zero-age helium star * mc = 0.d0 + mc_he(kidx) = 0.d0 + mc_co(kidx) = 0.d0 mass = mt kw = 7 CALL star(kw,mass,mt,tm,tn,tscls,lums,GB,zpars) @@ -198,6 +206,8 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, * Zero-age helium white dwarf. * mc = mt + mc_he(kidx) = 0.d0 + mc_co(kidx) = 0.d0 mass = mt kw = 10 endif @@ -235,6 +245,8 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, if(mass.le.zpars(2))then * Star has a degenerate He core which grows on the GB mc = mcgbf(lum,GB,lums(6)) + mc_he(kidx) = mc + mc_co(kidx) = 0.0 else * Star has a non-degenerate He core which may grow, but * only slightly, on the GB @@ -242,6 +254,11 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, mcx = mcheif(mass,zpars(2),zpars(9)) mcy = mcheif(mass,zpars(2),zpars(10)) mc = mcx + (mcy - mcx)*tau +* WRITE(*,*)'hrdiag249: k=',kidx,'mass=',mass,'kw=',kw +* WRITE(*,*)'mc=',mc,' mcx=',mcx,' mcy-mcx=',mcy-mcx +* WRITE(*,*) + mc_he(kidx) = mc + mc_co(kidx) = 0.0 endif r = rgbf(mt,lum) rg = r @@ -252,6 +269,8 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, * Zero-age helium star * mc = 0.d0 + mc_he(kidx) = mc + mc_co(kidx) = 0.0 mass = mt kw = 7 CALL star(kw,mass,mt,tm,tn,tscls,lums,GB,zpars) @@ -260,6 +279,8 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, * Zero-age helium white dwarf. * mc = mt + mc_he(kidx) = mc + mc_co(kidx) = 0.0 mass = mt kw = 10 endif @@ -280,7 +301,12 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, mcx = mcheif(mass,zpars(2),zpars(10)) endif tau = (aj - tscls(2))/tscls(3) +* here, mcx is the He core mass, and the other part is the C core mc = mcx + (mcagbf(mass) - mcx)*tau + WRITE(*,*)'hrdiag: mc=',mc,' mcx=',mcx, 'kw=',kw,' k=',kidx + WRITE(*,*)'hrdiag: mc_co=',(mcagbf(mass) - mcx)*tau + mc_he(kidx) = mcx + (mcagbf(mass) - mcx)*tau + mc_co(kidx) = 0.0 * if(mass.le.zpars(2))then lx = lums(5) @@ -401,6 +427,10 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, if(aj.lt.tscls(13))then mcx = mcgbtf(aj,GB(8),GB,tscls(7),tscls(8),tscls(9)) mc = mcbagb + WRITE(*,*)'hrdiag 430: mc=',mc,' mcx=',mcx,' kw=',kw + WRITE(*,*)'hrdiag: mcbagb=',mcbagb + mc_co(kidx) = mcx + mc_he(kidx) = mcbagb lum = lmcgbf(mcx,GB) if(mt.le.mc)then * @@ -412,6 +442,8 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, mt = mc mass = mt mc = mcx + mc_co(kidx) = mc + mc_he(kidx) = 0.d0 CALL star(kw,mass,mt,tm,tn,tscls,lums,GB,zpars) if(mc.le.GB(7))then aj = tscls(4) - (1.d0/((GB(5)-1.d0)*GB(8)*GB(4)))* @@ -438,6 +470,8 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, mcy = mc mc = mc - lambdahrdiag*(mcy-mcx) mcx = mc + mc_co(kidx) = mcx + mc_he(kidx) = mcy mcmax = MIN(mt,mcmax) endif r = ragbf(mt,lum,zpars(2)) @@ -477,6 +511,8 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, * which is why we subject mass and mt to mass loss for * this phase. mc = 0.d0 + mc_he(kidx) = mc + mc_co(kidx) = 0.0 if(mt.lt.zpars(10)) kw = 10 else * @@ -491,6 +527,8 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, r = rg endif mc = mcgbf(lum,GB,lums(6)) + mc_he(kidx) = 0.d0 + mc_co(kidx) = mc mtc = MIN(mt,1.45d0*mt-0.31d0) mcmax = MIN(mtc,MAX(mch,0.773d0*mass-0.35d0)) if(mcmax-mc.lt.tiny)then From 988e9eb2f93aae443ea70d31d85d7f1519f3a913 Mon Sep 17 00:00:00 2001 From: Katie Breivik Date: Mon, 11 Aug 2025 12:00:47 -0400 Subject: [PATCH 10/19] fixing pointers in docs that didnt carry over right... sorry for adding here.. --- docs/pages/config/config_insert_bse.html | 4 ++-- docs/pages/inifile.rst | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/pages/config/config_insert_bse.html b/docs/pages/config/config_insert_bse.html index f47ce4f5a..25da09dc1 100644 --- a/docs/pages/config/config_insert_bse.html +++ b/docs/pages/config/config_insert_bse.html @@ -404,12 +404,12 @@

kickflag

Sets the particular natal kick prescription to use. Note that sigmadiv, bhflag, bhsigmafrac, aic, and ussn, which are described below, are only used when abs(kickflag)=1. Positive values use the Pfahl+2002 prescription for handling natal kicks.

Default: 1

-
+
Option details

-
  • 1: The standard COSMIC kick prescription, where kicks are drawn from a bimodal distribution with standard FeCCSN getting a kick drawn from a Maxwellian distribution with dispersion parameter sigma and ECSN/USSN are drawn according to sigmadiv. This setting has additional possible options for bhflag, bhsigmafrac, aic and ussn.
  • 2: Natal kicks are drawn according to sigma and scaled by the ejecta mass and remnant mass following Eq. 1 of Giacobbo & Mapelli 2020 with their default parameters (\(m_{\rm NS = 1.2 {\rm M_\odot}\), \(m_{\rm ej = 9 {\rm M_\odot}\))
  • 3: Natal kicks are drawn according to sigma and scaled by just the ejecta mass following Eq. 2 of Giacobbo & Mapelli 2020, which does not scale the kick by (\(m_{\rm NS\)
  • 4: Natal kicks are drawn according to Eq. 1 of Bray & Eldridge 2016, with their default parameters (\(\alpha=70 \, {\rm km/s}, \beta = 120 \, {\rm km/s)}
  • negative values: Same as above settings but using the old Kiel & Hurley 2009 prescription for changing the orbital configuration of the binary, available for reproducibility purposes but not recommended for new work
+
  • 1: The standard COSMIC kick prescription, where kicks are drawn from a bimodal distribution with standard FeCCSN getting a kick drawn from a Maxwellian distribution with dispersion parameter sigma and ECSN/USSN are drawn according to sigmadiv. This setting has additional possible options for bhflag, bhsigmafrac, aic and ussn.
  • 2: Natal kicks are drawn according to sigma and scaled by the ejecta mass and remnant mass following Eq. 1 of Giacobbo & Mapelli 2020 with their default parameters (\(m_{\rm NS = 1.2 {\rm M_\odot}\), \(m_{\rm ej = 9 {\rm M_\odot}\))
  • 3: Natal kicks are drawn according to sigma and scaled by just the ejecta mass following Eq. 2 of Giacobbo & Mapelli 2020, which does not scale the kick by (\(m_{\rm NS\)
  • 4: Natal kicks are drawn according to Eq. 1 of Bray & Eldridge 2016, with their default parameters (\(\alpha=70 \, {\rm km/s}, \beta = 120 \, {\rm km/s)}
  • 5: Follows the same prescription as 1, but uses the kick prescription described in Disberg & Mandel 2025 for CCSN.
  • negative values: Same as above settings but using the old Kiel & Hurley 2009 prescription for changing the orbital configuration of the binary, available for reproducibility purposes but not recommended for new work
diff --git a/docs/pages/inifile.rst b/docs/pages/inifile.rst index 602486066..b3f1916a1 100644 --- a/docs/pages/inifile.rst +++ b/docs/pages/inifile.rst @@ -17,7 +17,7 @@ The buttons below link to the most recent stable and unstable default inifiles f .. raw:: html
-
Latest stable INIFILE
+
Latest stable INIFILE
Latest development INIFILE
From 9c7e81d8a6c0e7d1fbcf1cbd02d6c51a83b8a149 Mon Sep 17 00:00:00 2001 From: Katie Breivik Date: Mon, 8 Dec 2025 16:14:19 -0800 Subject: [PATCH 11/19] some more debugging on helium/carbon, still working on AGB phases --- src/cosmic/src/evolv2.f | 2 +- src/cosmic/src/hrdiag.f | 22 +++++++++++++--------- 2 files changed, 14 insertions(+), 10 deletions(-) diff --git a/src/cosmic/src/evolv2.f b/src/cosmic/src/evolv2.f index c48e19e37..0ab588748 100644 --- a/src/cosmic/src/evolv2.f +++ b/src/cosmic/src/evolv2.f @@ -156,7 +156,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, INTEGER loop,iter,intpol,k,ip,jp,j1,j2 INTEGER bcm_index_out, bpp_index_out INTEGER kcomp1,kcomp2,formation(2) - PARAMETER(loop=20000) + PARAMETER(loop=200000) INTEGER kstar(2),kw,kst,kw1,kw2,kmin,kmax INTEGER kstar1_bpp,kstar2_bpp * diff --git a/src/cosmic/src/hrdiag.f b/src/cosmic/src/hrdiag.f index 080a9bf41..ca4e83072 100644 --- a/src/cosmic/src/hrdiag.f +++ b/src/cosmic/src/hrdiag.f @@ -301,11 +301,11 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, mcx = mcheif(mass,zpars(2),zpars(10)) endif tau = (aj - tscls(2))/tscls(3) -* here, mcx is the He core mass, and the other part is the C core +* here, mcx is the helium core mass at helium ignition mc = mcx + (mcagbf(mass) - mcx)*tau WRITE(*,*)'hrdiag: mc=',mc,' mcx=',mcx, 'kw=',kw,' k=',kidx - WRITE(*,*)'hrdiag: mc_co=',(mcagbf(mass) - mcx)*tau - mc_he(kidx) = mcx + (mcagbf(mass) - mcx)*tau + WRITE(*,*)'hrdiag: mc_he=',(mcagbf(mass) - mcx)*tau + mc_he(kidx) = mc mc_co(kidx) = 0.0 * if(mass.le.zpars(2))then @@ -466,12 +466,12 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, * lambdahrdiag = MIN(0.9d0,0.3d0+0.001d0*mass**5) tau = tscls(13) - mcx = mcgbtf(tau,GB(2),GB,tscls(10),tscls(11),tscls(12)) - mcy = mc - mc = mc - lambdahrdiag*(mcy-mcx) + mcy = mcgbtf(tau,GB(2),GB,tscls(10),tscls(11),tscls(12)) + mcx = mc + mc = mcy - lambdahrdiag*(mcx-mcy) mcx = mc mc_co(kidx) = mcx - mc_he(kidx) = mcy + mc_he(kidx) = mcy - mcx mcmax = MIN(mt,mcmax) endif r = ragbf(mt,lum,zpars(2)) @@ -510,9 +510,10 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, * Star has no core mass and hence no memory of its past * which is why we subject mass and mt to mass loss for * this phase. +*KB: no helium core mass for stripped stars; no CO core yet since He MS mc = 0.d0 - mc_he(kidx) = mc - mc_co(kidx) = 0.0 + mc_he(kidx) = 0.d0 + mc_co(kidx) = 0.d0 if(mt.lt.zpars(10)) kw = 10 else * @@ -527,6 +528,9 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, r = rg endif mc = mcgbf(lum,GB,lums(6)) +* +*KB: helium core mass is always 0 for stripped stars; now calculate CO core mass +* mc_he(kidx) = 0.d0 mc_co(kidx) = mc mtc = MIN(mt,1.45d0*mt-0.31d0) From e587b50a59022bcd8191750a69d3ea72e0b5852e Mon Sep 17 00:00:00 2001 From: Katie Breivik Date: Mon, 15 Dec 2025 14:37:13 -0800 Subject: [PATCH 12/19] some prints and one switch to he vs co --- src/cosmic/src/assign_remnant.f | 10 +++++----- src/cosmic/src/evolv2.f | 3 +++ src/cosmic/src/hrdiag.f | 2 +- 3 files changed, 9 insertions(+), 6 deletions(-) diff --git a/src/cosmic/src/assign_remnant.f b/src/cosmic/src/assign_remnant.f index 68a3464d9..4842bdecf 100644 --- a/src/cosmic/src/assign_remnant.f +++ b/src/cosmic/src/assign_remnant.f @@ -4,8 +4,8 @@ SUBROUTINE assign_remnant(zpars,mc,mcbagb,mass,mt,kw,bhspin,kidx) INCLUDE 'const_bse.h' common /fall/fallback - REAL*8 mass_preSN, mHe_preSN, massc_preSN - COMMON mass_preSN, mHe_preSN, massc_preSN +* REAL*8 mass_preSN, mHe_preSN, massc_preSN +* COMMON mass_preSN, mHe_preSN, massc_preSN REAL*8 fallback REAL ran3 EXTERNAL ran3 @@ -73,9 +73,9 @@ SUBROUTINE assign_remnant(zpars,mc,mcbagb,mass,mt,kw,bhspin,kidx) * else * Store values in common block - mass_preSN = mt - mHe_preSN = mcbagb - mc - massc_preSN = mc +* mass_preSN = mt +* mHe_preSN = mc_he(kidx) +* massc_preSN = mc_co(kidx) if(ecsn.gt.0.d0.and.mcbagb.lt.ecsn_mlow)then * * Star is not massive enough to ignite C burning. diff --git a/src/cosmic/src/evolv2.f b/src/cosmic/src/evolv2.f index 0ab588748..402495548 100644 --- a/src/cosmic/src/evolv2.f +++ b/src/cosmic/src/evolv2.f @@ -1580,6 +1580,9 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, * if((tphys.lt.tiny.and.ABS(dtm).lt.tiny.and. & (mass2i.lt.0.1d0.or..not.sgl)).or.snova)then + if(kstar(1).eq.14)then + WRITE(*,*)'BH???' + endif evolve_type = 1.d0 rrl1 = rad(1)/rol(1) rrl2 = rad(2)/rol(2) diff --git a/src/cosmic/src/hrdiag.f b/src/cosmic/src/hrdiag.f index ca4e83072..f1988934d 100644 --- a/src/cosmic/src/hrdiag.f +++ b/src/cosmic/src/hrdiag.f @@ -430,7 +430,7 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, WRITE(*,*)'hrdiag 430: mc=',mc,' mcx=',mcx,' kw=',kw WRITE(*,*)'hrdiag: mcbagb=',mcbagb mc_co(kidx) = mcx - mc_he(kidx) = mcbagb + mc_he(kidx) = mcbagb - mcx lum = lmcgbf(mcx,GB) if(mt.le.mc)then * From b08ef98f3ccf812752c814aeafe60c0a0ac4addf Mon Sep 17 00:00:00 2001 From: Katie Breivik Date: Mon, 15 Dec 2025 17:05:52 -0800 Subject: [PATCH 13/19] more debugging --- src/cosmic/src/evolv2.f | 9 +++++---- src/cosmic/src/hrdiag.f | 10 ++++++---- 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/src/cosmic/src/evolv2.f b/src/cosmic/src/evolv2.f index 402495548..095ce35df 100644 --- a/src/cosmic/src/evolv2.f +++ b/src/cosmic/src/evolv2.f @@ -156,7 +156,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, INTEGER loop,iter,intpol,k,ip,jp,j1,j2 INTEGER bcm_index_out, bpp_index_out INTEGER kcomp1,kcomp2,formation(2) - PARAMETER(loop=200000) + PARAMETER(loop=20000) INTEGER kstar(2),kw,kst,kw1,kw2,kmin,kmax INTEGER kstar1_bpp,kstar2_bpp * @@ -1580,10 +1580,11 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, * if((tphys.lt.tiny.and.ABS(dtm).lt.tiny.and. & (mass2i.lt.0.1d0.or..not.sgl)).or.snova)then - if(kstar(1).eq.14)then - WRITE(*,*)'BH???' - endif evolve_type = 1.d0 + if(snova)then +* We should capture to evol_type change for SN as an evolutionary change + evolve_type = 2.d0 + endif rrl1 = rad(1)/rol(1) rrl2 = rad(2)/rol(2) teff1 = 1000.d0*((1130.d0*lumin(1)/ diff --git a/src/cosmic/src/hrdiag.f b/src/cosmic/src/hrdiag.f index f1988934d..af56901d1 100644 --- a/src/cosmic/src/hrdiag.f +++ b/src/cosmic/src/hrdiag.f @@ -109,6 +109,8 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, * Main sequence star. * mc = 0.d0 + mc_he(kidx) = 0.d0 + mc_co(kidx) = 0.d0 tau = aj/tm thook = thookf(mass)*tscls(1) zeta = 0.01d0 @@ -303,8 +305,8 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, tau = (aj - tscls(2))/tscls(3) * here, mcx is the helium core mass at helium ignition mc = mcx + (mcagbf(mass) - mcx)*tau - WRITE(*,*)'hrdiag: mc=',mc,' mcx=',mcx, 'kw=',kw,' k=',kidx - WRITE(*,*)'hrdiag: mc_he=',(mcagbf(mass) - mcx)*tau +* WRITE(*,*)'hrdiag: mc=',mc,' mcx=',mcx, 'kw=',kw,' k=',kidx +* WRITE(*,*)'hrdiag: mc_he=',(mcagbf(mass) - mcx)*tau mc_he(kidx) = mc mc_co(kidx) = 0.0 * @@ -427,8 +429,8 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, if(aj.lt.tscls(13))then mcx = mcgbtf(aj,GB(8),GB,tscls(7),tscls(8),tscls(9)) mc = mcbagb - WRITE(*,*)'hrdiag 430: mc=',mc,' mcx=',mcx,' kw=',kw - WRITE(*,*)'hrdiag: mcbagb=',mcbagb +* WRITE(*,*)'hrdiag 430: mc=',mc,' mcx=',mcx,' kw=',kw +* WRITE(*,*)'hrdiag: mcbagb=',mcbagb mc_co(kidx) = mcx mc_he(kidx) = mcbagb - mcx lum = lmcgbf(mcx,GB) From f3a754777365fbd4679dac3b1cdf76ce5c4a4537 Mon Sep 17 00:00:00 2001 From: Katie Breivik Date: Tue, 16 Dec 2025 09:06:00 -0800 Subject: [PATCH 14/19] dont make metallicity specified by default --- src/cosmic/evolve.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/cosmic/evolve.py b/src/cosmic/evolve.py index 6c9d105e9..2d9140ffe 100644 --- a/src/cosmic/evolve.py +++ b/src/cosmic/evolve.py @@ -60,7 +60,7 @@ INTEGER_COLUMNS = ["bin_state", "bin_num", "kstar_1", "kstar_2", "SN_1", "SN_2", "evol_type"] -BPP_COLUMNS = ['tphys', 'metallicity', 'mass_1', 'mass_2', 'kstar_1', 'kstar_2', +BPP_COLUMNS = ['tphys', 'mass_1', 'mass_2', 'kstar_1', 'kstar_2', 'sep', 'porb', 'ecc', 'RRLO_1', 'RRLO_2', 'evol_type', 'aj_1', 'aj_2', 'tms_1', 'tms_2', 'massc_he_1', 'massc_he_2', 'massc_co_1', 'massc_co_2', 'rad_1', 'rad_2', @@ -70,7 +70,7 @@ 'tacc_1', 'tacc_2', 'epoch_1', 'epoch_2', 'bhspin_1', 'bhspin_2'] -BCM_COLUMNS = ['tphys', 'metallicity', 'kstar_1', 'mass0_1', 'mass_1', 'lum_1', 'rad_1', +BCM_COLUMNS = ['tphys', 'kstar_1', 'mass0_1', 'mass_1', 'lum_1', 'rad_1', 'teff_1', 'massc_he_1', 'massc_co_1', 'radc_1', 'menv_1', 'renv_1', 'epoch_1', 'omega_spin_1', 'deltam_1', 'RRLO_1', 'kstar_2', 'mass0_2', 'mass_2', 'lum_2', 'rad_2', 'teff_2', 'massc_he_2', 'massc_co_2', 'radc_2', 'menv_2', From 04f0dd5282eb311bcb0ba79fbdc9e81c95c79ccc Mon Sep 17 00:00:00 2001 From: Katie Breivik Date: Tue, 16 Dec 2025 09:40:30 -0800 Subject: [PATCH 15/19] fixing ordering of core mass calcs on TP-AGB --- src/cosmic/src/evolv2.f | 3 +++ src/cosmic/src/hrdiag.f | 15 +++++++++------ 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/src/cosmic/src/evolv2.f b/src/cosmic/src/evolv2.f index 095ce35df..48f48942d 100644 --- a/src/cosmic/src/evolv2.f +++ b/src/cosmic/src/evolv2.f @@ -1335,6 +1335,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, b02_bcm = B(2) endif * Load preSN values for the SN writetab +* KB fix this mass0(k) = mass_preSN m0 = mass_preSN menv(k) = mHe_preSN @@ -1381,6 +1382,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, endif * Load preSN values for the SN writetab +* KB fix this too mass0(k) = mass_preSN m0 = mass_preSN menv(k) = mHe_preSN @@ -3659,6 +3661,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, b02_bcm = B(2) endif * Load preSN values for the SN writetab +* KB fix this too mass0(k) = mass_preSN m0 = mass_preSN menv(k) = mHe_preSN diff --git a/src/cosmic/src/hrdiag.f b/src/cosmic/src/hrdiag.f index af56901d1..91de4c2df 100644 --- a/src/cosmic/src/hrdiag.f +++ b/src/cosmic/src/hrdiag.f @@ -467,13 +467,16 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, * Approximate 3rd Dredge-up on AGB by limiting Mc. * lambdahrdiag = MIN(0.9d0,0.3d0+0.001d0*mass**5) +* Tau is the time at the start of the TP-AGB tau = tscls(13) - mcy = mcgbtf(tau,GB(2),GB,tscls(10),tscls(11),tscls(12)) - mcx = mc - mc = mcy - lambdahrdiag*(mcx-mcy) - mcx = mc - mc_co(kidx) = mcx - mc_he(kidx) = mcy - mcx +* mcx is M_c,DU in the equation *between* 73 and 74 of Hurley et al. 2000 + mcx = mcgbtf(tau,GB(2),GB,tscls(10),tscls(11),tscls(12)) +* mcy is M_c' in the same equation; it is defined in line 464 above for the current age. + mcy = mc +* The current core mass is then M_c' - lambda*(M_c' - M_c,DU) + mc = mcy - lambdahrdiag*(mcy-mcx) + mc_co(kidx) = mc + mc_he(kidx) = mc mcmax = MIN(mt,mcmax) endif r = ragbf(mt,lum,zpars(2)) From 6adfa19024aca01396e94e2053ee8309f5ab3688 Mon Sep 17 00:00:00 2001 From: Katie Breivik Date: Wed, 17 Dec 2025 12:05:47 -0800 Subject: [PATCH 16/19] change name to reflect core mass --- src/cosmic/evolve.py | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/src/cosmic/evolve.py b/src/cosmic/evolve.py index 2d9140ffe..40d1162c4 100644 --- a/src/cosmic/evolve.py +++ b/src/cosmic/evolve.py @@ -49,7 +49,7 @@ ALL_COLUMNS = ['tphys', 'mass_1', 'mass_2', 'kstar_1', 'kstar_2', 'sep', 'porb', 'ecc', 'RRLO_1', 'RRLO_2', 'evol_type', 'aj_1', 'aj_2', 'tms_1', - 'tms_2', 'massc_he_1', 'massc_he_2', 'massc_co_1', 'massc_co_2', + 'tms_2', 'massc_he_layer_1', 'massc_he_layer_2', 'massc_co_layer_1', 'massc_co_layer_2', 'rad_1', 'rad_2', 'mass0_1', 'mass0_2', 'lum_1', 'lum_2', 'teff_1', 'teff_2', 'radc_1', 'radc_2', 'menv_1', 'menv_2', 'renv_1', 'renv_2', 'omega_spin_1', @@ -63,7 +63,7 @@ BPP_COLUMNS = ['tphys', 'mass_1', 'mass_2', 'kstar_1', 'kstar_2', 'sep', 'porb', 'ecc', 'RRLO_1', 'RRLO_2', 'evol_type', 'aj_1', 'aj_2', 'tms_1', 'tms_2', - 'massc_he_1', 'massc_he_2', 'massc_co_1', 'massc_co_2', 'rad_1', 'rad_2', + 'massc_he_layer_1', 'massc_he_layer_2', 'massc_co_layer_1', 'massc_co_layer_2', 'rad_1', 'rad_2', 'mass0_1', 'mass0_2', 'lum_1', 'lum_2', 'teff_1', 'teff_2', 'radc_1', 'radc_2', 'menv_1', 'menv_2', 'renv_1', 'renv_2', 'omega_spin_1', 'omega_spin_2', 'B_1', 'B_2', 'bacc_1', 'bacc_2', @@ -71,9 +71,9 @@ 'bhspin_1', 'bhspin_2'] BCM_COLUMNS = ['tphys', 'kstar_1', 'mass0_1', 'mass_1', 'lum_1', 'rad_1', - 'teff_1', 'massc_he_1', 'massc_co_1', 'radc_1', 'menv_1', 'renv_1', 'epoch_1', + 'teff_1', 'massc_he_layer_1', 'massc_co_layer_1', 'radc_1', 'menv_1', 'renv_1', 'epoch_1', 'omega_spin_1', 'deltam_1', 'RRLO_1', 'kstar_2', 'mass0_2', 'mass_2', - 'lum_2', 'rad_2', 'teff_2', 'massc_he_2', 'massc_co_2', 'radc_2', 'menv_2', + 'lum_2', 'rad_2', 'teff_2', 'massc_he_layer_2', 'massc_co_layer_2', 'radc_2', 'menv_2', 'renv_2', 'epoch_2', 'omega_spin_2', 'deltam_2', 'RRLO_2', 'porb', 'sep', 'ecc', 'B_1', 'B_2', 'SN_1', 'SN_2', 'bin_state', 'merger_type'] @@ -369,18 +369,14 @@ def evolve(self, initialbinarytable, pool=None, bpp_columns=None, bcm_columns=No for i in range(len(initial_conditions)): initial_conditions[i]["n_col_bpp"] = len(bpp_columns) initial_conditions[i]["col_inds_bpp"] = col_inds_bpp - print("n_col_bpp: ", initial_conditions[0]["n_col_bpp"]) - print("col_inds_bpp: ", initial_conditions[0]["col_inds_bpp"], len(initial_conditions[0]["col_inds_bpp"])) - + # same for bcm col_inds_bcm = np.zeros(len(ALL_COLUMNS), dtype=int) col_inds_bcm[:len(bcm_columns)] = [ALL_COLUMNS.index(col) + 1 for col in bcm_columns] for i in range(len(initial_conditions)): initial_conditions[i]["n_col_bcm"] = len(bcm_columns) initial_conditions[i]["col_inds_bcm"] = col_inds_bcm - print("n_col_bcm: ", initial_conditions[0]["n_col_bcm"]) - print("col_inds_bcm: ", initial_conditions[0]["col_inds_bcm"], len(initial_conditions[0]["col_inds_bcm"])) - + # check if a pool was passed if pool is None: with MultiPool(processes=nproc) as pool: From 4514e770160d38b2a8e42310d09d55d37311c8d5 Mon Sep 17 00:00:00 2001 From: Katie Breivik Date: Wed, 17 Dec 2025 12:07:03 -0800 Subject: [PATCH 17/19] remove preSN variables in favor of new const_bse variables; also fix logging so that core masses are available at core collapse --- src/cosmic/src/evolv2.f | 52 +++++++++++++++++++---------------------- src/cosmic/src/hrdiag.f | 14 ++++------- 2 files changed, 28 insertions(+), 38 deletions(-) diff --git a/src/cosmic/src/evolv2.f b/src/cosmic/src/evolv2.f index 48f48942d..604494bb2 100644 --- a/src/cosmic/src/evolv2.f +++ b/src/cosmic/src/evolv2.f @@ -195,8 +195,6 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, REAL*8 deltam1_bcm,deltam2_bcm,b01_bcm,b02_bcm REAL*8 B(2),Bbot,omdot,b_mdot,b_mdot_lim,evolve_type COMMON /fall/fallback - REAL*8 mass_preSN, mHe_preSN, massc_preSN - COMMON mass_preSN, mHe_preSN, massc_preSN REAL ran3 EXTERNAL ran3 * @@ -414,7 +412,6 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, CALL star(kstar(k),mass0(k),mass(k),tm,tn,tscls,lums,GB,zpars) CALL hrdiag(mass0(k),age,mass(k),tm,tn,tscls,lums,GB,zpars, & rm,lum,kstar(k),mc,rc,me,re,k2,bhspin(k),k) - WRITE(*,*)'k=',k,' kstar=',kstar(k),'HRDIAG ran' aj(k) = age epoch(k) = tphys - age rad(k) = rm @@ -1217,7 +1214,6 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, * goto 140 endif * -* WRITE(*,*)'hrdiag call 1220' CALL star(kw,m0,mt,tm,tn,tscls,lums,GB,zpars) CALL hrdiag(m0,age,mt,tm,tn,tscls,lums,GB,zpars, & rm,lum,kw,mc,rc,me,re,k2,bhspin(k),k) @@ -1225,7 +1221,6 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, if(kw.ne.15)then ospin(k) = jspin(k)/(k2*(mt-mc)*rm*rm+k3*mc*rc*rc) endif -* WRITE(*,*)'hrdiag call finished' * * At this point there may have been a supernova. * @@ -1334,12 +1329,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, else b02_bcm = B(2) endif -* Load preSN values for the SN writetab -* KB fix this - mass0(k) = mass_preSN - m0 = mass_preSN - menv(k) = mHe_preSN - massc(k) = massc_preSN + CALL writetab(jp,tphys,evolve_type, & mass(1),mass(2),kstar(1),kstar(2), & sep,tb,ecc,rrl1,rrl2, @@ -1381,12 +1371,6 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, b02_bcm = B(2) endif -* Load preSN values for the SN writetab -* KB fix this too - mass0(k) = mass_preSN - m0 = mass_preSN - menv(k) = mHe_preSN - massc(k) = massc_preSN CALL writetab(jp,tphys,evolve_type, & mass(1),mass(2),kstar(1),kstar(2), & sep,tb,ecc,rrl1,rrl2, @@ -1498,7 +1482,6 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, * if(s.ge.-2.0457d0.or.s.le.-2.53d0) goto 175 if(s.ge.-1.6457d0.or.s.le.-2.53d0) goto 175 ospin(k) = (twopi*yearsc)/(10.d0**s)!have commented this out to keeps same spin -* write(*,*)'P=',s 176 u1 = ran3(idum1) u2 = ran3(idum1) if(u1.gt.0.9999d0) u1 = 0.9999d0 @@ -1586,7 +1569,17 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, if(snova)then * We should capture to evol_type change for SN as an evolutionary change evolve_type = 2.d0 - endif + endif + +* KB: set core masses to zero for remnants + if(kstar(1).ge.10)then + mc_he(1) = 0 + mc_co(1) = 0 + endif + if(kstar(2).ge.10)then + mc_he(2) = 0 + mc_co(2) = 0 + endif rrl1 = rad(1)/rol(1) rrl2 = rad(2)/rol(2) teff1 = 1000.d0*((1130.d0*lumin(1)/ @@ -1623,7 +1616,6 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, & deltam1_bcm,deltam2_bcm,formation(1), & formation(2),binstate,mergertype,z,'bpp') if(snova)then - bpp(jp,11) = 2.0 dtm = 0.d0 goto 4 endif @@ -1697,7 +1689,6 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, * If not interpolating set the next timestep. * if(intpol.eq.0)then -* WRITE(*,*)'you should see this to advance the time' if(output) write(*,*)'nxt t, prior:',tphys,dtm,dtmi(1),dtmi(2) dtm = MAX(1.0d-07*tphys,MIN(dtmi(1),dtmi(2))) dtm = MIN(dtm,tsave-tphys) @@ -1801,6 +1792,17 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, if(change)then change = .false. evolve_type = 2.d0 + +* KB: set core masses to zero for remnants + if(kstar(1).ge.10)then + mc_he(1) = 0 + mc_co(1) = 0 + endif + if(kstar(2).ge.10)then + mc_he(2) = 0 + mc_co(2) = 0 + endif + rrl1 = rad(1)/rol(1) rrl2 = rad(2)/rol(2) teff1 = 1000.d0*((1130.d0*lumin(1)/ @@ -3156,7 +3158,6 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, do 602 , k = 1,2 * dms(k) = km*dms(k) -* WRITE(*,*)dme/tb,dms(j2)/tb/km,dmt(j2),dms(j1)/tb/km,dmr(j1) if(kstar(k).lt.10) dms(k) = MIN(dms(k),mass(k) - massc(k)) * @@ -3660,12 +3661,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, else b02_bcm = B(2) endif -* Load preSN values for the SN writetab -* KB fix this too - mass0(k) = mass_preSN - m0 = mass_preSN - menv(k) = mHe_preSN - massc(k) = massc_preSN + CALL writetab(jp,tphys,evolve_type, & mass(1),mass(2),kstar(1),kstar(2), & sep,tb,ecc,rrl1,rrl2, diff --git a/src/cosmic/src/hrdiag.f b/src/cosmic/src/hrdiag.f index 91de4c2df..85e0ec90b 100644 --- a/src/cosmic/src/hrdiag.f +++ b/src/cosmic/src/hrdiag.f @@ -180,8 +180,6 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, endif eta = mctmsf(mass) tau = (aj - tm)/thg -* WRITE(*,*)'hrdiag: 184, k=',kidx,' mass=',mass -* WRITE(*,*)'mc=',mc,' mcx=',mcx mc = ((1.d0 - tau)*eta + tau)*mc mc = MAX(mc,mcx) @@ -256,9 +254,7 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, mcx = mcheif(mass,zpars(2),zpars(9)) mcy = mcheif(mass,zpars(2),zpars(10)) mc = mcx + (mcy - mcx)*tau -* WRITE(*,*)'hrdiag249: k=',kidx,'mass=',mass,'kw=',kw -* WRITE(*,*)'mc=',mc,' mcx=',mcx,' mcy-mcx=',mcy-mcx -* WRITE(*,*) + mc_he(kidx) = mc mc_co(kidx) = 0.0 endif @@ -305,8 +301,6 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, tau = (aj - tscls(2))/tscls(3) * here, mcx is the helium core mass at helium ignition mc = mcx + (mcagbf(mass) - mcx)*tau -* WRITE(*,*)'hrdiag: mc=',mc,' mcx=',mcx, 'kw=',kw,' k=',kidx -* WRITE(*,*)'hrdiag: mc_he=',(mcagbf(mass) - mcx)*tau mc_he(kidx) = mc mc_co(kidx) = 0.0 * @@ -429,8 +423,6 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, if(aj.lt.tscls(13))then mcx = mcgbtf(aj,GB(8),GB,tscls(7),tscls(8),tscls(9)) mc = mcbagb -* WRITE(*,*)'hrdiag 430: mc=',mc,' mcx=',mcx,' kw=',kw -* WRITE(*,*)'hrdiag: mcbagb=',mcbagb mc_co(kidx) = mcx mc_he(kidx) = mcbagb - mcx lum = lmcgbf(mcx,GB) @@ -475,8 +467,9 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, mcy = mc * The current core mass is then M_c' - lambda*(M_c' - M_c,DU) mc = mcy - lambdahrdiag*(mcy-mcx) + mcx = mc mc_co(kidx) = mc - mc_he(kidx) = mc + mc_he(kidx) = 0.0 mcmax = MIN(mt,mcmax) endif r = ragbf(mt,lum,zpars(2)) @@ -649,6 +642,7 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, k2 = 0.1d0 endif endif + * C if(mass.gt.99.99d0)then C mass = mass0 From fd37b3cb29ff09d21e35b8803817f183277094da Mon Sep 17 00:00:00 2001 From: Katie Breivik Date: Mon, 5 Jan 2026 13:49:58 -0700 Subject: [PATCH 18/19] this should fix things but Im still gettin local errors --- meson.build | 8 +- src/cosmic/src/comenv.f | 152 +++++++++--------- src/cosmic/src/comprad.f | 6 + src/cosmic/tests/data/unit_tests_results.hdf5 | Bin 8567500 -> 9620252 bytes 4 files changed, 88 insertions(+), 78 deletions(-) diff --git a/meson.build b/meson.build index 266c9aa4a..c3dd26277 100644 --- a/meson.build +++ b/meson.build @@ -37,20 +37,16 @@ lib_source = [ 'src/cosmic/src/mrenv.f', 'src/cosmic/src/ran3.f', 'src/cosmic/src/rl.f', + 'src/cosmic/src/hrdiag.f', + 'src/cosmic/src/star.f', 'src/cosmic/src/concatkstars.f', 'src/cosmic/src/comprad.f', 'src/cosmic/src/bpp_array.f', 'src/cosmic/src/checkstate.f', 'src/cosmic/src/deltat.f', 'src/cosmic/src/mlwind.f', - 'src/cosmic/src/hrdiag.f', - 'src/cosmic/src/star.f', 'src/cosmic/src/zcnsts.f', 'src/cosmic/src/deltat.f', - 'src/cosmic/src/mlwind.f', - 'src/cosmic/src/hrdiag.f', - 'src/cosmic/src/star.f', - 'src/cosmic/src/zcnsts.f', 'src/cosmic/src/zfuncs.f'] # Detect operating system and set appropriate linker flags diff --git a/src/cosmic/src/comenv.f b/src/cosmic/src/comenv.f index 25b17eb59..ad1bcfe37 100644 --- a/src/cosmic/src/comenv.f +++ b/src/cosmic/src/comenv.f @@ -321,15 +321,16 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, & kstar2_bpp,SEP_postCE,TB,ECC, & rrl1_bpp,rrl2_bpp, & aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp, - & massc1_bpp,massc_preSN,rad1_bpp,rad2_bpp, - & M02,mass_preSN,lumin(1),lumin(2), - & teff1,teff2, - & RC2,RC1,menv_bpp(1),mHe_preSN,renv_bpp(1), - & renv_bpp(2),OSPIN2,OSPIN1,B_0(1),B_0(2), - & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), - & epoch(2),bhspin2,bhspin1, - & deltam_2,deltam_1,formation2,formation1, - & binstate,mergertype,'bpp') + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad1_bpp,rad2_bpp, + & M02,mass_preSN,lumin(1),lumin(2), + & teff1,teff2, + & RC2,RC1,menv_bpp(1),mHe_preSN,renv_bpp(1), + & renv_bpp(2),OSPIN2,OSPIN1,B_0(1),B_0(2), + & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), + & epoch(2),bhspin2,bhspin1, + & deltam_2,deltam_1,formation2,formation1, + & binstate,mergertype,zpars(14)**2.d5,'bpp') else teff1 = 1000.d0*((1130.d0*lumin(1)/ & (rad1_bpp**2.d0))**(1.d0/4.d0)) @@ -340,15 +341,16 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, & kstar2_bpp,SEP_postCE,TB,ECC, & rrl1_bpp,rrl2_bpp, & aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp, - & massc_preSN,massc2_bpp,rad1_bpp,rad2_bpp, - & mass_preSN,M02,lumin(1),lumin(2), - & teff1,teff2, - & RC1,RC2,mHe_preSN,menv_bpp(2),renv_bpp(1), - & renv_bpp(2),OSPIN1,OSPIN2,B_0(1),B_0(2), - & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), - & epoch(2),bhspin1,bhspin2, - & deltam_1,deltam_2,formation1,formation2, - & binstate,mergertype,'bpp') + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad1_bpp,rad2_bpp, + & mass_preSN,M02,lumin(1),lumin(2), + & teff1,teff2, + & RC1,RC2,mHe_preSN,menv_bpp(2),renv_bpp(1), + & renv_bpp(2),OSPIN1,OSPIN2,B_0(1),B_0(2), + & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), + & epoch(2),bhspin1,bhspin2, + & deltam_1,deltam_2,formation1,formation2, + & binstate,mergertype,zpars(14)**2.d5,'bpp') endif endif CALL kick(KW1,M_postCE,M1,M2,ECC,SEP_postCE, @@ -628,15 +630,16 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, & kstar2_bpp,SEP_postCE,TB,ECC, & rrl1_bpp,rrl2_bpp, & aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp, - & massc1_bpp,massc_preSN,rad1_bpp,rad2_bpp, - & M02,mass_preSN,lumin(1),lumin(2), - & teff1,teff2, - & RC2,RC1,menv_bpp(1),mHe_preSN,renv_bpp(1), - & renv_bpp(2),OSPIN2,OSPIN1,B_0(1),B_0(2), - & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), - & epoch(2),bhspin2,bhspin1, - & deltam_2,deltam_1,formation2,formation1, - & binstate,mergertype,'bpp') + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad1_bpp,rad2_bpp, + & M02,mass_preSN,lumin(1),lumin(2), + & teff1,teff2, + & RC2,RC1,menv_bpp(1),mHe_preSN,renv_bpp(1), + & renv_bpp(2),OSPIN2,OSPIN1,B_0(1),B_0(2), + & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), + & epoch(2),bhspin2,bhspin1, + & deltam_2,deltam_1,formation2,formation1, + & binstate,mergertype,zpars(14)**2.d5,'bpp') else teff1 = 1000.d0*((1130.d0*lumin(1)/ & (rad1_bpp**2.d0))**(1.d0/4.d0)) @@ -647,15 +650,16 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, & kstar2_bpp,SEP_postCE,TB,ECC, & rrl1_bpp,rrl2_bpp, & aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp, - & massc_preSN,massc2_bpp,rad1_bpp,rad2_bpp, - & mass_preSN,M02,lumin(1),lumin(2), - & teff1,teff2, - & RC1,RC2,mHe_preSN,menv_bpp(2),renv_bpp(1), - & renv_bpp(2),OSPIN1,OSPIN2,B_0(1),B_0(2), - & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), - & epoch(2),bhspin1,bhspin2, - & deltam_1,deltam_2,formation1,formation2, - & binstate,mergertype,'bpp') + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad1_bpp,rad2_bpp, + & mass_preSN,M02,lumin(1),lumin(2), + & teff1,teff2, + & RC1,RC2,mHe_preSN,menv_bpp(2),renv_bpp(1), + & renv_bpp(2),OSPIN1,OSPIN2,B_0(1),B_0(2), + & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), + & epoch(2),bhspin1,bhspin2, + & deltam_1,deltam_2,formation1,formation2, + & binstate,mergertype,zpars(14)**2.d5,'bpp') endif endif * USSN: if ussn flag is set, have reduced kicks for stripped He stars (SN=8) @@ -799,15 +803,16 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, & kstar2_bpp,SEP_postCE,TB,ECC, & rrl1_bpp,rrl2_bpp, & aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp, - & massc_preSN,massc2_bpp,rad1_bpp,rad2_bpp, - & mass_preSN,M01,lumin(1),lumin(2), - & teff1,teff2, - & RC2,RC1,mHe_preSN,menv_bpp(2),renv_bpp(1), - & renv_bpp(2),OSPIN2,OSPIN1,B_0(1),B_0(2), - & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), - & epoch(2),bhspin2,bhspin1, - & deltam_2,deltam_1,formation2,formation1, - & binstate,mergertype,'bpp') + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad1_bpp,rad2_bpp, + & mass_preSN,M01,lumin(1),lumin(2), + & teff1,teff2, + & RC2,RC1,mHe_preSN,menv_bpp(2),renv_bpp(1), + & renv_bpp(2),OSPIN2,OSPIN1,B_0(1),B_0(2), + & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), + & epoch(2),bhspin2,bhspin1, + & deltam_2,deltam_1,formation2,formation1, + & binstate,mergertype,zpars(14)**2.d5,'bpp') else teff1 = 1000.d0*((1130.d0*lumin(1)/ & (rad1_bpp**2.d0))**(1.d0/4.d0)) @@ -818,15 +823,16 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, & kstar2_bpp,SEP_postCE,TB,ECC, & rrl1_bpp,rrl2_bpp, & aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp, - & massc1_bpp,massc_preSN,rad1_bpp,rad2_bpp, - & M01,mass_preSN,lumin(1),lumin(2), - & teff1,teff2, - & RC1,RC2,menv_bpp(1),mHe_preSN,renv_bpp(1), - & renv_bpp(2),OSPIN1,OSPIN2,B_0(1),B_0(2), - & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), - & epoch(2),bhspin1,bhspin2, - & deltam_1,deltam_2,formation1,formation2, - & binstate,mergertype,'bpp') + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad1_bpp,rad2_bpp, + & M01,mass_preSN,lumin(1),lumin(2), + & teff1,teff2, + & RC1,RC2,menv_bpp(1),mHe_preSN,renv_bpp(1), + & renv_bpp(2),OSPIN1,OSPIN2,B_0(1),B_0(2), + & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), + & epoch(2),bhspin1,bhspin2, + & deltam_1,deltam_2,formation1,formation2, + & binstate,mergertype,zpars(14)**2.d5,'bpp') endif endif CALL kick(KW2,M_postCE,M2,M1,ECC,SEP_postCE, @@ -1033,15 +1039,16 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, & kstar2_bpp,-1.d0,TB,0.d0, & rrl1_bpp,rrl2_bpp, & aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp, - & massc1_bpp,massc_preSN,rad1_bpp,rad2_bpp, - & M02,mass_preSN,lumin(1),lumin(2), - & teff1,teff2, - & RC2,RC1,menv_bpp(1),mHe_preSN,renv_bpp(1), - & renv_bpp(2),OSPIN2,OSPIN1,B_0(1),B_0(2), - & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), - & epoch(2),bhspin2,bhspin1, - & deltam_2,deltam_1,formation2,formation1, - & binstate,mergertype,'bpp') + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad1_bpp,rad2_bpp, + & M02,mass_preSN,lumin(1),lumin(2), + & teff1,teff2, + & RC2,RC1,menv_bpp(1),mHe_preSN,renv_bpp(1), + & renv_bpp(2),OSPIN2,OSPIN1,B_0(1),B_0(2), + & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), + & epoch(2),bhspin2,bhspin1, + & deltam_2,deltam_1,formation2,formation1, + & binstate,mergertype,zpars(14)**2.d5,'bpp') else teff1 = 1000.d0*((1130.d0*lumin(1)/ & (rad1_bpp**2.d0))**(1.d0/4.d0)) @@ -1052,15 +1059,16 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, & kstar2_bpp,-1.d0,TB,0.d0, & rrl1_bpp,rrl2_bpp, & aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp, - & massc_preSN,massc2_bpp,rad1_bpp,rad2_bpp, - & mass_preSN,M02,lumin(1),lumin(2), - & teff1,teff2, - & RC1,RC2,mHe_preSN,menv_bpp(2),renv_bpp(1), - & renv_bpp(2),OSPIN1,OSPIN2,B_0(1),B_0(2), - & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), - & epoch(2),bhspin1,bhspin2, - & deltam_1,deltam_2,formation1,formation2, - & binstate,mergertype,'bpp') + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad1_bpp,rad2_bpp, + & mass_preSN,M02,lumin(1),lumin(2), + & teff1,teff2, + & RC1,RC2,mHe_preSN,menv_bpp(2),renv_bpp(1), + & renv_bpp(2),OSPIN1,OSPIN2,B_0(1),B_0(2), + & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), + & epoch(2),bhspin1,bhspin2, + & deltam_1,deltam_2,formation1,formation2, + & binstate,mergertype,zpars(14)**2.d5,'bpp') endif endif CALL kick(KW,MF,M1,0.d0,0.d0,-1.d0,0.d0,vk,star1, diff --git a/src/cosmic/src/comprad.f b/src/cosmic/src/comprad.f index bfbf3dc5c..287cab16e 100644 --- a/src/cosmic/src/comprad.f +++ b/src/cosmic/src/comprad.f @@ -34,7 +34,13 @@ SUBROUTINE compute_r(mass,z,num,rad) *** * Then just loop through everything *** + + if(num.gt.loop) then + print *, 'Error in compute_r: num > 100000' + stop + end if do 10 , k = 1,num + WRITE(*,*)' Computing radius for star ',k,' of ',num age = 0.0 mc = 0.d0 tm = 0.d0 diff --git a/src/cosmic/tests/data/unit_tests_results.hdf5 b/src/cosmic/tests/data/unit_tests_results.hdf5 index 138c564f5ad91af20ccfcac4f9c2c76a64bc345c..7ed44481526185871e2015b58202410ef56b9fe0 100644 GIT binary patch delta 4476 zcmeH~eN0=|6~KL;`FOU0Jd7z0;luIPgb+R+CLw_Ud3*&gO9LfY<8-XXU{liA*chi; zq%7ScRYYqwVQ!Smq-dA`A<0M+&z)>gHF@18b%k}F7HyeSZAf#asN1A$(l*vLV&{E1 zuU?{Q+V;n$J?Y2uIq#l(?>X6+i9OYMv zV^ulQy%PU|ZqpvWDlQV=LfARVDA-0g4Xwjm7p=QqCk}=M$iWb`kDQO1ci>4iiy1Y8 zHw?wNvc7>-#Vnz@`-R06+f8cDp{J%;E4CV@@qmHB>joOdg5lfzba?U@t$qMg%`Eq| zU{hTawpd;uZ@LkqzT`&1a4BXpXsTjQFkD7*Jcy<)a_|#>8I-mLUfItpt)09w0m^}P zUbzHHb;s{^Nhs+q%Z!QL=wa926GHACfV?~ICrOq3WtejNyXj2vH}m`Hpd zLr?Y4_&#kUoiPiim#28pA)47#q-=rOhF#j-c%Yuf@gAc_Ue1gi9R-LHUaNL(#V=SW z*8zXtMO4CvmoiO)gL7;H()xOpQ(N&$PZ_p#^ve-=9$|e~1I6@$4|UXfH~vyrx%e#9 zHGm-hE=1JbQ5rYZ?ub0-wG@;W36U|1lR^Z<7+4LWA{6|}PR*o6S~48|Gm4B1e?K#q zSDYRXIQ(|VAo6So6;$&ok7InC*T|Rqj7KS|{0^rmV~4xo6@CiDkLyK;ecX3yR8%@$?tnwwRI9MdzP=+pit&`|E+BYOMBd7xRaoT)Oc3ZX@pWnzDN~ zPWrMfmH70N1}%B5YA91;KZ6Q`aU-4|JAi8&8*zn>A?`DX2Ho8J*hhTb=iRFprx&Y% zG_4Fa;mNvCs>*grtM5`)naM3`Cl%EURc1q^bDpXo9|sT}y)ufDY%@PDJ(&er5z$os zOG}BeMToEvT`A<)-cl$`n4W6sugLYX*KH4+^|_=89WP0*Cy`?iDlXuWs*51Vy=UZ2 z{(b$V>M1?%((8-QLurwpG#8NLA;gr*(<-w*Rmfu2JKc}Z`t)jjY?dp|hsS2U;(2-% zD6+mp7^YANZ$^>O#4wA*T@E^ihLp7#;AO=sf~1RlK?G27RI#2juNpflizJf;2 zc=$*7L|!h`4@0iZB3m%d;6#ZjYPgC_D7ef-_sk(q_z{>#EGl~A2eq5oKT~nMA%9K( z@^WonizfQb1=NW0-`V`wz7*LwC_aN#$p_a^xmagaC z)(r=IH;$^?Qr4?CmDXQFZyMu^n<+o*e*E6$Qy*q<^4hgMn-Q|`TFP3t1lC*S9me0f z^Os!fpHkM-r1};r+9GX~lg=P=>=s08h7o^k!&KzPq4*np*$l-YU>U)uJ9K01QAO@Dxx4 z7=f)oF|Z9N0ZM@~pd2s(+oM&B=KYN3oYeSOCp=Lz@G0Qaz|+7EU?)%k>;fu*-9Qyk z4b%X&01MO+Pt@WSJ`EdE?;?I0aBHC-C$x;FMoN(y+}2ivS=xl}=!_F;d)y z4ZuENKhOv?0SAC)pap0J4gws|2DAemKqt@z90IxlEAScMFz{L62=EMW6zBnZfoI8D ztnVhdw%pf1?#@vRHWivC-e2wmjseGk&jHT?Cjc985_lf?JkSq(0eAth0|S7AyuUnX zBv%&nG`cwPaCs0I0-S&g7zSPhMu1bmOTcNs4S0Z2zzg_*G4gQP{|X;e;*6S$$@%5J z9r;J4Kn8pxXUF-N{+G6Yq8sVoWd0+s{P3kd`u}6?93760$!D*`I{$ zPeS%5A^Ve%{Yl9FBxHXQvOfvgpM>mBLiQ(PTr?SLd*>lGxpPYRBzRJ@#J8L2%}^C@ zkRJH-X6XOT&%fOK4DfOOPc}bziSb3)?u<-4Y0;ufWMakNMVvE~X^N%sdXGiJj`3Si rsdS`vj@l~!k7`4aJHx)rU_6#R4-bF)mS1T3BMJZAjeqZ-dH>%4y#M-6 delta 2255 zcmb`|drVVT90%}wE(LB&#Zst%ye^1P6{Uc@lt)4AI>#u(%_kG)lSvV%Y{?R`na%~G zSyDKQ$3Gm(OVK!Ot;>bZKbW#C29jaIqDz)&7MHk(i(639VA!|4Ejs^LvSe*OJ-6rn z?z!jO-#O==I5w`Bcz8tN&<9O;>|HB9A)PMQ9Eo#s^5?j}G)Kt?gwwgl`7e28#&EW; zrjxqzlIgkN6h&pSxq(gw&kzT_9+FNv-G7GMBoTUl18RS%r#}c~s>8iR{=bkJ`uP~& zwQi$Q24Aize1%C5yTZ=12{v@in8# zkz1?gIM&=m^4W({N*Z+pWwZB3Y6|GntO8cz*_p?t8_6kZ$gXD-SM4hHhm{EY6k&I+ z+T*FFyqZ=E0&mP^7aDQvr#Y7t2QebIb9Fu8_4#z0P(rKoR4lEXXz2dDMG>NWPmcGO zqHko=O@V58WZpKYQ8_~mMuD%l&^5*?>QEWoqxl+v?LNti%0N3Jbr~axn%%IHgrMnW zj=6gIq(C2{2y>3vD`1KRppkUB5~x9e7NJR1TGSgB$#<3|qYQ z&K!Z)toH^FScyiR?my`*nK))|qlPen&9)IiI&py9HPEr}P2s)=;0vs&skFF<)kDPl z$gxo4MQOT^h(t;qAVwm!4H5$hT~f0v`lwPTeRY}?5>m^8he%Q&18l_*DHNvWX6w$o zXG)oMh}`=3h%0Bwt@J?3Ov*j0ejOK#QeDZ`tlIk}s?p6eVq&tFjkJz;f7~QEW&L<1Zk4l;1BALt#ToIUFfxZXImzt0wog{+%*Q?kr7jaS1sz$zI z=yZ5m*|k*6m^LV5J|X9+J939=8^t|T|;%)DNo!x1EQK^k2R$r<<*))WU?r@AOfWze-XhsM zd}~da;a_VRv^O`J8U#J9)*bJd7|w!xq>I z+h9A?!rQO|>Y!dyc;0z1r=0P>>K}AG8`%lFU^nc6y|51&;9b}cjbMc)H~`Ji0 Date: Mon, 5 Jan 2026 15:47:16 -0700 Subject: [PATCH 19/19] found the killer -- kidx was becoming huge becuase in comprad it is the iterator that goes up to 1e5 even though mc_co and mc_he are of size 2, so memory gets overwritten and eventually busses. rip --- src/cosmic/sample/sampler/independent.py | 4 +--- src/cosmic/src/comprad.f | 7 +++++-- src/cosmic/src/const_bse.h | 7 ++++--- 3 files changed, 10 insertions(+), 8 deletions(-) diff --git a/src/cosmic/sample/sampler/independent.py b/src/cosmic/sample/sampler/independent.py index fb6a052c5..bd926edb7 100644 --- a/src/cosmic/sample/sampler/independent.py +++ b/src/cosmic/sample/sampler/independent.py @@ -1158,10 +1158,8 @@ def set_reff(self, mass, metallicity, zsun=0.02): of length 10^5. If your masses are more than that, you'll need to divide it into chunks """ - from cosmic import _evolvebin - max_array_size = 100000 total_length = len(mass) radii = np.zeros(total_length) @@ -1183,7 +1181,7 @@ def set_reff(self, mass, metallicity, zsun=0.02): length_remaining = total_length - ## if smaller than 10^5, need to pad out the array + # if smaller than 10^5, need to pad out the array temp_mass = np.zeros(max_array_size) temp_mass[:length_remaining] = mass[-length_remaining:] diff --git a/src/cosmic/src/comprad.f b/src/cosmic/src/comprad.f index 287cab16e..ce7abb404 100644 --- a/src/cosmic/src/comprad.f +++ b/src/cosmic/src/comprad.f @@ -39,8 +39,11 @@ SUBROUTINE compute_r(mass,z,num,rad) print *, 'Error in compute_r: num > 100000' stop end if + mc_he(1) = 0.d0 + mc_he(2) = 0.d0 + mc_co(1) = 0.d0 + mc_co(2) = 0.d0 do 10 , k = 1,num - WRITE(*,*)' Computing radius for star ',k,' of ',num age = 0.0 mc = 0.d0 tm = 0.d0 @@ -55,7 +58,7 @@ SUBROUTINE compute_r(mass,z,num,rad) rc = 0.d0 CALL star(kstar,mass0,mt,tm,tn,tscls,lums,GB,zpars) CALL hrdiag(mass0,age,mt,tm,tn,tscls,lums,GB,zpars, - & rad(k),lum,kstar,mc,rc,me,re,k2,bhspin,k) + & rad(k),lum,kstar,mc,rc,me,re,k2,bhspin,1) 10 continue diff --git a/src/cosmic/src/const_bse.h b/src/cosmic/src/const_bse.h index 2327a95d6..1eab290d7 100644 --- a/src/cosmic/src/const_bse.h +++ b/src/cosmic/src/const_bse.h @@ -20,8 +20,7 @@ INTEGER ceflag,cekickflag,cemergeflag,cehestarflag,ussn COMMON /CEFLAGS/ ceflag,cekickflag,cemergeflag,cehestarflag,ussn INTEGER pisn_track(2) - REAL*8 mc_he(2),mc_co(2) - COMMON /TRACKERS/ pisn_track,mc_he,mc_co + COMMON /TRACKERS/ pisn_track * REAL*8 zsun COMMON /METVARS/ zsun @@ -39,9 +38,11 @@ REAL*8 polar_kick_angle REAL*8 ecsn,ecsn_mlow,bhspinmag,rembar_massloss REAL*8 natal_kick_array(2,5) + REAL*8 mc_he(2),mc_co(2) COMMON /SNVARS/ natal_kick_array,sigma,sigmadiv,bhsigmafrac, & polar_kick_angle,pisn,ecsn,ecsn_mlow, - & bhspinmag,mxns,rembar_massloss,kickflag + & bhspinmag,mxns,rembar_massloss,kickflag, + & mc_he,mc_co REAL*8 fprimc_array(16) COMMON /TIDALVARS/ fprimc_array REAL*8 rejuv_fac