Add Salim+18 modified Calzetti attenuation law (δ, B) with end-to-end ProSpect integration

R/ProSpect.R

@@ -47,6 +47,10 @@ ProSpectSED = function(SFH = SFHfunc,
                           Eb = 0,
                           L0 = 2175.8,
                           LFWHM = 470,
+                          dust_law = 'CF',
+                          delta_screen = 0,
+                          B_screen = 0,
+                          Rv_screen = 4.05,
                           IGMabsorb = 0,
                           Inoue14_LAFcoef = NULL,
                           Inoue14_DLAcoef = NULL,
@@ -74,6 +78,10 @@ ProSpectSED = function(SFH = SFHfunc,
     alpha_SF_screen = .interval(alpha_SF_screen, 0.0625, 4, reflect = FALSE)
   }
 
+  if(!dust_law %in% c('CF', 'Salim18')){
+    stop("dust_law must be one of 'CF' or 'Salim18'")
+  }
+
   if(!is.null(AGN)){
     tau_AGN = .interval(tau_AGN, 0, 10, reflect = FALSE)
     pow_AGN = .interval(pow_AGN, -2, 0, reflect = FALSE)
@@ -96,6 +104,10 @@ ProSpectSED = function(SFH = SFHfunc,
       Eb = Eb,
       L0 = L0,
       LFWHM = LFWHM,
+      dust_law = dust_law,
+      delta_screen = delta_screen,
+      B_screen = B_screen,
+      Rv_screen = Rv_screen,
       ...
     )
 
@@ -213,7 +225,11 @@ ProSpectSED = function(SFH = SFHfunc,
         waveout = waveout,
         Eb = Eb,
         L0 = L0,
-        LFWHM = LFWHM
+        LFWHM = LFWHM,
+        dust_law = dust_law,
+        delta = delta_screen,
+        B = B_screen,
+        Rv = Rv_screen
       )
       if (!is.null(Dale_M2L_func) & returnall) {
         dustlum_screen = dustlum_screen + AGN$total_atten
@@ -280,7 +296,11 @@ ProSpectSED = function(SFH = SFHfunc,
         waveout = waveout,
         Eb = Eb,
         L0 = L0,
-        LFWHM = LFWHM
+        LFWHM = LFWHM,
+        dust_law = dust_law,
+        delta = delta_screen,
+        B = B_screen,
+        Rv = Rv_screen
       )
 
       if (!is.null(Dale_M2L_func) & returnall) {

R/SFH.R

@@ -33,13 +33,21 @@ SFHfunc = function(massfunc = massfunc_b5,
                    Eb = 0,
                    L0 = 2175.8,
                    LFWHM = 470,
+                   dust_law = 'CF',
+                   delta_screen = 0,
+                   B_screen = 0,
+                   Rv_screen = 4.05,
                    SMstar = FALSE,
                    ...) {
   #Ly_limit should be 911.75 Ang (the actual ionisation limit) or sometimes 1215.67 Ang (Lyman alpha)
 
   dots = list(...)
   massfunc_args = dots[names(dots) %in% names(formals(massfunc))]
 
+  if(!dust_law %in% c('CF', 'Salim18')){
+    stop("dust_law must be one of 'CF' or 'Salim18'")
+  }
+
   if (is.null(speclib)) {
     if (stellpop == 'BC03lr') {
       BC03lr = NULL
@@ -413,13 +421,17 @@ SFHfunc = function(massfunc = massfunc_b5,
   }
 
   if (tau_screen != 0) {
-    lum = lum * CF_screen(
+    lum = lum * screen_atten(
       wave_lum,
       tau = tau_screen,
       pow = pow_screen,
       Eb = Eb,
       L0 = L0,
-      LFWHM = LFWHM
+      LFWHM = LFWHM,
+      dust_law = dust_law,
+      delta = delta_screen,
+      B = B_screen,
+      Rv = Rv_screen
     )
     lumtot_screen = (lumtot_unatten - lumtot_birth) - sum(.qdiff(wave_lum) * lum)
 
@@ -808,9 +820,17 @@ SFHburst = function(burstmass = 1e8,
                     Eb = 0,
                     L0 = 2175.8,
                     LFWHM = 470,
+                    dust_law = 'CF',
+                    delta_screen = 0,
+                    B_screen = 0,
+                    Rv_screen = 4.05,
                     ...) {
   burstmass = .interval(burstmass, 0, Inf, reflect = FALSE)
 
+  if(!dust_law %in% c('CF', 'Salim18')){
+    stop("dust_law must be one of 'CF' or 'Salim18'")
+  }
+
   if (stellpop == 'BC03lr') {
     if (is.null(speclib)) {
       BC03lr = NULL
@@ -1025,13 +1045,17 @@ SFHburst = function(burstmass = 1e8,
   }
 
   if (tau_screen != 0) {
-    lum = lum * CF_screen(
+    lum = lum * screen_atten(
       wave_lum,
       tau = tau_screen,
       pow = pow_screen,
       Eb = Eb,
       L0 = L0,
-      LFWHM = LFWHM
+      LFWHM = LFWHM,
+      dust_law = dust_law,
+      delta = delta_screen,
+      B = B_screen,
+      Rv = Rv_screen
     )
     lumtot_screen = (lumtot_unatten - lumtot_birth) - sum(.qdiff(wave_lum) * lum)
   } else{

R/dust.R

@@ -17,6 +17,49 @@ CF_screen = function(wave, tau=0.3, pow=-0.7, pivot=5500, Eb=0, L0=2175.8, LFWHM
   }
 }
 
+.k_calzetti = function(wave, Rv = 4.05){
+  wave_um = wave / 1e4
+  k = rep(NA_real_, length(wave_um))
+
+  blue_sel = wave_um >= 0.12 & wave_um < 0.63
+  red_sel = wave_um >= 0.63 & wave_um <= 2.2
+  low_sel = wave_um < 0.12
+  high_sel = wave_um > 2.2
+
+  if(any(blue_sel)){
+    w = wave_um[blue_sel]
+    k[blue_sel] = 2.659 * (-2.156 + 1.509 / w - 0.198 / w^2 + 0.011 / w^3) + Rv
+  }
+  if(any(red_sel)){
+    w = wave_um[red_sel]
+    k[red_sel] = 2.659 * (-1.857 + 1.040 / w) + Rv
+  }
+  if(any(low_sel)){
+    k[low_sel] = 2.659 * (-2.156 + 1.509 / 0.12 - 0.198 / 0.12^2 + 0.011 / 0.12^3) + Rv
+  }
+  if(any(high_sel)){
+    k[high_sel] = 2.659 * (-1.857 + 1.040 / 2.2) + Rv
+  }
+  return(k)
+}
+
+Salim18_screen = function(wave, tau = 0.3, delta = 0, B = 0, pivot = 5500, Rv = 4.05, L0 = 2175.8, LFWHM = 350){
+  curve = (.k_calzetti(wave = wave, Rv = Rv) + .drude(wave = wave, Eb = B, L0 = L0, LFWHM = LFWHM)) *
+    (wave / pivot)^delta
+  curve_pivot = (.k_calzetti(wave = pivot, Rv = Rv) + .drude(wave = pivot, Eb = B, L0 = L0, LFWHM = LFWHM))
+  return(exp(-tau * curve / curve_pivot))
+}
+
+screen_atten = function(wave, tau=0.3, pow=-0.7, pivot=5500, Eb=0, L0=2175.8, LFWHM=470, dust_law='CF', delta=0, B=0, Rv=4.05){
+  if(dust_law == 'Salim18'){
+    if(B == 0 && Eb != 0){
+      B = Eb
+    }
+    return(Salim18_screen(wave = wave, tau = tau, delta = delta, B = B, pivot = pivot, Rv = Rv, L0 = L0, LFWHM = LFWHM))
+  }
+  return(CF_screen(wave = wave, tau = tau, pow = pow, pivot = pivot, Eb = Eb, L0 = L0, LFWHM = LFWHM))
+}
+
 CF_birth_atten=function(wave, flux, tau=1.0, pow=-0.7, pivot=5500){
   if(tau == 0){return(list(flux=flux, total_atten=0, attenfrac=0))}
   flux_atten = CF_birth(wave, tau=tau, pow=pow, pivot=pivot)*flux
@@ -26,9 +69,9 @@ CF_birth_atten=function(wave, flux, tau=1.0, pow=-0.7, pivot=5500){
   return(list(flux=flux_atten, total_atten=total_atten, attenfrac=atten/unatten))
 }
 
-CF_screen_atten=function(wave, flux, tau=0.3, pow=-0.7, pivot=5500, Eb=0, L0=2175.8, LFWHM=470){
+CF_screen_atten=function(wave, flux, tau=0.3, pow=-0.7, pivot=5500, Eb=0, L0=2175.8, LFWHM=470, dust_law='CF', delta=0, B=0, Rv=4.05){
   if(tau == 0){return(list(flux=flux, total_atten=0, attenfrac=0))}
-  flux_atten = CF_screen(wave, tau=tau, pow=pow, pivot=pivot, Eb=Eb, L0=L0, LFWHM=LFWHM)*flux
+  flux_atten=screen_atten(wave, tau=tau, pow=pow, pivot=pivot, Eb=Eb, L0=L0, LFWHM=LFWHM, dust_law=dust_law, delta=delta, B=B, Rv=Rv)*flux
   unatten = sum(flux*c(0,diff(wave)))
   atten = sum(flux_atten*c(0,diff(wave)))
   total_atten = max(unatten - atten, 0)
@@ -48,11 +91,11 @@ CF_atten=function(wave, flux, tau=0.3, pow=-0.7, pivot=5500){
   return(wave^(-beta)/(850e4^(-beta)))
 }
 
-atten_emit = function(wave, flux, tau=0.3, pow=-0.7, alpha_SF=1.5, Dale=NULL, Dale_M2L_func=NULL, waveout=NULL, Eb=0, L0=2175.8, LFWHM=470){
-  atten = CF_screen_atten(wave=wave, flux=flux, tau=tau, pow=pow, Eb=Eb, L0=L0, LFWHM=LFWHM)
-  emit = Dale_interp(alpha_SF=alpha_SF, AGNfrac = 0, Dale=Dale)
-  emit$Aspec = emit$Aspec*atten$total_atten
-  final = addspec(wave1=wave, flux1=atten$flux, wave2=emit$Wave, flux2=emit$Aspec, extrap=0, waveout=waveout)
+atten_emit=function(wave, flux, tau=0.3, pow=-0.7, alpha_SF=1.5, Dale=NULL, Dale_M2L_func=NULL, waveout=NULL, Eb=0, L0=2175.8, LFWHM=470, dust_law='CF', delta=0, B=0, Rv=4.05){
+  atten=CF_screen_atten(wave=wave, flux=flux, tau=tau, pow=pow, Eb=Eb, L0=L0, LFWHM=LFWHM, dust_law=dust_law, delta=delta, B=B, Rv=Rv)
+  emit=Dale_interp(alpha_SF=alpha_SF, AGNfrac = 0, Dale=Dale)
+  emit$Aspec=emit$Aspec*atten$total_atten
+  final=addspec(wave1=wave, flux1=atten$flux, wave2=emit$Wave, flux2=emit$Aspec, extrap=0, waveout=waveout)
   if(!is.null(Dale_M2L_func)){
     dustmass=atten$total_atten/Dale_M2L_func(alpha_SF)
   }else{

man/ProSpectSED.Rd

@@ -17,9 +17,10 @@ ProSpectSED(SFH = SFHfunc, z = 0.1, tau_birth = 1, tau_screen = 0.3, tau_AGN = 1
   waveout = seq(2, 9.35, by = 0.01), ref, unimax = 1.38e+10, agemax = NULL,
   LumDist_Mpc = NULL, addradio_SF = FALSE, addradio_AGN = FALSE, Te_SF = 10000,
   ff_frac_SF = 0.1, ff_power_SF = -0.1, sy_power_SF = -0.8, Te_AGN = 10000,
-  ff_frac_AGN = 0.1, ff_power_AGN = -0.1, sy_power_AGN = -0.8, AGNct = 40, AGNrm = 60,
-  AGNan = 30, AGNta = 1, AGNal = 4, AGNbe = -0.5, AGNp = 1, AGNq = 1, Eb = 0, L0 = 2175.8,
-  LFWHM = 470, IGMabsorb = 0, Inoue14_LAFcoef = NULL, Inoue14_DLAcoef = NULL, ...)
+   ff_frac_AGN = 0.1, ff_power_AGN = -0.1, sy_power_AGN = -0.8, AGNct = 40, AGNrm = 60,
+   AGNan = 30, AGNta = 1, AGNal = 4, AGNbe = -0.5, AGNp = 1, AGNq = 1, Eb = 0, L0 = 2175.8,
+   LFWHM = 470, dust_law = "CF", delta_screen = 0, B_screen = 0, Rv_screen = 4.05,
+   IGMabsorb = 0, Inoue14_LAFcoef = NULL, Inoue14_DLAcoef = NULL, ...)
 
 ProSpectSEDlike(parm = c(8,9,10,10,0,-0.5,0.2), Data)
 }
@@ -178,6 +179,18 @@ Numeric scalar; location of the 2175.8 Ang dust bump in Angstroms (probably do n
 }
   \item{LFWHM}{
 Numeric scalar; width of the 2175.8 Ang dust bump in Angstroms (probably do not adjust this).
+}
+  \item{dust_law}{
+Character scalar; attenuation law for screen dust. One of \option{'CF'} (default) or \option{'Salim18'}.
+}
+  \item{delta_screen}{
+Numeric scalar; power-law tilt of the Salim et al. (2018) modified Calzetti screen attenuation law.
+}
+  \item{B_screen}{
+Numeric scalar; UV bump amplitude of the Salim et al. (2018) modified Calzetti screen attenuation law.
+}
+  \item{Rv_screen}{
+Numeric scalar; V-band normalisation term for the Salim et al. (2018) modified Calzetti screen attenuation law.
 }
   \item{IGMabsorb}{
 Absorbing fraction for the intervening IGM. This should be thought of as absorption due to material between the galaxy of interest an the observer. This suppressing fraction is applied to wavelength in the Lyman regime (911.8 - 1215.7 Angstrom). Based on Songaila (2004) a resonable approximation for this is given by \option{IGMabsorb} = pnorm(z, mean=3.8, sd=1.2), meaning it is almost 0 (complete transmission) at low redshift and almost 1 (complete absorption) by redshift 8.

man/SFH.Rd

@@ -21,7 +21,8 @@ SFHfunc(massfunc = massfunc_b5, forcemass = FALSE, agescale = 1,
  Ly_limit = 911.75, LKL10 = NULL, disp_stars = FALSE, LSF = NULL, z = 0.1, H0 = 67.8,
  OmegaM = 0.308, OmegaL = 1 - OmegaM, ref, outtype = "mag", sparse = 5, intSFR = FALSE,
  unimax = 1.38e+10, agemax = NULL, LumDist_Mpc = NULL, Eb = 0, L0 = 2175.8,
- LFWHM = 470, SMstar = FALSE, ...)
+ LFWHM = 470, dust_law = "CF", delta_screen = 0, B_screen = 0, Rv_screen = 4.05,
+ SMstar = FALSE, ...)
 
 SMstarfunc(massfunc = massfunc_b5, forcemass = FALSE, agescale = 1,
  burstage = c(0, 1e+08), youngage = c(1e+08, 1e+09), midage = c(1e+09, 5e+09),
@@ -35,7 +36,8 @@ SFHburst(burstmass = 1e8, burstage = 0, stellpop = 'BC03lr', speclib = NULL,
  emission_scale = "FUV", escape_frac = 1 - emission, Ly_limit = 911.75, LKL10 = NULL,
  disp_stars = FALSE, LSF = NULL, z = 0.1, H0 = 67.8, OmegaM = 0.308, OmegaL = 1 - OmegaM,
  ref, outtype = 'mag', sparse = 5, unimax=  13.8e9, agemax = NULL, LumDist_Mpc = NULL,
- Eb = 0, L0 = 2175.8, LFWHM = 470, ...)
+ Eb = 0, L0 = 2175.8, LFWHM = 470, dust_law = "CF", delta_screen = 0, B_screen = 0,
+ Rv_screen = 4.05, ...)
 }
 %- maybe also 'usage' for other objects documented here.
 \arguments{
@@ -161,6 +163,18 @@ Numeric scalar; location of the 2175.8 Ang dust bump in Angstroms (probably do n
 }
   \item{LFWHM}{
 Numeric scalar; width of the 2175.8 Ang dust bump in Angstroms (probably do not adjust this).
+}
+  \item{dust_law}{
+Character scalar; attenuation law for screen dust. One of \option{'CF'} (default) or \option{'Salim18'}.
+}
+  \item{delta_screen}{
+Numeric scalar; power-law tilt of the Salim et al. (2018) modified Calzetti screen attenuation law.
+}
+  \item{B_screen}{
+Numeric scalar; UV bump amplitude of the Salim et al. (2018) modified Calzetti screen attenuation law.
+}
+  \item{Rv_screen}{
+Numeric scalar; V-band normalisation term for the Salim et al. (2018) modified Calzetti screen attenuation law.
 }
   \item{SMstar}{
 Logical; if you have the \code{ParmOff} package (on Github asgr/ParmOff, required for correct argument matching between the various functions used) setting this to TRUE will also compute the various stellar mass calculations of \code{\link{SMstarfunc}} including stellar mass remaining (TotSMstar) which is the reason this function is usually run. The output appears in the list object 'SMstar'.

man/dust.Rd

@@ -5,6 +5,8 @@
 \alias{CF_birth_atten}
 \alias{CF_screen}
 \alias{CF_screen_atten}
+\alias{Salim18_screen}
+\alias{screen_atten}
 \alias{atten_emit}
 \alias{dust}
 %- Also NEED an '\alias' for EACH other topic documented here.
@@ -21,10 +23,15 @@ CF_birth(wave, tau = 1.0, pow = -0.7, pivot = 5500)
 CF_birth_atten(wave, flux, tau = 1.0, pow = -0.7, pivot = 5500)
 CF_screen(wave, tau = 0.3, pow = -0.7, pivot = 5500, Eb = 0, L0 = 2175.8,
   LFWHM = 470)
-CF_screen_atten(wave, flux, tau = 0.3, pow = -0.7, pivot = 5500, Eb = 0,
-  L0 = 2175.8, LFWHM = 470)
+Salim18_screen(wave, tau = 0.3, delta = 0, B = 0, pivot = 5500, Rv = 4.05,
+  L0 = 2175.8, LFWHM = 350)
+screen_atten(wave, tau = 0.3, pow = -0.7, pivot = 5500, Eb = 0, L0 = 2175.8,
+  LFWHM = 470, dust_law = "CF", delta = 0, B = 0, Rv = 4.05)
+CF_screen_atten(wave, flux, tau = 0.3, pow = -0.7, pivot = 5500, Eb = 0, L0 = 2175.8,
+  LFWHM = 470, dust_law = "CF", delta = 0, B = 0, Rv = 4.05)
 atten_emit(wave, flux, tau = 0.3, pow = -0.7, alpha_SF = 1.5, Dale = NULL,
-Dale_M2L_func = NULL, waveout = NULL, Eb = 0, L0 = 2175.8, LFWHM = 470) 
+Dale_M2L_func = NULL, waveout = NULL, Eb = 0, L0 = 2175.8, LFWHM = 470,
+dust_law = "CF", delta = 0, B = 0, Rv = 4.05) 
 }
 %- maybe also 'usage' for other objects documented here.
 \arguments{
@@ -51,6 +58,18 @@ Numeric scalar; location of the 2175.8 Ang dust bump in Angstroms (probably do n
 }
   \item{LFWHM}{
 Numeric scalar; width of the 2175.8 Ang dust bump in Angstroms (probably do not adjust this). 
+}
+  \item{delta}{
+Numeric scalar; power-law tilt for the modified Calzetti attenuation law (Salim+18).
+}
+  \item{B}{
+Numeric scalar; amplitude of the UV bump for the modified Calzetti attenuation law (Salim+18).
+}
+  \item{Rv}{
+Numeric scalar; V-band normalisation term in the Calzetti law (default 4.05).
+}
+  \item{dust_law}{
+Character scalar; attenuation law for screen dust. One of \option{'CF'} (default) or \option{'Salim18'}.
 }
   \item{alpha_SF}{
 Numeric scalar; desired interpolated alpha slope of the star forming population. Lower values mean hotter dust.  
@@ -71,13 +90,16 @@ We use the simple dust models for birth clouds and dust screens as given by Char
 \deqn{A=\exp(-\tau*(\lambda/\lambda_{piv})^n)}{A=exp(-tau*(wave/pivot)^pow)}
 
 The defaults should be reasonable in the regime they are used over, where birth cloud dust should only be applied to stellar populations younger than 10 Myrs. Note that BC03 has 70 stellar population spectra which are younger than this, whilst EMILES has none. This means that EMILES cannot realistically capture the birth dust attenuation phase, so for highly star bursting populations BC03 should probably be preferred.
+
+For screen dust, the Salim et al. (2018) modified Calzetti law is also available via \code{Salim18_screen} and \code{screen_atten(..., dust_law='Salim18')}, with free slope (\option{delta}) and UV bump strength (\option{B}).
+For this law \option{LFWHM} defaults to 350 Ang in \code{Salim18_screen}, matching common modified Calzetti/Noll-style bump choices, while the existing Charlot \& Fall bump-related defaults remain unchanged.
 }
 \value{
-For \code{CF}, \code{CF_birth} and \code{CF_screen}, a numeric vector; the attenuation curve that fluxes should be multiplied by.
+For \code{CF}, \code{CF_birth}, \code{CF_screen}, \code{Salim18_screen} and \code{screen_atten}, a numeric vector; the attenuation curve that fluxes should be multiplied by.
 
 For \code{CF_atten}, \code{CF_birth_atten} and \code{CF_screen_atten}, a list where the first element is the attenuated flux [flux], the second is the total energy of the flux attenuated (i.e. this is the amount that should be re-radiated elsewhere) [total_atten], and the third is the ratio of attenuated to unattenuated stellar light (i.e. 1 would mean no flux has been attenuated by dust, and 0 would mean 100\% attenuation).
 
-For \code{atten_emit} the input spectrum is attenuated by the specified Charlot and Fall mode, and then re-emitted with the specified Dale dust template. Returns a list containing:
+For \code{atten_emit} the input spectrum is attenuated by the specified screen dust law, and then re-emitted with the specified Dale dust template. Returns a list containing:
 
 \describe{
   \item{final}{Two column data.frame; Wavelength in Ang, attenuated and re-emitted flux in user units.}
@@ -90,7 +112,8 @@ For \code{atten_emit} the input spectrum is attenuated by the specified Charlot
 }
 \references{
 Charlot & Fall, 2000, ApJ, 539, 718 \cr
-da Cunha et al, 2008, MNRAS, 388, 1595
+da Cunha et al, 2008, MNRAS, 388, 1595 \cr
+Salim et al., 2018, ApJ, 859, 11
 }
 \author{
 Aaron Robotham