From 82e49c4e64de4ee7496dba666e29dbedbf91a35d Mon Sep 17 00:00:00 2001 From: hskobe Date: Fri, 10 Apr 2026 14:13:18 -0400 Subject: [PATCH 01/15] Updated comenv.f with alpha choice for binding energy --- src/cosmic/src/comenv.f | 22 +++++++++++++++++++--- 1 file changed, 19 insertions(+), 3 deletions(-) diff --git a/src/cosmic/src/comenv.f b/src/cosmic/src/comenv.f index 83e7ee744..aecc3328f 100644 --- a/src/cosmic/src/comenv.f +++ b/src/cosmic/src/comenv.f @@ -77,6 +77,7 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, REAL*8 KW1_TEMP, KW2_TEMP REAL*8 rad(2),tms(2),lumin(2),B_0(2),bacc(2),tacc(2),epoch(2) REAL*8 menv_bpp(2),renv_bpp(2) + REAL*8 ALPHA_CE * * Initialize * @@ -154,7 +155,12 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, * * Calculate the final orbital energy without coalescence. * - EORBF = EORBI + EBINDI/ALPHA1 + IF(switchedCE)THEN + ALPHA_CE = ALPHA1(2) + ELSE + ALPHA_CE = ALPHA1(1) + ENDIF + EORBF = EORBI + EBINDI/ALPHA_CE * * If the secondary is on the main sequence see if it fills its Roche lobe. * @@ -205,7 +211,12 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, * Coalescence - calculate final binding energy. * EORBF = MAX(MC1*M2/(2.D0*SEPL),EORBI) - EBINDF = EBINDI - ALPHA1*(EORBF - EORBI) + IF(switchedCE)THEN + ALPHA_CE = ALPHA1(2) + ELSE + ALPHA_CE = ALPHA1(1) + ENDIF + EBINDF = EBINDI - ALPHA_CE*(EORBF - EORBI) KW1_TEMP = KW KW2_TEMP = 15 ELSE @@ -441,7 +452,12 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, * Calculate the final envelope binding energy. * EORBF = MAX(MC1*MC2/(2.D0*SEPL),EORBI) - EBINDF = EBINDI - ALPHA1*(EORBF - EORBI) + IF(switchedCE)THEN + ALPHA_CE = ALPHA1(2) + ELSE + ALPHA_CE = ALPHA1(1) + ENDIF + EBINDF = EBINDI - ALPHA_CE*(EORBF - EORBI) if(output) write(*,*)'In dg or giant 1:',M01,M1,R1,M02,M2, & R2,MC1,MC2,MC3,KW1,KW2,KW,EORBF,EBINDF * From 147a4240730e9e41c32d2e45860454dfabebd928 Mon Sep 17 00:00:00 2001 From: hskobe Date: Fri, 10 Apr 2026 14:58:07 -0400 Subject: [PATCH 02/15] Updated alpha/acclim in const_bse.h --- src/cosmic/src/const_bse.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/cosmic/src/const_bse.h b/src/cosmic/src/const_bse.h index 997e97fed..11588b856 100644 --- a/src/cosmic/src/const_bse.h +++ b/src/cosmic/src/const_bse.h @@ -15,7 +15,7 @@ & qcflag,eddlimflag,bhspinflag,aic,rejuvflag, & htpmb,ST_cr,ST_tide,bdecayfac,grflag, & bhms_coll_flag,wd_mass_lim,rtmsflag,maltsev_mode - REAL*8 don_lim,acc_lim,Mbh_initial + REAL*8 don_lim,acc_lim(2),Mbh_initial COMMON /MTVARS/ don_lim,acc_lim,Mbh_initial INTEGER ceflag,cekickflag,cemergeflag,cehestarflag,ussn COMMON /CEFLAGS/ ceflag,cekickflag,cemergeflag,cehestarflag,ussn @@ -29,7 +29,7 @@ INTEGER LBV_flag COMMON /WINDVARS/ neta,bwind,hewind,beta,xi,acc2,epsnov, & eddfac,gamma,LBV_flag - REAL*8 alpha1,lambdaf + REAL*8 alpha1(2),lambdaf REAL*8 qcrit_array(16) COMMON /CEVARS/ qcrit_array,alpha1,lambdaf REAL*8 bconst,CK @@ -60,7 +60,7 @@ COMMON /POINTS/ pts1,pts2,pts3 REAL*8 dmmax,drmax COMMON /TSTEPC/ dmmax,drmax - REAL*8 scm(50000,14),spp(20,3) + REAL*8 scm(50000,16),spp(25,20) COMMON /SINGLE/ scm,spp REAL*8 bcm(50000,52),bpp(1000,52) COMMON /BINARY/ bcm,bpp From 9e7ada68f45d7c9f0ba38703ca203258cec8f00f Mon Sep 17 00:00:00 2001 From: hskobe Date: Mon, 13 Apr 2026 12:20:18 -0400 Subject: [PATCH 03/15] Finished updates for evolv1.f --- src/cosmic/src/evolv1.f | 368 +++++++++++++++++++++++++++++++--------- 1 file changed, 291 insertions(+), 77 deletions(-) diff --git a/src/cosmic/src/evolv1.f b/src/cosmic/src/evolv1.f index 615ccfac4..780a65e75 100644 --- a/src/cosmic/src/evolv1.f +++ b/src/cosmic/src/evolv1.f @@ -34,16 +34,17 @@ SUBROUTINE evolv1(kw,mass,mt,r,lum,mc,rc,menv,renv,ospin, c c-------------------------------------------------------------c implicit none + INCLUDE 'const_bse.h' * - integer kw,it,ip,jp,j,kwold,rflag,fb - integer nv + integer kw,it,ip,jp,j,kwold,fb,k2str,pulsar + integer nv,intpol,k parameter(nv=50000) * real*8 mass,z,aj real*8 epoch,tphys,tphys2,tmold,tbgold real*8 mt,tm,tn,tphysf,dtp,tsave real*8 tscls(20),lums(10),GB(10),zpars(20) - real*8 r,lum,mc,teff,rc,menv,renv,kick_info(12) + real*8 r,lum,mc,teff,rc,menv,renv,kick_info(2,17) real*8 ospin,jspin,djt,djmb,k2,k3 parameter(k3=0.21d0) real*8 m0,r1,lum1,mc1,rc1,menv1,renv1,k21 @@ -54,17 +55,22 @@ SUBROUTINE evolv1(kw,mass,mt,r,lum,mc,rc,menv,renv,ospin, parameter(eps=1.0d-06,alpha2=0.09d0) real*8 mlwind,vrotf external mlwind,vrotf - logical iplot,isave + logical iplot,isave,disrupt common /fall/fallback -* - REAL*8 fallback,vk !PDK -* - REAL*8 neta,bwind,hewind,mxns - COMMON /VALUE1/ neta,bwind,hewind,mxns - REAL*8 pts1,pts2,pts3 - COMMON /POINTS/ pts1,pts2,pts3 - REAL scm(50000,14),spp(20,3) - COMMON /SINGLE/ scm,spp + REAL*8 kw3,wsun,wx + PARAMETER(kw3=619.2d0,wsun=9.46d+07,wx=9.46d+08) + REAL*8 yeardy,yearsc,aursun + PARAMETER(yeardy=365.24d0,aursun=214.95d0,yearsc=3.1557d+07) + REAL*8 twopi,u1,u2,mew,tacc + REAL*8 fallback,vk,B_0,B,ospbru,b01_bcm,b_mdot,b_mdot_lim + REAL*8 bacc,Bbot,convradcomp,deltam,dspint,eqspin + REAL*8 dtj,evolve_type + REAL*8 omdot,sigmahold,sn,bhspin,s + REAL*8 djtx,djspint,jspbru,Kconst + REAL*8 bkick(20) +* + REAL ran3 + EXTERNAL ran3 * dtm = 0.d0 r = 0.d0 @@ -78,15 +84,28 @@ SUBROUTINE evolv1(kw,mass,mt,r,lum,mc,rc,menv,renv,ospin, jspin = 1.0d-10 endif k2 = 0.15d0 - rflag = 0 + * PDK additions fb = 1 !turns kick velocity limit owning to fallback on (1) or off (0). vk = 0.d0 !gets passed to and from kick(), is kick mag. can be used to set initial pulsar particulars. + B_0 = 0.0 + k2str = 0 + if(kw.ne.13)then + bacc = 0.d0 + tacc = 0.d0 + endif + aic = 1 + b_mdot_lim = -1.0e-11 + Bbot = 5e+7 + Kconst = 2.5d-49 + + twopi = 2.d0*ACOS(-1.d0) + * * Setup variables which control the output (if it is required). * - ip = 0 - jp = 0 + ip = 1 + jp = 1 tsave = tphys isave = .true. iplot = .false. @@ -125,20 +144,146 @@ SUBROUTINE evolv1(kw,mass,mt,r,lum,mc,rc,menv,renv,ospin, dms = 0.01d0*mt endif * -* Calculate the rate of angular momentum loss due to magnetic braking -* and/or mass loss. +* Evaluate convective/radiative limits for a variety of stars as based +* on the work of Belczynski et al. (2008). As option. +* Note only certain stars of type k = 0, 1, 2, 4, 5, 6, 9 **double check +* 0, 1 have range of 0.35-Mms,conv. Mms,conv is function of metalicity. +* 2 & 4 have a temperature dependence. Convective if Teff < 10**3.73 +* 3, 5 & 6 are giants and have no radiative envelopes. +* 9's envelope is convective when M < Mhe,conv, Mhe,conv = 3.0Msun. +* + if(kw.le.1)then +* Main sequence mass limit varying with metallicity + convradcomp = 1.25d0 + if(z.gt.0.001d0.and.z.lt.zsun)then + convradcomp = 0.747d0 + 55.73d0*z - 1532*z*z + elseif(z.le.0.001d0)then + convradcomp = 0.8d0 + endif + elseif(kw.eq.2.or.kw.eq.4)then +* H-rich HG and CHeB temperature limit + convradcomp = 10.d0**3.73d0 + elseif(kw.le.6)then +* No limit or all other giant stars +* (compare this large value to mass). + convradcomp = 99999999.d0 + elseif(kw.eq.9)then +* Mass limit for evolved He star + convradcomp = 3.d0 + endif +* +*Magnetic Braking * - if(j.gt.1)then - djt = (2.d0/3.d0)*(dms/(1.0d+06*dtm))*r*r*ospin + djmb = 0.d0 + if(htpmb.eq.0)then if(mt.gt.0.35d0.and.kw.lt.10)then djmb = 5.83d-16*menv*(r*ospin)**3/mt - djt = djt + djmb + djspint = djspint + djmb + endif + else + if((menv.gt.0.d0.and. + & ((kw.le.1.and.mt.gt.0.35d0.and. + & mt.le.convradcomp).or. + & (kw.eq.2.and.teff.le.convradcomp).or. + & (kw.eq.4.and.teff.le.convradcomp).or. + & ((kw.eq.3).or.(kw.eq.5).or. + & (kw.eq.6)))))then +* MB given in Ivanova & Taam (2002) + if(ospin.le.wx) djmb = kw3 * r**4.0d0 * + & (ospin/wsun)**3.0d0 + if(ospin.gt.wx) djmb = kw3 * r**4.0d0 * + & (ospin/wsun)**1.3d0 * + & (wx/wsun)**1.7d0 + djspint = djspint + djmb + endif + endif + if(djmb.gt.tiny)then + dtj = 0.03d0*jspin/ABS(djmb) + dt = MIN(dt,dtj) + endif + if(kw.eq.13.and.pulsar.gt.0)then +* +* NS magnetic braking. PK. +* Single NS evolution. +* + if(bdecayfac.eq.0)then + if(B_0.eq.0.d0)then + B = 0.d0 + elseif((tphys-epoch-tacc).lt.tiny)then + B = B_0*EXP(-CK*bacc) + Bbot + else + B = B_0*EXP(-(tphys-epoch-tacc)/bconst)* + & EXP(-CK*bacc) + Bbot + endif + else + if(B_0.eq.0.d0)then + B = 0.d0 + elseif((tphys-epoch-tacc).lt.tiny.and. + & bacc.eq.0.d0)then + B = B_0 + Bbot + elseif((tphys-epoch-tacc).lt.tiny)then + B = B_0/(1.d0 + (bacc/1.0d-6)) + Bbot + elseif(bacc.eq.0.d0)then + B = B_0*EXP(-(tphys-epoch-tacc)/bconst) + & + Bbot + else + B = B_0*EXP(-(tphys-epoch-tacc)/bconst)/ + & (1.d0 + (bacc/1.0d-6)) + Bbot + endif + endif + omdot = Kconst*B*B*ospin**3 + djmb = 0.4d0*mass*r*r*omdot + djspint = djspint + djmb +* Consider update of time-stepping due to dj, i.e. dt = dj/(dj/dt). + if(djmb.gt.tiny)then + dtj = 0.1d0*(jspin/ABS(djmb)) + dt = MIN(dt,dtj) + endif + endif + dtm = dt/1.0d+06 + +* +* Update mass and intrinsic spin (checking that the star is not spun +* past the equilibrium) and reset epoch for a MS (and possibly a HG) star. +* +* + if(eqspin.gt.0.d0.and.ABS(dspint).gt.tiny)then + if(intpol.eq.0)then + if(dspint.ge.0.d0)then + dspint = MIN(dspint,(eqspin-ospin)/dt) + else + dspint = MAX(dspint,(eqspin-ospin)/dt) + endif + djt = (k2str*(mt-mc)*r*r + + & k3*mc*rc*rc)*dspint + djspint = djspint - djt + endif + endif +* + jspin = MAX(1.0d-10,jspin - djspint*dt) +* +* Ensure that the star does not spin up beyond break-up. +* + ospbru = twopi*SQRT(mt*aursun**3/r**3) + jspbru = (k2str*(mt-mc)*r*r + + & k3*mc*rc*rc)*ospbru + if((jspin.gt.jspbru.or. + & (jspin.eq.1.0d-10.and.intpol.gt.0)).and. + & ABS(dtm).gt.tiny)then !PDK add check for jspin intpol issues here. +* If rapidly spinning star (generally NS) pushed over the edge (spins up so +* much so that jspin would have become negative according to bse, so it sets +* it to 1.d0-10) in intpol then give it the stars maximum spin. + mew = 1.d0 + if(djtx.gt.0.d0.and.jspin.gt.1.d0-10)then !dont go in here if jspin hit wall. + mew = MIN(mew,(jspin - jspbru)/djtx) endif + jspin = jspbru endif + * -* Update mass and time and reset epoch for a MS (and possibly a HG) star. +* Update mass * - if(dms.gt.0.d0)then + if(ABS(dms).gt.tiny)then mt = mt - dms if(kw.le.2.or.kw.eq.7)then m0 = mass @@ -146,9 +291,9 @@ SUBROUTINE evolv1(kw,mass,mt,r,lum,mc,rc,menv,renv,ospin, mass = mt tmold = tm tbgold = tscls(1) - CALL star(kw,mass,mt,tm,tn,tscls,lums,GB,zpars) + CALL star(kw,m0,mass,tm,tn,tscls,lums,GB,zpars) if(kw.eq.2)then - if(GB(9).lt.mc1.or.m0.gt.zpars(3))then + if(GB(9).lt.mc.or.m0.gt.zpars(3))then mass = m0 else epoch = tm + (tscls(1) - tm)*(ajhold-tmold)/ @@ -159,6 +304,20 @@ SUBROUTINE evolv1(kw,mass,mt,r,lum,mc,rc,menv,renv,ospin, epoch = tphys - ajhold*tm/tmold endif endif +* Update NS magnetic field owing to accretion, as a function of mass accreted. PK. + if(kw.eq.13.and.pulsar.gt.0)then + if(dms.lt.0.d0)then !negative dms is mass gained. +* When propeller ev. include .not.prop here... + b_mdot = dms/dt + if(b_mdot_lim.gt.0.d0.and.b_mdot.gt.b_mdot_lim)then + bacc = bacc - dms + tacc = tacc + dtm + elseif(b_mdot_lim.le.0.d0)then + bacc = bacc - dms + tacc = tacc + dtm + endif + endif + endif endif tphys2 = tphys tphys = tphys + dtm @@ -166,20 +325,21 @@ SUBROUTINE evolv1(kw,mass,mt,r,lum,mc,rc,menv,renv,ospin, * Find the landmark luminosities and timescales as well as setting * the GB parameters. * - write(*,*) aj,tphys,epoch + aj = tphys - epoch CALL star(kw,mass,mt,tm,tn,tscls,lums,GB,zpars) * * Find the current radius, luminosity, core mass and stellar type * given the initial mass, current mass, metallicity and age kwold = kw + bhspin = 0.0 CALL hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, - & r,lum,kw,mc,rc,menv,renv,k2,1) + & r,lum,kw,mc,rc,menv,renv,k2,bhspin,1) * * If mass loss has occurred and no type change then check that we * have indeed limited the radius change to 10%. * - if(kw.eq.kwold.and.dms.gt.0.d0.and.rflag.ne.0)then + if(kw.eq.kwold.and.dms.gt.0.d0)then mt2 = mt + dms dml = dms/dtm it = 0 @@ -192,7 +352,7 @@ SUBROUTINE evolv1(kw,mass,mt,r,lum,mc,rc,menv,renv,ospin, CALL FLUSH(99) WRITE(*,*)' STOP: EVOLV1 FATAL ERROR ' CALL exit(0) - STOP + STOP endif dtdr = dtm/ABS(dr) dtm = alpha2*MAX(rm0,r)*dtdr @@ -230,7 +390,6 @@ SUBROUTINE evolv1(kw,mass,mt,r,lum,mc,rc,menv,renv,ospin, * if(j.eq.1)then if(tphys.lt.tiny.and.ospin.lt.0.001d0)then -* write (*,*) vrotf(mt,0),mt,r ospin = 45.35d0*vrotf(mt,0)/r endif jspin = ospin*(k2*r*r*(mt-mc)+k3*rc*rc*mc) @@ -242,59 +401,114 @@ SUBROUTINE evolv1(kw,mass,mt,r,lum,mc,rc,menv,renv,ospin, * Test for changes in evolution type. * if(j.eq.1.or.kw.ne.kwold)then + if((kwold.le.12.and.(kw.eq.13.or.kw.eq.14)))then + if(kw.eq.13.and.ecsn.gt.0.d0)then + if(kwold.le.6)then + if(mass.le.zpars(5))then + if(sigma.gt.0.d0.and.sigmadiv.gt.0.d0)then + sigma = -sigmahold/sigmadiv + elseif(sigma.gt.0.d0.and.sigmadiv.lt.0.d0)then + sigma = sigmadiv + endif + endif + elseif(kwold.ge.7.and.kwold.le.9)then + if(mt.gt.ecsn_mlow.and.mt.le.ecsn)then +* BSE orgi: 1.6-2.25, Pod: 1.4-2.5, StarTrack: 1.83-2.25 (all in Msun) + if(sigma.gt.0.d0.and.sigmadiv.gt.0.d0)then + sigma = -sigmahold/sigmadiv + elseif(sigma.gt.0.d0.and.sigmadiv.lt.0.d0)then + sigma = sigmadiv + endif + endif + elseif(kwold.ge.10.or.kwold.eq.12)then +* AIC formation, will never happen here but... + if(sigma.gt.0.d0.and.sigmadiv.gt.0.d0)then + sigma = -sigmahold/sigmadiv + elseif(sigma.gt.0.d0.and.sigmadiv.lt.0.d0)then + sigma = sigmadiv + endif + endif + elseif(kw.eq.13.and.aic.gt.0)then + if(kwold.ge.10.or.kwold.eq.12)then +* AIC formation, will never happen here but... + if(sigma.gt.0.d0.and.sigmadiv.gt.0.d0)then + sigma = -sigmahold/sigmadiv + elseif(sigma.gt.0.d0.and.sigmadiv.lt.0.d0)then + sigma = sigmadiv + endif + endif + endif + evolve_type = 14.d0 + FLOAT(k) + teff = 1000.d0*((1130.d0*lum/ + & (r**2.d0))**(1.d0/4.d0)) + if(B_0.eq.0.d0)then !PK. + b01_bcm = 0.d0 + elseif(B_0.gt.0.d0.and.B.eq.0.d0)then + b01_bcm = B_0 + else + b01_bcm = B + endif + CALL kick(kw,mass,mt,0.d0,0.d0,-1.d0,0.d0,vk,k, + & 0.d0,fallback,sigmahold,kick_info, + & disrupt,bkick) + endif * -* Force new NS or BH to have a one second period. -* - if((kw.eq.13.or.kw.eq.14).and.kw.ne.kwold)then - ospin = 2.0d+08 +* Force new NS or BH to have a birth spin peirod and magnetic field. +* + if(kw.eq.13.or.kw.eq.14)then + if(tphys-epoch.lt.tiny)then + if(kw.eq.13)then +* + 170 u1 = ran3(idum1) + if(u1.ge.1.d0) goto 170 + u2 = ran3(idum1) + s = sqrt(-2.d0*LOG(1.d0-u1))*COS(twopi*u2) + s = 0.7d0*s - 0.6d0 + if(s.ge.0.013d0.or.s.le.-1.5d0) goto 170 + ospin = (twopi*yearsc)/(10.d0**s) + 174 u1 = ran3(idum1) + if(u1.ge.1.d0) goto 174 + u2 = ran3(idum1) + s = sqrt(-2.d0*LOG(1.d0-u1))*COS(twopi*u2) + s = 0.68d0*s + 12.6d0 + if(s.lt.11.5d0.or.s.gt.13.8d0) goto 174 + B_0 = 10.d0**s + bacc = 0.d0 ! If it has been a NS before reset + tacc = 0.d0 + endif + else + ospin = 2.0d+08 + endif + bacc = 0.d0 !make sure if its been a NS before its now a new one... + tacc = 0.d0 jspin = k3*rc*rc*mc*ospin - CALL kick(kw,mass,mc,mt,0.d0,0.d0,-1.d0,0.d0,vk,1, - & 0.d0,fallback,kick_info) - endif - jp = jp + 1 - spp(jp,1) = tphys - spp(jp,2) = float(kw) - if(kw.eq.15)then - spp(jp,3) = mass - goto 90 - else - spp(jp,3) = mt + sigma = sigmahold !reset sigma after possible ECSN kick dist. endif - endif + * * Record values for plotting and reset epoch. * - epoch = tphys - aj - if((isave.and.tphys.ge.tsave).or.iplot)then + epoch = tphys - aj + + CALL WRITESPP(jp,tphys,evolve_type, + & mass,kw,aj,tm,mc,r,m0,lum, + & teff,rc,menv,renv,ospin,b01_bcm, + & bacc,tacc,epoch,bhspin) + CALL WRITESCM(ip,tphys,kw,m0,mass,lum,r,teff, + & mc,rc,menv,renv,epoch,deltam, + & ospin,b01_bcm,SN) + jp = jp + 1 ip = ip + 1 - scm(ip,1) = tphys - scm(ip,2) = float(kw) - scm(ip,3) = mass - scm(ip,4) = mt - scm(ip,5) = log10(lum) - scm(ip,6) = log10(r) - teff = 1000.d0*((1130.d0*lum/(r**2.d0))**(1.d0/4.d0)) - scm(ip,7) = log10(teff) - scm(ip,8) = mc - scm(ip,9) = rc - scm(ip,10) = menv - scm(ip,11) = renv - scm(ip,12) = epoch - scm(ip,13) = ospin - if(isave) tsave = tsave + dtp - if(tphysf.lt.tiny)then - ip = ip + 1 - do 35 , it = 1,13 - scm(ip,it) = scm(ip-1,it) - 35 continue + if(kw.eq.15)then + goto 90 endif endif * if(tphys.ge.tphysf)then - jp = jp + 1 - spp(jp,1) = tphys - spp(jp,2) = float(kw) - spp(jp,3) = mt + CALL WRITESPP(jp,tphys,evolve_type, + & mass,kw,aj,tm,mc,r,m0,lum, + & teff,rc,menv,renv,ospin,B_0, + & bacc,tacc,epoch,bhspin) goto 90 endif * @@ -314,7 +528,7 @@ SUBROUTINE evolv1(kw,mass,mt,r,lum,mc,rc,menv,renv,ospin, * Check final radius for too large a jump. * aj = MAX(aj,aj*(1.d0-eps)+dtr) - mc1 = mc + mc1 = mc CALL hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, & r1,lum1,kw,mc1,rc1,menv1,renv1,k21,1) dr = r1 - rm0 @@ -334,7 +548,7 @@ SUBROUTINE evolv1(kw,mass,mt,r,lum,mc,rc,menv,renv,ospin, * will not change the radius. * 40 aj = ajhold + dtm - mc1 = mc + mc1 = mc CALL hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, & r1,lum1,kw,mc1,rc1,menv1,renv1,k21,1) dr = r1 - rm0 @@ -345,7 +559,7 @@ SUBROUTINE evolv1(kw,mass,mt,r,lum,mc,rc,menv,renv,ospin, CALL FLUSH(99) WRITE(*,*)' STOP: EVOLV1 FATAL ERROR ' CALL exit(0) - STOP + STOP endif if(ABS(dr).gt.0.1d0*rm0)then dtdr = dtm/ABS(dr) @@ -356,8 +570,8 @@ SUBROUTINE evolv1(kw,mass,mt,r,lum,mc,rc,menv,renv,ospin, * 50 continue * -* Ensure that change of type has not occurred during radius check. -* This is rare but may occur for HG stars of ZAMS mass > 50 Msun. +* Ensure that change of type has not occurred during radius check. +* This is rare but may occur for HG stars of ZAMS mass > 50 Msun. * if(kw.ne.kwold)then kw = kwold From af46cc80518300173adb0b5260bb09e36c39c655 Mon Sep 17 00:00:00 2001 From: hskobe Date: Mon, 13 Apr 2026 12:27:49 -0400 Subject: [PATCH 04/15] Finished updates for evolv2.f --- src/cosmic/src/evolv2.f | 78 ++++++++++++++++++++++++----------------- 1 file changed, 46 insertions(+), 32 deletions(-) diff --git a/src/cosmic/src/evolv2.f b/src/cosmic/src/evolv2.f index 9084d5265..600b00c3e 100644 --- a/src/cosmic/src/evolv2.f +++ b/src/cosmic/src/evolv2.f @@ -154,7 +154,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, *** * INTEGER loop,iter,intpol,k,ip,j1,j2 - INTEGER bcm_index_out, bpp_index_out + INTEGER bcm_index_out, bpp_index_out, kstar1, kstar2 INTEGER kcomp1,kcomp2,formation(2) PARAMETER(loop=40000) INTEGER kstar(2),kw,kst,kw1,kw2,kmin,kmax @@ -217,6 +217,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, REAL*8 qc_fixed LOGICAL switchedCE,disrupt + Cf2py intent(in) kstar Cf2py intent(in) mass Cf2py intent(in) tb @@ -246,6 +247,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, Cf2py intent(out) bcm_index_out Cf2py intent(out) kick_info_out + if(using_cmc.eq.0)then CALL instar endif @@ -255,6 +257,11 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, * * CE2flag = 0 + +* +* Get merger type from julia call +* + kstar1_bpp = 0 kstar2_bpp = 0 @@ -405,7 +412,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, if(ospin(1).lt.0.d0) ospin(1) = oorb if(ospin(2).lt.0.d0) ospin(2) = oorb endif -* + do 500 , k = kmin,kmax age = tphys - epoch(k) mc = massc(k) @@ -447,6 +454,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, endif * 500 continue +* * if(output) write(*,*)'Init:',mass(1),mass(2),massc(1),massc(2), & rad(1),rad(2),kstar(1),kstar(2),sep,ospin(1),ospin(2),jspin(1), @@ -1574,6 +1582,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, rol(k) = 10000.d0*rad(k) 508 continue endif + * if((tphys.lt.tiny.and.ABS(dtm).lt.tiny.and. & (mass2i.lt.0.1d0.or..not.sgl)).or.snova)then @@ -2379,6 +2388,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, endif coel = .true. binstate = 1 + if(mass(j2).gt.0.d0)then mass(j1) = 0.d0 kstar(j1) = 15 @@ -2774,6 +2784,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, endif coel = .true. binstate = 1 + goto 135 elseif(kstar(j1).eq.13)then * @@ -2788,11 +2799,13 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, kstar(j2) = 14 coel = .true. binstate = 1 + goto 135 elseif(kstar(j1).eq.14)then * * Both stars are black holes. Let them merge quietly. * + CALL CONCATKSTARS(kstar(j1), kstar(j2), mergertype) dm1 = mass(j1) mass(j1) = 0.d0 @@ -3015,12 +3028,12 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, if(kstar(j2).le.2.or.kstar(j2).eq.4)then - if(acc_lim.eq.-1.or.acc_lim.eq.-3)then + if(acc_lim(j2).eq.-1.or.acc_lim(j2).eq.-3)then dm2 = MIN(1.d0,10.d0*taum/tkh(j2))*dm1 - elseif(acc_lim.eq.-2.or.acc_lim.eq.-4)then + elseif(acc_lim(j2).eq.-2.or.acc_lim(j2).eq.-4)then dm2 = MIN(1.d0,taum/tkh(j2))*dm1 - elseif(acc_lim.ge.0.d0)then - dm2 = acc_lim*dm1 + elseif(acc_lim(j2).ge.0.d0)then + dm2 = acc_lim(j2)*dm1 endif elseif(kstar(j2).ge.7.and.kstar(j2).le.9)then * @@ -3028,18 +3041,18 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, * or SAGB star unless the primary is also a helium star. * if(kstar(j1).ge.7)then - if(acc_lim.eq.-1.or.acc_lim.eq.-3)then + if(acc_lim(j2).eq.-1.or.acc_lim(j2).eq.-3)then dm2 = MIN(1.d0,10.d0*taum/tkh(j2))*dm1 - elseif(acc_lim.eq.-2.or.acc_lim.eq.-4)then + elseif(acc_lim(j2).eq.-2.or.acc_lim(j2).eq.-4)then dm2 = MIN(1.d0,taum/tkh(j2))*dm1 - elseif(acc_lim.ge.0.d0)then - dm2 = acc_lim*dm1 + elseif(acc_lim(j2).ge.0.d0)then + dm2 = acc_lim(j2)*dm1 endif else - if(acc_lim.lt.0.d0)then + if(acc_lim(j2).lt.0.d0)then dm2 = dm1 - elseif(acc_lim.ge.0.d0)then - dm2 = acc_lim*dm1 + elseif(acc_lim(j2).ge.0.d0)then + dm2 = acc_lim(j2)*dm1 endif dmchk = dm2 - 1.05d0*dms(j2) if(dmchk.gt.0.d0.and.dm2/mass(j2).gt.1.0d-04)then @@ -3066,32 +3079,32 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, * Accrete until a nova explosion blows away most of the accreted material. * novae = .true. - if(acc_lim.lt.0.d0)then + if(acc_lim(j2).lt.0.d0)then dm2 = MIN(dm1,dme) if(dm2.lt.dm1) supedd = .true. - elseif(acc_lim.ge.0.d0)then - dm2 = MIN(dm2,acc_lim*dm1) - if(dm2.lt.acc_lim*dm1) supedd = .true. + elseif(acc_lim(j2).ge.0.d0)then + dm2 = MIN(dm2,acc_lim(j2)*dm1) + if(dm2.lt.acc_lim(j2)*dm1) supedd = .true. endif dm22 = epsnov*dm2 else * * Steady burning at the surface * - if(acc_lim.lt.0.d0)then + if(acc_lim(j2).lt.0.d0)then dm2 = dm1 - elseif(acc_lim.ge.0.d0)then - dm2 = acc_lim*dm1 + elseif(acc_lim(j2).ge.0.d0)then + dm2 = acc_lim(j2)*dm1 endif endif else * * Make a new giant envelope. * - if(acc_lim.lt.0.d0)then + if(acc_lim(j2).lt.0.d0)then dm2 = dm1 - elseif(acc_lim.ge.0.d0)then - dm2 = MIN(dm2,acc_lim*dm1) + elseif(acc_lim(j2).ge.0.d0)then + dm2 = MIN(dm2,acc_lim(j2)*dm1) endif * * Check for planets or low-mass WDs. @@ -3111,14 +3124,14 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, elseif(kstar(j2).eq.3.or.kstar(j2).eq.5.or.kstar(j2).eq.6)then * We have a giant w/ kstar(j2) = 3,5,6 * - if(acc_lim.eq.-1.or.acc_lim.eq.-2)then + if(acc_lim(j2).eq.-1.or.acc_lim(j2).eq.-2)then dm2 = dm1 - elseif(acc_lim.eq.-3)then + elseif(acc_lim(j2).eq.-3)then dm2 = MIN(1.d0,10*taum/tkh(j2))*dm1 - elseif(acc_lim.eq.-4)then + elseif(acc_lim(j2).eq.-4)then dm2 = MIN(1.d0,taum/tkh(j2))*dm1 - elseif(acc_lim.ge.0.d0)then - dm2 = MIN(dm2,acc_lim*dm1) + elseif(acc_lim(j2).ge.0.d0)then + dm2 = MIN(dm2,acc_lim(j2)*dm1) endif endif @@ -3127,7 +3140,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, * Impose the Eddington limit. * if(kstar(j2).ge.10)then - if(acc_lim.lt.0.d0)then + if(acc_lim(j2).lt.0.d0)then * * If there is wind accretion the total amount of mass change is * dms(j2) = dmr(j2) - dmt(j2), where dmt(j2) is the accretion @@ -3141,20 +3154,20 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, * if(supedd.eqv..true.) dm2 = 0.d0 if(dm2.lt.dm1) supedd = .true. - elseif(acc_lim.ge.0.d0)then + elseif(acc_lim(j2).ge.0.d0)then * * If there is wind accretion the total amount of mass change is * dms(j2) = dmr(j2) - dmt(j2), where dmt(j2) is the accretion * from the companion. We should limit to the Eddington limit minus * the amount of accretion that is already coming in from Winds * - dm2 = MIN(acc_lim*dm1,dme + dms(j2)) + dm2 = MIN(acc_lim(j2)*dm1,dme + dms(j2)) * * If we already hit supereddington wind accretion, don't add * any more mass through RLO * if(supedd.eqv..true.) dm2 = 0.d0 - if(dm2.lt.acc_lim*dm1) supedd = .true. + if(dm2.lt.acc_lim(j2)*dm1) supedd = .true. endif * @@ -4389,6 +4402,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf, 135 continue * sgl = .true. + if(kstar(1).eq.13.and.mergemsp.eq.1.and. & notamerger.eq.0)then s = (twopi*yearsc)/ospin(1) From 2cc01a2d5068c64fbc48b6d1ecf46ff2d48b55b4 Mon Sep 17 00:00:00 2001 From: hskobe Date: Mon, 13 Apr 2026 12:34:07 -0400 Subject: [PATCH 05/15] Added alpha1/acclim updates to evolve.py --- src/cosmic/evolve.py | 28 +++++++++++++++++++++++++++- 1 file changed, 27 insertions(+), 1 deletion(-) diff --git a/src/cosmic/evolve.py b/src/cosmic/evolve.py index 68a1e6b0d..5e5a9b0e0 100644 --- a/src/cosmic/evolve.py +++ b/src/cosmic/evolve.py @@ -123,7 +123,7 @@ else: INITIAL_BINARY_TABLE_SAVE_COLUMNS = INITIAL_CONDITIONS_PASS_COLUMNS.copy() -for col in ['natal_kick_array', 'qcrit_array', 'fprimc_array']: +for col in ['natal_kick_array', 'qcrit_array', 'fprimc_array', 'alpha1', 'acc_lim']: INITIAL_BINARY_TABLE_SAVE_COLUMNS.remove(col) NATAL_KICK_COLUMNS = ['natal_kick', @@ -139,10 +139,14 @@ QCRIT_COLUMNS = ['qcrit_{0}'.format(kstar) for kstar in range(0, 16)] FPRIMC_COLUMNS = ['fprimc_{0}'.format(kstar) for kstar in range(0, 16)] +ALPHA_COLUMNS = ['alpha1_{0}'.format(star) for star in range(0, 2)] +ACCLIM_COLUMNS = ['acc_lim_{0}'.format(star) for star in range(0, 2)] INITIAL_BINARY_TABLE_SAVE_COLUMNS.extend(FLATTENED_NATAL_KICK_COLUMNS) INITIAL_BINARY_TABLE_SAVE_COLUMNS.extend(QCRIT_COLUMNS) INITIAL_BINARY_TABLE_SAVE_COLUMNS.extend(FPRIMC_COLUMNS) +INITIAL_BINARY_TABLE_SAVE_COLUMNS.extend(ALPHA_COLUMNS) +INITIAL_BINARY_TABLE_SAVE_COLUMNS.extend(ACCLIM_COLUMNS) # BSE doesn't need the binary fraction, so just add to columns for saving INITIAL_BINARY_TABLE_SAVE_COLUMNS.insert(7, 'binfrac') @@ -343,6 +347,22 @@ def evolve(self, initialbinarytable, pool=None, bpp_columns=None, bcm_columns=No initialbinarytable[col] = BSEDict['fprimc_array'][kstar] else: new_cols[col] = BSEDict['fprimc_array'][kstar] + elif k == 'alpha1': + columns_values = [BSEDict['alpha1']] * len(initialbinarytable) + initialbinarytable = initialbinarytable.assign(alpha1=columns_values) + for kstar in range(0,2): + columns_values = pd.Series([BSEDict['alpha1'][kstar]] * len(initialbinarytable), + index=initialbinarytable.index, + name='alpha1_{0}'.format(kstar)) + initialbinarytable.loc[:, 'alpha1_{0}'.format(kstar)] = columns_values + elif k == 'acc_lim': + columns_values = [BSEDict['acc_lim']] * len(initialbinarytable) + initialbinarytable = initialbinarytable.assign(acc_lim=columns_values) + for kstar in range(0,2): + columns_values = pd.Series([BSEDict['acc_lim'][kstar]] * len(initialbinarytable), + index=initialbinarytable.index, + name='acc_lim_{0}'.format(kstar)) + initialbinarytable.loc[:, 'acc_lim_{0}'.format(kstar)] = columns_values else: # base case: if it's present, overwrite, if not, add to a list of new columns (see below) if k in initialbinarytable.columns: @@ -386,6 +406,12 @@ def evolve(self, initialbinarytable, pool=None, bpp_columns=None, bcm_columns=No if (pd.Series(FPRIMC_COLUMNS).isin(initialbinarytable.keys()).all()) and ('fprimc_array' not in BSEDict): initialbinarytable = initialbinarytable.assign(fprimc_array=initialbinarytable[FPRIMC_COLUMNS].values.tolist()) + if (pd.Series(ALPHA_COLUMNS).isin(initialbinarytable.keys()).all()) and ('alpha1' not in BSEDict): + initialbinarytable = initialbinarytable.assign(alpha1=initialbinarytable[ALPHA_COLUMNS].values.tolist()) + + if (pd.Series(ACCLIM_COLUMNS).isin(initialbinarytable.keys()).all()) and ('acc_lim' not in BSEDict): + initialbinarytable = initialbinarytable.assign(acc_lim=initialbinarytable[ACCLIM_COLUMNS].values.tolist()) + # update timesteps based on mass modifier mass_modifier_cols = set(['pts1', 'pts2', 'pts3']).difference(columns_in_passed_initC) if dt_mass_modifiers and len(mass_modifier_cols) != 0: From 1c5e177f25de26cc93c441e52df17fe3f6b33b8a Mon Sep 17 00:00:00 2001 From: hskobe Date: Mon, 13 Apr 2026 16:49:37 -0400 Subject: [PATCH 06/15] Updated utils.py with alpha1/acc_lim raised ValueError --- src/cosmic/utils.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/cosmic/utils.py b/src/cosmic/utils.py index 3aaf10ed2..7d1ce1c6a 100644 --- a/src/cosmic/utils.py +++ b/src/cosmic/utils.py @@ -1266,6 +1266,20 @@ def error_check(BSEDict, filters=None, convergence=None, sampling=None): f"qcrit_array values must be >= 0 and there must be 16 values " f'(you set them to {BSEDict["qcrit_array"]}], length={len(BSEDict["qcrit_array"])})' ) + + if "alpah1" in BSEDict.keys(): + if np.any(np.array(BSEDict["alpha1"]) < 0.0): + raise ValueError( + f"alpha1 values must be >= 0" + f'(you set them to {BSEDict["alpha1"]}])' + ) + + if "acc_lim" in BSEDict.keys(): + if np.any(np.array(BSEDict["acc_lim"]) not in [-1, -2, -3, -4]) or np.any(np.array(BSEDict["acc_lim"]) < 0.0): + raise ValueError( + f"acc_lim values must be set to -1, -2, -3, -4 or be >=0 " + f'(you set them to {BSEDict["acc_lim"]}])' + ) return From 1bc18bed0ec7044690027a9c246ad3e617883556 Mon Sep 17 00:00:00 2001 From: hskobe Date: Tue, 14 Apr 2026 17:48:51 -0400 Subject: [PATCH 07/15] Updated to have multi alpha1/acc_lim --- examples/Params.ini | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/Params.ini b/examples/Params.ini index 25541a3ea..751706a78 100644 --- a/examples/Params.ini +++ b/examples/Params.ini @@ -357,7 +357,7 @@ LBV_flag = 1 ; positive values - Sets the common-envelope efficiency parameter, \(\alpha\) ; 1.0 - Default value ; Default: 1.0 -alpha1 = 1.0 +alpha1 = [1.0, 1.0] ; lambdaf ; Binding energy factor for common envelope evolution. The initial binding energy of the stellar envelope is proportional to \(1 / \lambda\). See Hurley+2002, Eq. 69. @@ -451,7 +451,7 @@ qcrit_array = [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] ; -5 - As 5, 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 ; -6 - As 6, 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 ; Default: 5 -kickflag = 5 +kickflag = 1 ; sigma ; Sets the dispersion in the Maxwellian for the SN kick velocity in km/s @@ -774,7 +774,7 @@ don_lim = -1 ; positive values - sets overall fraction of donor material that is accreted, with the rest being lost from the system (``acc_lim = 0.5`` assumes 50% accretion efficiency as in Belczynski+2008) ; 0.0 - assume all mass is lost from the system (fully nonconservative) ; Default: -1 -acc_lim = -1 +acc_lim = [-1, -1] ;;;;;;;;;;;;; ;;; Tides ;;; From 5760eb1b2e579b9877683595a5c675d6ccb97d7d Mon Sep 17 00:00:00 2001 From: hskobe Date: Tue, 14 Apr 2026 17:50:28 -0400 Subject: [PATCH 08/15] Fixed issues with acc_lim raised error; test passed when generating/evolving population --- src/cosmic/utils.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/src/cosmic/utils.py b/src/cosmic/utils.py index 7d1ce1c6a..978519120 100644 --- a/src/cosmic/utils.py +++ b/src/cosmic/utils.py @@ -1146,7 +1146,7 @@ def error_check(BSEDict, filters=None, convergence=None, sampling=None): settings_path = io_files("cosmic.data").joinpath('cosmic-settings.json') settings = json.loads(settings_path.read_text(encoding='utf-8')) - handle_separately = ['qcrit_array', 'natal_kick_array', 'fprimc_array'] + handle_separately = ['qcrit_array', 'natal_kick_array', 'fprimc_array', 'alpha1', 'acc_lim'] # go through the different categories in the settings file for cat in settings: @@ -1267,7 +1267,7 @@ def error_check(BSEDict, filters=None, convergence=None, sampling=None): f'(you set them to {BSEDict["qcrit_array"]}], length={len(BSEDict["qcrit_array"])})' ) - if "alpah1" in BSEDict.keys(): + if "alpha1" in BSEDict.keys(): if np.any(np.array(BSEDict["alpha1"]) < 0.0): raise ValueError( f"alpha1 values must be >= 0" @@ -1275,11 +1275,13 @@ def error_check(BSEDict, filters=None, convergence=None, sampling=None): ) if "acc_lim" in BSEDict.keys(): - if np.any(np.array(BSEDict["acc_lim"]) not in [-1, -2, -3, -4]) or np.any(np.array(BSEDict["acc_lim"]) < 0.0): - raise ValueError( - f"acc_lim values must be set to -1, -2, -3, -4 or be >=0 " - f'(you set them to {BSEDict["acc_lim"]}])' - ) + for f in BSEDict["acc_lim"]: + if f not in [-1, -2, -3, -4]: + if f < 0.0: + raise ValueError( + f"acc_lim values must be set to -1, -2, -3, -4 or be >=0 " + f'(you set them to {BSEDict["acc_lim"]})' + ) return From 4d17f40fc139b9583f7190d5361d2db4d81d62ad Mon Sep 17 00:00:00 2001 From: hskobe Date: Tue, 14 Apr 2026 18:04:12 -0400 Subject: [PATCH 09/15] Fixed typo in raised ValueError message for alpha --- src/cosmic/utils.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/cosmic/utils.py b/src/cosmic/utils.py index 978519120..50be36f70 100644 --- a/src/cosmic/utils.py +++ b/src/cosmic/utils.py @@ -1270,8 +1270,8 @@ def error_check(BSEDict, filters=None, convergence=None, sampling=None): if "alpha1" in BSEDict.keys(): if np.any(np.array(BSEDict["alpha1"]) < 0.0): raise ValueError( - f"alpha1 values must be >= 0" - f'(you set them to {BSEDict["alpha1"]}])' + f"alpha1 values must be >= 0 " + f'(you set them to {BSEDict["alpha1"]})' ) if "acc_lim" in BSEDict.keys(): From 113c876ecc9472e8570d55c81c5ca79ed6844bfa Mon Sep 17 00:00:00 2001 From: hskobe Date: Thu, 16 Apr 2026 11:03:35 -0400 Subject: [PATCH 10/15] Fixed alpha1 index --- 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 5e5a9b0e0..b1ccc0fa9 100644 --- a/src/cosmic/evolve.py +++ b/src/cosmic/evolve.py @@ -139,8 +139,8 @@ QCRIT_COLUMNS = ['qcrit_{0}'.format(kstar) for kstar in range(0, 16)] FPRIMC_COLUMNS = ['fprimc_{0}'.format(kstar) for kstar in range(0, 16)] -ALPHA_COLUMNS = ['alpha1_{0}'.format(star) for star in range(0, 2)] -ACCLIM_COLUMNS = ['acc_lim_{0}'.format(star) for star in range(0, 2)] +ALPHA_COLUMNS = ['alpha1_{0}'.format(star) for star in range(1, 3)] +ACCLIM_COLUMNS = ['acc_lim_{0}'.format(star) for star in range(1, 3)] INITIAL_BINARY_TABLE_SAVE_COLUMNS.extend(FLATTENED_NATAL_KICK_COLUMNS) INITIAL_BINARY_TABLE_SAVE_COLUMNS.extend(QCRIT_COLUMNS) From 9056a650d7e6703e53db3380d63d631463b22a0e Mon Sep 17 00:00:00 2001 From: hskobe Date: Thu, 16 Apr 2026 11:04:15 -0400 Subject: [PATCH 11/15] Added columns for split alpha1/acclim' --- src/cosmic/tests/test_evolve.py | 7 ++++++- src/cosmic/tests/test_kick.py | 7 ++++++- 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/src/cosmic/tests/test_evolve.py b/src/cosmic/tests/test_evolve.py index 11acb19de..b074f85d8 100644 --- a/src/cosmic/tests/test_evolve.py +++ b/src/cosmic/tests/test_evolve.py @@ -26,6 +26,10 @@ INIT_CONDITIONS_NO_BSE_COLUMNS = INIT_CONDITIONS[init_conds_columns] BPP_DF = pd.read_hdf(os.path.join(TEST_DATA_DIR, 'unit_tests_results.hdf5'), key='bpp') BCM_DF = pd.read_hdf(os.path.join(TEST_DATA_DIR, 'unit_tests_results.hdf5'), key='bcm') +INIT_CONDITIONS['alpha1_1'] = 1.0 +INIT_CONDITIONS['alpha1_2'] = 1.0 +INIT_CONDITIONS['acc_lim_1'] = -1 +INIT_CONDITIONS['acc_lim_2'] = -1 BSEFlag_columns = list(set(evolve.INITIAL_BINARY_TABLE_SAVE_COLUMNS) - set(initialbinarytable.INITIAL_CONDITIONS_COLUMNS_ALL)) BSEDict = INIT_CONDITIONS[BSEFlag_columns].to_dict(orient='index')[0] BSEDict['qcrit_array'] = [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] @@ -35,9 +39,10 @@ 2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0] BSEDict['grflag'] = 1 BSEDict['don_lim'] = -1 -BSEDict['acc_lim'] = -1 +BSEDict['acc_lim'] = [-1, -1] BSEDict['wd_mass_lim'] = 0 BSEDict['kick_flag'] = -1 +BSEDict['alpha1'] = [1.0, 1.0] class TestEvolve(unittest.TestCase): """`TestCase` for the cosmic diff --git a/src/cosmic/tests/test_kick.py b/src/cosmic/tests/test_kick.py index 66c1c69da..db6d1a919 100644 --- a/src/cosmic/tests/test_kick.py +++ b/src/cosmic/tests/test_kick.py @@ -18,6 +18,10 @@ TEST_DATA_DIR = os.path.join(os.path.split(__file__)[0], 'data') INIT_CONDITIONS = pd.read_hdf(os.path.join(TEST_DATA_DIR, 'initial_conditions_for_testing.hdf5'), key='initC') +INIT_CONDITIONS['alpha1_1'] = 1.0 +INIT_CONDITIONS['alpha1_2'] = 1.0 +INIT_CONDITIONS['acc_lim_1'] = -1 +INIT_CONDITIONS['acc_lim_2'] = -1 BSEFlag_columns = list(set(INITIAL_BINARY_TABLE_SAVE_COLUMNS) - set(INITIAL_CONDITIONS_COLUMNS_ALL)) BSEDict = INIT_CONDITIONS[BSEFlag_columns].to_dict(orient='index')[0] BSEDict['qcrit_array'] = [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] @@ -27,8 +31,9 @@ 2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0] BSEDict['grflag'] = 1 BSEDict['don_lim'] = -1 -BSEDict['acc_lim'] = -1 +BSEDict['acc_lim'] = [-1, -1] BSEDict['wd_mass_lim'] = 0 +BSEDict['alpha1'] = [1.0, 1.0] # avoid using the same randomseed del BSEDict["bin_num"], BSEDict["randomseed"] From 4a19ac7466ff39d1b63d8589117c019dd364a6c9 Mon Sep 17 00:00:00 2001 From: hskobe Date: Thu, 16 Apr 2026 11:19:35 -0400 Subject: [PATCH 12/15] Updated split alpha/acc_lim for docs --- src/cosmic/data/cosmic-settings.json | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/cosmic/data/cosmic-settings.json b/src/cosmic/data/cosmic-settings.json index a1b9ca26a..e1111fdf2 100644 --- a/src/cosmic/data/cosmic-settings.json +++ b/src/cosmic/data/cosmic-settings.json @@ -725,10 +725,10 @@ "options": [ { "name": "positive values", - "description": "Sets the common-envelope efficiency parameter, \\(\\alpha\\)" + "description": "Sets the common-envelope efficiency parameter, \\(\\alpha\\), for each binary component" }, { - "name": 1.0, + "name": "[1.0, 1.0]", "description": "Default value", "default": true } @@ -1565,14 +1565,13 @@ }, { "name": "acc_lim", - "description": "Limits the amount of mass accreted during Roche-lobe overflow", + "description": "Limits the amount of mass accreted during Roche-lobe overflow where each binary component must have its own limit specified in a list", "type": "number", "options-preface": "", "options": [ { "name": -1, - "description": "limited to 10x the thermal rate of the accretor for MS/HG/CHeB and unlimited for GB/EAGB/AGB stars", - "default": true + "description": "limited to 10x the thermal rate of the accretor for MS/HG/CHeB and unlimited for GB/EAGB/AGB stars" }, { "name": -2, @@ -1593,7 +1592,13 @@ { "name": 0.0, "description": "assume all mass is lost from the system (fully nonconservative)" + }, + { + "name": "[-1, -1]", + "description": "Default choice", + "default": true } + ] }, { From bd8f9b498fa6a5f50ef09cca0794f3145ee7ef40 Mon Sep 17 00:00:00 2001 From: Katie Breivik Date: Thu, 30 Apr 2026 12:54:59 -0400 Subject: [PATCH 13/15] tyring to get tests to pass --- src/cosmic/tests/data/Params.ini | 4 ++-- src/cosmic/tests/test_evolve.py | 9 ++++++++- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/src/cosmic/tests/data/Params.ini b/src/cosmic/tests/data/Params.ini index 5a0f1fbdc..aefc27ac6 100644 --- a/src/cosmic/tests/data/Params.ini +++ b/src/cosmic/tests/data/Params.ini @@ -157,7 +157,7 @@ LBV_flag=2 ; alpha1 is the common-envelope efficiency parameter ; default=1.0 -alpha1 = 1.0 +alpha1 = [1.0, 1.0] ; lambdaf is the binding energy factor for common envelope evolution ; lambdaf>0.0 uses variable lambda prescription written by Onno Pols @@ -387,7 +387,7 @@ don_lim = -1 ; acc_lim = -2: assumes the accretion is limited to tkh_acc for MS/CHeB ; acc_lim = -3: assumes the accretion is limited by 10*tkh_acc for all fusing stars ; acc_lim = -4: assumes the accretion is limited by tkh_for all fusing stars -acc_lim = -1 +acc_lim = [-1, -1] ;;;;;;;;;;;;;;;;;;; ;;; TIDES FLAGS ;;; diff --git a/src/cosmic/tests/test_evolve.py b/src/cosmic/tests/test_evolve.py index b074f85d8..1ed61cc0a 100644 --- a/src/cosmic/tests/test_evolve.py +++ b/src/cosmic/tests/test_evolve.py @@ -26,10 +26,16 @@ INIT_CONDITIONS_NO_BSE_COLUMNS = INIT_CONDITIONS[init_conds_columns] BPP_DF = pd.read_hdf(os.path.join(TEST_DATA_DIR, 'unit_tests_results.hdf5'), key='bpp') BCM_DF = pd.read_hdf(os.path.join(TEST_DATA_DIR, 'unit_tests_results.hdf5'), key='bcm') + INIT_CONDITIONS['alpha1_1'] = 1.0 INIT_CONDITIONS['alpha1_2'] = 1.0 INIT_CONDITIONS['acc_lim_1'] = -1 INIT_CONDITIONS['acc_lim_2'] = -1 +KICK_INITC['alpha1_1'] = 1.0 +KICK_INITC['alpha1_2'] = 1.0 +KICK_INITC['acc_lim_1'] = -1 +KICK_INITC['acc_lim_2'] = -1 + BSEFlag_columns = list(set(evolve.INITIAL_BINARY_TABLE_SAVE_COLUMNS) - set(initialbinarytable.INITIAL_CONDITIONS_COLUMNS_ALL)) BSEDict = INIT_CONDITIONS[BSEFlag_columns].to_dict(orient='index')[0] BSEDict['qcrit_array'] = [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] @@ -44,6 +50,7 @@ BSEDict['kick_flag'] = -1 BSEDict['alpha1'] = [1.0, 1.0] + class TestEvolve(unittest.TestCase): """`TestCase` for the cosmic """ @@ -109,5 +116,5 @@ def test_multi_evolve_with_dict_and_table(self): def test_ejection_velocity_pfahl(self): EvolvedBinaryBPP, EvolvedBinaryBCM, initCond, kick_info = Evolve.evolve( initialbinarytable=KICK_INITC) - + self.assertAlmostEqual(kick_info['vsys_2_total'].iloc[0], 482.346136, places=5) From f2bdb76d6aa95681e7d9a51e63a815962cf88380 Mon Sep 17 00:00:00 2001 From: Katie Breivik Date: Thu, 30 Apr 2026 13:12:01 -0400 Subject: [PATCH 14/15] acutally fixed the tests; needed 0 and 1 not 1 and 2; this is consistent with the metisse-integrate-merge convention --- src/cosmic/evolve.py | 4 ++-- src/cosmic/tests/test_evolve.py | 8 ++++---- src/cosmic/tests/test_kick.py | 4 ++-- src/cosmic/tests/test_utils.py | 2 +- 4 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/cosmic/evolve.py b/src/cosmic/evolve.py index f63ec327a..6a5d70fc0 100644 --- a/src/cosmic/evolve.py +++ b/src/cosmic/evolve.py @@ -140,8 +140,8 @@ QCRIT_COLUMNS = ['qcrit_{0}'.format(kstar) for kstar in range(0, 16)] FPRIMC_COLUMNS = ['fprimc_{0}'.format(kstar) for kstar in range(0, 16)] -ALPHA_COLUMNS = ['alpha1_{0}'.format(star) for star in range(1, 3)] -ACCLIM_COLUMNS = ['acc_lim_{0}'.format(star) for star in range(1, 3)] +ALPHA_COLUMNS = ['alpha1_{0}'.format(star) for star in range(0, 2)] +ACCLIM_COLUMNS = ['acc_lim_{0}'.format(star) for star in range(0, 2)] INITIAL_BINARY_TABLE_SAVE_COLUMNS.extend(FLATTENED_NATAL_KICK_COLUMNS) INITIAL_BINARY_TABLE_SAVE_COLUMNS.extend(QCRIT_COLUMNS) diff --git a/src/cosmic/tests/test_evolve.py b/src/cosmic/tests/test_evolve.py index 1ed61cc0a..c1e593941 100644 --- a/src/cosmic/tests/test_evolve.py +++ b/src/cosmic/tests/test_evolve.py @@ -27,14 +27,14 @@ BPP_DF = pd.read_hdf(os.path.join(TEST_DATA_DIR, 'unit_tests_results.hdf5'), key='bpp') BCM_DF = pd.read_hdf(os.path.join(TEST_DATA_DIR, 'unit_tests_results.hdf5'), key='bcm') +INIT_CONDITIONS['alpha1_0'] = 1.0 INIT_CONDITIONS['alpha1_1'] = 1.0 -INIT_CONDITIONS['alpha1_2'] = 1.0 +INIT_CONDITIONS['acc_lim_0'] = -1 INIT_CONDITIONS['acc_lim_1'] = -1 -INIT_CONDITIONS['acc_lim_2'] = -1 +KICK_INITC['alpha1_0'] = 1.0 KICK_INITC['alpha1_1'] = 1.0 -KICK_INITC['alpha1_2'] = 1.0 +KICK_INITC['acc_lim_0'] = -1 KICK_INITC['acc_lim_1'] = -1 -KICK_INITC['acc_lim_2'] = -1 BSEFlag_columns = list(set(evolve.INITIAL_BINARY_TABLE_SAVE_COLUMNS) - set(initialbinarytable.INITIAL_CONDITIONS_COLUMNS_ALL)) BSEDict = INIT_CONDITIONS[BSEFlag_columns].to_dict(orient='index')[0] diff --git a/src/cosmic/tests/test_kick.py b/src/cosmic/tests/test_kick.py index db6d1a919..1f5314c12 100644 --- a/src/cosmic/tests/test_kick.py +++ b/src/cosmic/tests/test_kick.py @@ -18,10 +18,10 @@ TEST_DATA_DIR = os.path.join(os.path.split(__file__)[0], 'data') INIT_CONDITIONS = pd.read_hdf(os.path.join(TEST_DATA_DIR, 'initial_conditions_for_testing.hdf5'), key='initC') +INIT_CONDITIONS['alpha1_0'] = 1.0 INIT_CONDITIONS['alpha1_1'] = 1.0 -INIT_CONDITIONS['alpha1_2'] = 1.0 +INIT_CONDITIONS['acc_lim_0'] = -1 INIT_CONDITIONS['acc_lim_1'] = -1 -INIT_CONDITIONS['acc_lim_2'] = -1 BSEFlag_columns = list(set(INITIAL_BINARY_TABLE_SAVE_COLUMNS) - set(INITIAL_CONDITIONS_COLUMNS_ALL)) BSEDict = INIT_CONDITIONS[BSEFlag_columns].to_dict(orient='index')[0] BSEDict['qcrit_array'] = [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] diff --git a/src/cosmic/tests/test_utils.py b/src/cosmic/tests/test_utils.py index da8c14192..e740538ae 100644 --- a/src/cosmic/tests/test_utils.py +++ b/src/cosmic/tests/test_utils.py @@ -147,7 +147,7 @@ def test_binwidth_selector(self): self.assertTrue(bw.round(3) == BW_KNUTH) def test_error_check(self): - BSEDict = {'xi': 0.5, 'bhflag': 1, 'neta': 0.5, 'windflag': 3, 'rtmsflag' : 0, 'wdflag': 0, 'alpha1': 1.0, 'pts1': 0.05, 'pts3': 0.02, 'pts2': 0.01, 'epsnov': 0.001, 'hewind': 1.0, 'ck': 1000, 'bwind': 0.0, 'lambdaf': 0.5, 'mxns': 3.0, 'beta': -1.0, 'tflag': 1, 'acc2': 1.5, 'nsflag': 4, 'ceflag': 0, 'eddfac': 1.0, 'ifflag': 0, 'bconst': 3000, 'sigma': 265.0, 'gamma': -2.0, 'pisn': 45.0, 'natal_kick_array' :[[-100.0,-100.0,-100.0,-100.0,0.0], [-100.0,-100.0,-100.0,-100.0,0.0]], 'bhsigmafrac' : 1.0, 'polar_kick_angle' : 90, 'qcrit_array' : [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], 'cekickflag' : 2, 'cehestarflag' : 0, 'cemergeflag' : 0, 'ecsn' : 2.5, 'ecsn_mlow' : 1.4, 'aic' : 1, 'ussn' : 0, 'sigmadiv' :-20.0, 'qcflag' : 3, 'eddlimflag' : 0, 'fprimc_array' : [2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0], 'rembar_massloss' : 0.5, 'zsun' : 0.02, 'kickflag' : 1, 'grflag' : 1, 'acc_lim' : -1, 'don_lim' : -1} + BSEDict = {'xi': 0.5, 'bhflag': 1, 'neta': 0.5, 'windflag': 3, 'rtmsflag' : 0, 'wdflag': 0, 'alpha1': [1.0, 1.0], 'pts1': 0.05, 'pts3': 0.02, 'pts2': 0.01, 'epsnov': 0.001, 'hewind': 1.0, 'ck': 1000, 'bwind': 0.0, 'lambdaf': 0.5, 'mxns': 3.0, 'beta': -1.0, 'tflag': 1, 'acc2': 1.5, 'nsflag': 4, 'ceflag': 0, 'eddfac': 1.0, 'ifflag': 0, 'bconst': 3000, 'sigma': 265.0, 'gamma': -2.0, 'pisn': 45.0, 'natal_kick_array' :[[-100.0,-100.0,-100.0,-100.0,0.0], [-100.0,-100.0,-100.0,-100.0,0.0]], 'bhsigmafrac' : 1.0, 'polar_kick_angle' : 90, 'qcrit_array' : [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], 'cekickflag' : 2, 'cehestarflag' : 0, 'cemergeflag' : 0, 'ecsn' : 2.5, 'ecsn_mlow' : 1.4, 'aic' : 1, 'ussn' : 0, 'sigmadiv' :-20.0, 'qcflag' : 3, 'eddlimflag' : 0, 'fprimc_array' : [2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0], 'rembar_massloss' : 0.5, 'zsun' : 0.02, 'kickflag' : 1, 'grflag' : 1, 'acc_lim' : [-1, -1], 'don_lim' : -1} filters = {'binary_state': [0], 'timestep_conditions' : 'dtp=None'} convergence = {'convergence_params': ['mass_1', 'mass_2', 'sep', 'ecc'], 'pop_select': 'formation',\ 'match': -5.0, 'convergence_limits' : {"sep" : [0,1000]}, 'match' : -3.0,\ From c349d8dd6263ddf7fb428a3d1b196c3d25fccec5 Mon Sep 17 00:00:00 2001 From: Katie Breivik Date: Thu, 30 Apr 2026 13:29:08 -0400 Subject: [PATCH 15/15] small fix to make the docs look right --- docs/pages/evolve/interface.rst | 2 +- src/cosmic/data/cosmic-settings.json | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/pages/evolve/interface.rst b/docs/pages/evolve/interface.rst index c6ead2f90..404b270a6 100644 --- a/docs/pages/evolve/interface.rst +++ b/docs/pages/evolve/interface.rst @@ -236,7 +236,7 @@ new settings. For example, let's say we want to see how changing the common enve .. ipython:: python ce_alpha_10 = output.rerun_with_settings( - new_settings={'alpha1': 10}, inplace=False + new_settings={'alpha1_0': 10, 'alpha1_1': 10}, inplace=False ) n_merger_original = len(output.final_bpp[output.final_bpp["sep"] == 0.0]) n_merger_ce_alpha_10 = len(ce_alpha_10.final_bpp[ce_alpha_10.final_bpp["sep"] == 0.0]) diff --git a/src/cosmic/data/cosmic-settings.json b/src/cosmic/data/cosmic-settings.json index cd286aeb3..b5e4cf497 100644 --- a/src/cosmic/data/cosmic-settings.json +++ b/src/cosmic/data/cosmic-settings.json @@ -736,7 +736,7 @@ "settings-section-description": "Note: there are cases where a common envelope is forced regardless of the critical mass ratio for unstable mass transfer. In the following cases, a common envelope occurs regardless of the choices below:
  • contact: the stellar radii go into contact (common for similar ZAMS systems)
  • periapse contact: the periapse distance is smaller than either of the stellar radii (common for highly eccentric systems)
  • core Roche overflow: either of the stellar radii overflow their component's Roche radius (in this case, mass transfer from the convective core is always dynamically unstable)
", "name": "alpha1", "description": "Common-envelope efficiency parameter which scales the efficiency of transferring orbital energy to the envelope. See Hurley+2002, Eq. 71.", - "type": "number", + "type": "string", "options-preface": "This value is equivalent to \\(\\alpha\\) in Hurley+2002, Eq. 71.", "options": [ { @@ -1628,7 +1628,7 @@ { "name": "acc_lim", "description": "Limits the amount of mass accreted during Roche-lobe overflow where each binary component must have its own limit specified in a list", - "type": "number", + "type": "string", "options-preface": "", "options": [ {