diff --git a/src/cosmic/_commit_hash.py b/src/cosmic/_commit_hash.py
new file mode 100644
index 000000000..765e68a4e
--- /dev/null
+++ b/src/cosmic/_commit_hash.py
@@ -0,0 +1 @@
+COMMIT_HASH = "c165703dafee2fb2402d6221f4c73c5e493c5ab9"
diff --git a/src/cosmic/data/cosmic-settings.json b/src/cosmic/data/cosmic-settings.json
index a1b9ca26a..743c3c8dd 100644
--- a/src/cosmic/data/cosmic-settings.json
+++ b/src/cosmic/data/cosmic-settings.json
@@ -567,6 +567,18 @@
"name": 3,
"description": "Same as 2, but LBV-like mass loss for giants and non-degenerate stars beyond the Humphreys-Davidson limit",
"default": true
+ },
+ {
+ "name": 5,
+ "description": "Same as 3, but multiplied by 0.33"
+ },
+ {
+ "name": 6,
+ "description": "Bjorklund 2023 https://doi.org/10.1051/0004-6361/202141948"
+ },
+ {
+ "name": 7,
+ "description": "Kritcka, Kubat, Kritckova 2023 10.48550/arXiv.2311.01257"
}
]
},
@@ -1226,6 +1238,10 @@
"name": 6,
"description": "follows the prescription from Maltsev+2025, with additions from Willcox+2025 for the definition of BH masses.",
"version_added": "3.7.0"
+ },
+ {
+ "name": 7,
+ "description": "Follows Fryer+2022 remnant mass prescription"
}
]
},
@@ -1388,6 +1404,10 @@
"description": "Uses the top-down approach from Renzo+2022, with adjustable parameters from Hendriks+2023",
"version_added": "3.7.2"
},
+ {
+ "name": -5,
+ "description": "Uses Belczynski 'weak' PISN pulsations"
+ },
{
"name": "positive values",
"description": "turns on pulsational pair instability and pair instability SNe, and sets the maximum mass of the allowed remnant (i.e., the bottom of the pair instability mass gap). He core masses between pisn and 65 Msun are assumed to go through pulsational pair instability and limit the He core mass to pisn, while He core masses from 65-135 Msun are assumed have a pair instability SN and leave no remnant."
diff --git a/src/cosmic/src/assign_remnant.f b/src/cosmic/src/assign_remnant.f
index 803eee8d5..fedf83e7d 100644
--- a/src/cosmic/src/assign_remnant.f
+++ b/src/cosmic/src/assign_remnant.f
@@ -11,6 +11,7 @@ SUBROUTINE assign_remnant(zpars,mc,mcbagb,mass,kidx,mt,kw,bhspin)
real*8 mc,mcbagb,mass,mt,mc_tot,met
real*8 frac,kappa,sappa,alphap,polyfit
real*8 m_proto,m_FeNi,m_fb,bhspin,mrem,mch,dMppi
+ real*8 fmix, mcritnsbh, mtemp1, mtemp2
integer kw,kidx
* Inputs
@@ -85,6 +86,37 @@ SUBROUTINE assign_remnant(zpars,mc,mcbagb,mass,kidx,mt,kw,bhspin)
*
kw = 15
else
+* Beginning of supernova block
+*
+* Chris Belczynski Evolutionary Roads Weak PPISN
+* This has to happen before the SNa, because it modifies
+* the properties of the star during explosion
+ if(pisn.eq.-5.and.mt.ge.45d0)then
+ if(mcbagb.ge.65d0) then
+ mt = 0.d0
+ kw = 15
+ else
+* PPISN
+ if(mcbagb.ge.60d0) then
+ mtemp1 = 938d0 - (14.3d0*mcbagb)
+ elseif(mcbagb.ge.40d0) then
+ mtemp1 = 55.6d0
+ else
+ mtemp1 = 6.0d0 + (0.83d0*mcbagb)
+ endif
+* Update mass
+ if(mt.gt.mtemp1) then
+ mt = mtemp1
+ endif
+ if(mcbagb.gt.mtemp1) then
+ mcbagb = mtemp1
+ endif
+ if(mc.gt.mtemp1) then
+ mc = mtemp1
+ endif
+ endif
+ endif
+* Carry on with the Supernovae
*
* Use remnant mass given by Hurley+2000
if(remnantflag.eq.0)then
@@ -217,6 +249,48 @@ SUBROUTINE assign_remnant(zpars,mc,mcbagb,mass,kidx,mt,kw,bhspin)
elseif(remnantflag.eq.6)then
met = 10**(LOG10(zpars(14))/0.4)
call assign_remnant_maltsev(mc,mc_tot,met,kidx,kw,mt)
+ elseif(remnantflag.eq.7)then
+*
+* Use the Fryer et al. 2022 SN Prescription
+*
+* For this, we just set the proto-core mass to one
+ if(mc.le.3.5d0)then
+ m_proto = 1.2d0
+ elseif(mc.le.6.d0)then
+ m_proto = 1.3d0
+ elseif(mc.le.11.d0)then
+ m_proto = 1.4d0
+ elseif(mc.gt.11.d0)then
+ m_proto = 1.6d0
+ endif
+
+ if(ecsn.gt.0.d0.and.mcbagb.le.ecsn.and.
+ & mcbagb.ge.ecsn_mlow)then
+ mt = 1.38d0 ! ECSN fixed mass, no fallback
+ else
+* Parameters of Fryer2022 model
+ fmix=1.0
+ mcritnsbh=5.75
+* We need mt in multiple places, so temp1 will be the working mt
+ mtemp1=mt
+* mtemp2 is the calculated value of the remnant mass
+ mtemp2=1.2 + (0.05*fmix) +
+ & (0.01*((mc/fmix)**2)) +
+ & EXP(fmix*(mc-mcritnsbh))
+* We don't care about mtemp2 if it's less than zero
+ if(mtemp2.lt.0.)then
+ mtemp1 = 0.
+ kw=15
+* We only care about mtemp2 if it is less than the total
+* mass of the star
+ elseif(mtemp2.lt.mt)then
+ mtemp1 = mtemp2
+* If mtemp2 is less, we also want to estimate the fallback fraction
+ fallback=(mtemp1-m_proto) /(mt-m_proto)
+ mt = m_proto + fallback*(mtemp1 - m_proto)
+ endif
+ endif
+ mc = mt
endif
* Assign the BH spin based on the chosen prescription
@@ -747,4 +821,4 @@ integer function first_mt_type_as_donor(star)
endif
return
- end
\ No newline at end of file
+ end
diff --git a/src/cosmic/src/comenv.f b/src/cosmic/src/comenv.f
index 83e7ee744..1acfe3e43 100644
--- a/src/cosmic/src/comenv.f
+++ b/src/cosmic/src/comenv.f
@@ -56,8 +56,6 @@ 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
@@ -323,9 +321,10 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1,
& aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp,
& mc_he(1),mc_he(2),mc_co(1),mc_co(2),
& rad1_bpp,rad2_bpp,
- & M02,mass_preSN,lumin(1),lumin(2),
+ & M02,M01,lumin(1),lumin(2),
& teff1,teff2,
- & RC2,RC1,menv_bpp(1),mHe_preSN,renv_bpp(1),
+ & 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,
@@ -343,9 +342,10 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1,
& aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp,
& mc_he(1),mc_he(2),mc_co(1),mc_co(2),
& rad1_bpp,rad2_bpp,
- & mass_preSN,M02,lumin(1),lumin(2),
+ & M01,M02,lumin(1),lumin(2),
& teff1,teff2,
- & RC1,RC2,mHe_preSN,menv_bpp(2),renv_bpp(1),
+ & RC1,RC2,menv_bpp(1),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,
@@ -632,9 +632,10 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1,
& aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp,
& mc_he(1),mc_he(2),mc_co(1),mc_co(2),
& rad1_bpp,rad2_bpp,
- & M02,mass_preSN,lumin(1),lumin(2),
+ & M02,M01,lumin(1),lumin(2),
& teff1,teff2,
- & RC2,RC1,menv_bpp(1),mHe_preSN,renv_bpp(1),
+ & 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,
@@ -652,9 +653,10 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1,
& aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp,
& mc_he(1),mc_he(2),mc_co(1),mc_co(2),
& rad1_bpp,rad2_bpp,
- & mass_preSN,M02,lumin(1),lumin(2),
+ & M01,M02,lumin(1),lumin(2),
& teff1,teff2,
- & RC1,RC2,mHe_preSN,menv_bpp(2),renv_bpp(1),
+ & RC1,RC2,menv_bpp(1),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,
@@ -805,9 +807,10 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1,
& aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp,
& mc_he(1),mc_he(2),mc_co(1),mc_co(2),
& rad1_bpp,rad2_bpp,
- & mass_preSN,M01,lumin(1),lumin(2),
+ & M02,M01,lumin(1),lumin(2),
& teff1,teff2,
- & RC2,RC1,mHe_preSN,menv_bpp(2),renv_bpp(1),
+ & 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,
@@ -825,9 +828,10 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1,
& aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp,
& mc_he(1),mc_he(2),mc_co(1),mc_co(2),
& rad1_bpp,rad2_bpp,
- & M01,mass_preSN,lumin(1),lumin(2),
+ & M01,M02,lumin(1),lumin(2),
& teff1,teff2,
- & RC1,RC2,menv_bpp(1),mHe_preSN,renv_bpp(1),
+ & RC1,RC2,menv_bpp(1),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,
@@ -1041,9 +1045,10 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1,
& aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp,
& mc_he(1),mc_he(2),mc_co(1),mc_co(2),
& rad1_bpp,rad2_bpp,
- & M02,mass_preSN,lumin(1),lumin(2),
+ & M02,M01,lumin(1),lumin(2),
& teff1,teff2,
- & RC2,RC1,menv_bpp(1),mHe_preSN,renv_bpp(1),
+ & 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,
@@ -1061,9 +1066,10 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1,
& aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp,
& mc_he(1),mc_he(2),mc_co(1),mc_co(2),
& rad1_bpp,rad2_bpp,
- & mass_preSN,M02,lumin(1),lumin(2),
+ & M01,M02,lumin(1),lumin(2),
& teff1,teff2,
- & RC1,RC2,mHe_preSN,menv_bpp(2),renv_bpp(1),
+ & RC1,RC2,menv_bpp(1),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,
@@ -1115,8 +1121,9 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1,
& jp,tphys,11.d0,M1,M2,KW1,KW2,-1.d0,-1.d0,-1.d0,0.d0,
& 0.d0,aj1_bpp,aj2_bpp,tms1_bpp,tms2_bpp,mc_he(1),
& mc_he(2),mc_co(1),mc_co(2),rad(1),rad(2),M01,M02,lumin(1),
- & lumin(2),teff1,teff2,RC1,RC2,MENV,mHe_preSN,renv_bpp(1),
- & renv_bpp(2),OSPIN1,OSPIN2,B_0(1),B_0(2),bacc(1),bacc(2),
+ & lumin(2),teff1,teff2,RC1,RC2,menv_bpp(1),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,2,-1,
& zpars(14)**2.d5,'bpp')
diff --git a/src/cosmic/src/evolv2.f b/src/cosmic/src/evolv2.f
index 9084d5265..04a0b535d 100644
--- a/src/cosmic/src/evolv2.f
+++ b/src/cosmic/src/evolv2.f
@@ -197,8 +197,6 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf,
COMMON /fall/fallback
REAL ran3
EXTERNAL ran3
-*
-
*
REAL*8 z,tm,tn,m0,mt,rm,lum,mc,rc,me,re,k2,age,dtm,dtr
REAL*8 tscls(20),lums(10),GB(10),zpars(20)
diff --git a/src/cosmic/src/mlwind.f b/src/cosmic/src/mlwind.f
index df1f16449..9d9334fe0 100644
--- a/src/cosmic/src/mlwind.f
+++ b/src/cosmic/src/mlwind.f
@@ -116,7 +116,8 @@ real*8 FUNCTION mlwind(kw,lum,r,mt,mc,rl,z)
endif
*
mlwind = dms
- elseif(windflag.eq.2.or.windflag.eq.3)then
+ elseif(windflag.eq.2.or.windflag.eq.3.or.windflag.eq.5
+ & .or.windflag.eq.6.or.windflag.eq.7)then
* Vink winds etc according to as implemented following
* Belczynski, Bulik, Fryer, Ruiter, Valsecchi, Vink & Hurley 2010.
*
@@ -159,22 +160,47 @@ real*8 FUNCTION mlwind(kw,lum,r,mt,mc,rl,z)
endif
* Apply Vink, de Koter & Lamers (2001) OB star winds.
* Next check if hot massive H-rich O/B star in appropriate temperature ranges.
- if(teff.ge.12500.and.teff.le.25000)then
- if(eddlimflag.eq.0) alpha = 0.85d0
- if(eddlimflag.eq.1) alpha = MLalpha(mt,lum,kw)
- dms = -6.688d0 + 2.210d0*LOG10(lum/1.0d+05) -
- & 1.339d0*LOG10(mt/30.d0) - 1.601d0*LOG10(1.3d0/2.d0) +
- & alpha*LOG10(z/zsun) + 1.07d0*LOG10(teff/2.0d+04)
- dms = 10.d0**dms
- elseif(teff.gt.25000.)then
-* Although Vink et al. formulae are only defined until Teff=50000K,
-* we follow the Dutch prescription of MESA, and extend to higher Teff
- dms = -6.697d0 + 2.194d0*LOG10(lum/1.0d+05) -
- & 1.313d0*LOG10(mt/30.d0) - 1.226d0*LOG10(2.6d0/2.d0) +
- & alpha*LOG10(z/zsun) +0.933d0*LOG10(teff/4.0d+04) -
- & 10.92d0*(LOG10(teff/4.0d+04)**2)
- dms = 10.d0**dms
+ if (windflag.eq.2.or.windflag.eq.3.or.windflag.eq.5) then
+ if(teff.ge.12500.and.teff.le.25000)then
+ if(eddlimflag.eq.0) alpha = 0.85d0
+ if(eddlimflag.eq.1) alpha = MLalpha(mt,lum,kw)
+ dms = -6.688d0 + 2.210d0*LOG10(lum/1.0d+05) -
+ & 1.339d0*LOG10(mt/30.d0) -
+ & 1.601d0*LOG10(1.3d0/2.d0) +
+ & alpha*LOG10(z/zsun) + 1.07d0*LOG10(teff/2.0d+04)
+ dms = 10.d0**dms
+ elseif(teff.gt.25000.)then
+* Although Vink et al. formulae are only defined until Teff=50000K,
+* we follow the Dutch prescription of MESA, and extend to higher Teff
+ dms = -6.697d0 + 2.194d0*LOG10(lum/1.0d+05) -
+ & 1.313d0*LOG10(mt/30.d0) -
+ & 1.226d0*LOG10(2.6d0/2.d0) +
+ & alpha*LOG10(z/zsun) +0.933d0*LOG10(teff/4.0d+04) -
+ & 10.92d0*(LOG10(teff/4.0d+04)**2)
+ dms = 10.d0**dms
+ endif
+* Apply fwind cut
+ if (windflag.eq.5) then
+ dms = dms*0.33
+ endif
+ endif
+* Bjorklund 2023 https://doi.org/10.1051/0004-6361/202141948
+ if (windflag.eq.6) then
+ dms = -5.52d0 + 2.39d0*LOG10(lum/1.0d+06) +
+ & 1.48*LOG10(mt/45.0d0) + 2.12d0 *LOG10(teff/4.5d+04) +
+ & (0.75d+0 - 1.87d0*LOG10(teff/4.5d+04))*LOG10(z/0.014d0)
+ dms = 10.0d0**dms
+ endif
+* Kritcka, Kubat, Kritckova 2023 10.48550/arXiv.2311.01257
+ if (windflag.eq.7) then
+ dms = -13.82d0 + 0.358d0*LOG10(z/0.0134d0) +
+ & (1.52d0 - 0.11d0*LOG10(z/0.0134d0))*LOG10(lum/1.0d+06) +
+ & 13.82d0*LOG10((1.0d0 + 0.73*LOG10(z/0.0134d0))*EXP(
+ & -1.0d0 * (teff - 1.416d+4)**2 / (3.58d+3)**2) +
+ & 3.84d0*EXP(-1.0d0*(teff-3.790d+04)**2 / (5.65d+04)**2))
+ dms = 10.0d0**dms
endif
+*
if((windflag.eq.3.or.kw.ge.2).and.kw.le.6)then
* LBV-like mass loss beyond the Humphreys-Davidson limit.