Fixing comenv SN parameters, adding new wind/remnant/ppisn options

src/cosmic/_commit_hash.py

@@ -0,0 +1 @@
+COMMIT_HASH = "c165703dafee2fb2402d6221f4c73c5e493c5ab9"

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 <a class='reference external' href='https://scixplorer.org/abs/2025A%26A...700A..20M/abstract'>Maltsev+2025</a>, with additions from <a class='reference external' href='https://scixplorer.org/abs/2025arXiv251007573W/abstract'>Willcox+2025</a> 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 <a class='reference external' href='https://ui.adsabs.harvard.edu/abs/2022RNAAS...6...25R/abstract'>Renzo+2022</a>, with adjustable parameters from <a class='reference external' href='https://scixplorer.org/abs/2023MNRAS.526.4130H/abstract'>Hendriks+2023</a>",
                         "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."

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
+      end

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')

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)

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.