diff --git a/R/AGN.R b/R/AGN.R
index acc73f1..a0ea34f 100644
--- a/R/AGN.R
+++ b/R/AGN.R
@@ -66,3 +66,79 @@ SKIRTOR_interp = function(lum = 1e+44, ct = 40, rm = 60, an = 30, ta = 1, p = 1,
   return(data.frame(wave=SKIRTOR$Wave, lum = agn_spectrum))
 
 }
+
+Fritz_anav_interp = function(lum = 1e+44, ct = 40, rm = 60, ta = 1, al = 4, be = -0.5,
+                             Fritz_anav = NULL){
+  if(is.null(Fritz_anav)){
+    Fritz_anav = NULL
+    data('Fritz_anav', envir = environment())
+    Fritz_anav = get('Fritz_anav', envir = environment())
+  }
+  if(!inherits(Fritz_anav, 'Fritz_anav')){
+    stop('Fritz_anav must be a Fritz_anav object.', call. = FALSE)
+  }
+
+  ctmix = interp_quick(ct, Fritz_anav$ct)
+  almix = interp_quick(al, Fritz_anav$al)
+  bemix = interp_quick(be, Fritz_anav$be)
+  tamix = interp_quick(ta, Fritz_anav$ta)
+  rmmix = interp_quick(rm, Fritz_anav$rm)
+
+  slice = Fritz_anav$Aspec[
+    c(ctmix[1:2]),
+    c(almix[1:2]),
+    c(bemix[1:2]),
+    c(tamix[1:2]),
+    c(rmmix[1:2]),
+  ]
+
+  weights = rep(1,32)
+  weights = weights * rep(ctmix[3:4], each=1, times=16)
+  weights = weights * rep(almix[3:4], each=2, times=8)
+  weights = weights * rep(bemix[3:4], each=4, times=4)
+  weights = weights * rep(tamix[3:4], each=8, times=2)
+  weights = weights * rep(rmmix[3:4], each=16, times=1)
+  tempmat = matrix(as.numeric(slice), 32, length(Fritz_anav$Wave))
+  return(data.frame(wave = Fritz_anav$Wave, lum = lum * colSums(tempmat * weights)))
+}
+
+
+SKIRTOR_anav_interp = function(lum = 1e+44, ct = 40, rm = 60, ta = 1, p = 1, q = 1,
+                               SKIRTOR_anav = NULL){
+
+  if(is.null(SKIRTOR_anav)){
+    SKIRTOR_anav = NULL
+    data('SKIRTOR_anav', envir = environment())
+    SKIRTOR_anav = get('SKIRTOR_anav', envir = environment())
+  }
+  if(!inherits(SKIRTOR_anav, 'SKIRTOR_anav')){
+    stop('SKIRTOR_anav must be a SKIRTOR_anav object.', call. = FALSE)
+  }
+
+  tmix = interp_quick(ta, SKIRTOR_anav$ta)
+  pmix = interp_quick(p, SKIRTOR_anav$p)
+  qmix = interp_quick(q, SKIRTOR_anav$q)
+  oamix = interp_quick(ct, SKIRTOR_anav$ct)
+  rmmix = interp_quick(rm, SKIRTOR_anav$rm)
+
+  slice = SKIRTOR_anav$Aspec[
+    c(tmix[1:2]),
+    c(pmix[1:2]),
+    c(qmix[1:2]),
+    c(oamix[1:2]),
+    c(rmmix[1:2]),
+    1,
+  ]
+
+  weights = rep(1,32)
+  weights = weights * rep(tmix[3:4], each=1, times=16)
+  weights = weights * rep(pmix[3:4], each=2, times=8)
+  weights = weights * rep(qmix[3:4], each=4, times=4)
+  weights = weights * rep(oamix[3:4], each=8, times=2)
+  weights = weights * rep(rmmix[3:4], each=16, times=1)
+  tempmat = matrix(as.numeric(slice), 32, length(SKIRTOR_anav$Wave))
+  agn_spectrum = colSums(tempmat * weights)/SKIRTOR_anav$Wave
+  agn_spectrum = (agn_spectrum/(3.839e33 * 1E11))*lum
+  return(data.frame(wave = SKIRTOR_anav$Wave, lum = agn_spectrum))
+
+}
diff --git a/R/ProSpect.R b/R/ProSpect.R
index f55c19d..2169e2c 100644
--- a/R/ProSpect.R
+++ b/R/ProSpect.R
@@ -50,9 +50,17 @@ ProSpectSED = function(SFH = SFHfunc,
                           IGMabsorb = 0,
                           Inoue14_LAFcoef = NULL,
                           Inoue14_DLAcoef = NULL,
+                          AGN_screen_model = 'viewing_angle',
+                          AGN_anav = NULL,
                           ...) {
   #call = match.call()
 
+  if (!is.character(AGN_screen_model) || length(AGN_screen_model) != 1L ||
+      is.na(AGN_screen_model) ||
+      !AGN_screen_model %in% c('viewing_angle', 'angle_averaged', 'None')) {
+    stop("AGN_screen_model must be 'viewing_angle', 'angle_averaged', or 'None'.", call. = FALSE)
+  }
+
   if(!is.null(waveout)){
     waveout_max = max(waveout)
   }else{
@@ -173,10 +181,21 @@ ProSpectSED = function(SFH = SFHfunc,
       stop('Dale cannot be FALSE when using an AGN model!')
     }
 
+    if (AGN_screen_model == 'angle_averaged' &&
+        !inherits(AGN, 'Fritz') && !inherits(AGN, 'SKIRTOR')) {
+      stop('Angle-averaged AGN emission requires a Fritz or SKIRTOR model')
+    }
+
     if(inherits(AGN, 'Fritz') |
        inherits(AGN, 'SKIRTOR')) {
       #Use new model
       if (inherits(AGN, 'Fritz')) {
+        if (AGN_screen_model == 'angle_averaged') {
+          AGN_mean = Fritz_anav_interp(
+            lum = AGNlum, ct = AGNct, rm = AGNrm, ta = AGNta,
+            al = AGNal, be = AGNbe, Fritz_anav = AGN_anav
+          )
+        }
         AGN = Fritz_interp(
           lum = AGNlum,
           ct = AGNct,
@@ -188,6 +207,12 @@ ProSpectSED = function(SFH = SFHfunc,
           Fritz = AGN
         )
       } else if (inherits(AGN, 'SKIRTOR')) {
+        if (AGN_screen_model == 'angle_averaged') {
+          AGN_mean = SKIRTOR_anav_interp(
+            lum = AGNlum, ta = AGNta, p = AGNp, q = AGNq,
+            ct = AGNct, rm = AGNrm, SKIRTOR_anav = AGN_anav
+          )
+        }
         AGN = SKIRTOR_interp(
           lum = AGNlum,
           ta = AGNta,
@@ -203,10 +228,18 @@ ProSpectSED = function(SFH = SFHfunc,
       dustlum_AGN = NA
       dustmass_AGN = NA
       #Second we re-attenuate the above by the screen (since it still has to pass out of the galaxy)
+      if (AGN_screen_model == 'angle_averaged') {
+        direct = CF_screen_atten(
+          wave = AGN$wave, flux = AGN$lum * .erg_to_lsol,
+          tau = tau_screen, pow = pow_screen, Eb = Eb, L0 = L0, LFWHM = LFWHM
+        )
+        selected = AGN
+        AGN = AGN_mean
+      }
       AGN = atten_emit(
         wave = AGN$wave,
         flux = AGN$lum * .erg_to_lsol,
-        tau = tau_screen,
+        tau = if (AGN_screen_model == 'None') 0 else tau_screen,
         pow = pow_screen,
         alpha_SF = alpha_SF_screen,
         Dale = Dale,
@@ -216,6 +249,14 @@ ProSpectSED = function(SFH = SFHfunc,
         L0 = L0,
         LFWHM = LFWHM
       )
+      if (AGN_screen_model == 'angle_averaged') {
+        #The sightline sets direct light; the angular mean sets dust heating.
+        AGN$final = addspec(
+          wave1 = selected$wave, flux1 = direct$flux,
+          wave2 = AGN$emit$wave, flux2 = AGN$emit$flux,
+          extrap = 0, waveout = waveout
+        )
+      }
       if (!is.null(Dale_M2L_func) & returnall) {
         dustlum_screen = dustlum_screen + AGN$total_atten
         dustmass_screen = dustmass_screen + AGN$dustmass
@@ -273,7 +314,7 @@ ProSpectSED = function(SFH = SFHfunc,
       AGN = atten_emit(
         wave = AGN$final$wave,
         flux = AGN$final$flux,
-        tau = tau_screen,
+        tau = if (AGN_screen_model == 'None') 0 else tau_screen,
         pow = pow_screen,
         alpha_SF = alpha_SF_screen,
         Dale = Dale,
@@ -500,6 +541,20 @@ ProSpectSED = function(SFH = SFHfunc,
 }
 
 ProSpectSEDlike = function(parm = c(8, 9, 10, 10, 0, -0.5, 0.2), Data) {
+  if (!is.null(Data$AGN_screen_model)) {
+    Data$arglist$AGN_screen_model = Data$AGN_screen_model
+  }
+  if (!is.null(Data$AGN_anav)) {
+    Data$arglist$AGN_anav = Data$AGN_anav
+  }
+  if (!is.null(Data$arglist$AGN_screen_model) &&
+      (!is.character(Data$arglist$AGN_screen_model) ||
+       length(Data$arglist$AGN_screen_model) != 1L ||
+       is.na(Data$arglist$AGN_screen_model) ||
+       !Data$arglist$AGN_screen_model %in% c('viewing_angle', 'angle_averaged', 'None'))) {
+    stop("Data$AGN_screen_model must be 'viewing_angle', 'angle_averaged', or 'None'.", call. = FALSE)
+  }
+
   if (is.null(Data$fit)) {
     Data$fit = 'optim'
   }
diff --git a/R/plots.R b/R/plots.R
index 11b51ce..2477b7a 100644
--- a/R/plots.R
+++ b/R/plots.R
@@ -12,7 +12,7 @@ plot.ProSpectSED = function(x,
     if (ylim[1] == 'auto') {
       if(!is.null(x$Stars)){
         ylim = c(quantile(x$FinalLum[, 2], 0.45),
-                 max(x$StarsUnAtten, na.rm = TRUE))
+                 max(x$StarsUnAtten$lum, na.rm = TRUE))
       }else{
         ylim = c(quantile(x$FinalLum[, 2], 0.45),
                  max(x$FinalLum[, 2], na.rm = TRUE))
diff --git a/man/AGNinterp.Rd b/man/AGNinterp.Rd
index 54191fc..8a3a2de 100644
--- a/man/AGNinterp.Rd
+++ b/man/AGNinterp.Rd
@@ -2,12 +2,18 @@
 \alias{AGNinterp}
 \alias{Fritz_interp}
 \alias{SKIRTOR_interp}
+\alias{Fritz_anav_interp}
+\alias{SKIRTOR_anav_interp}
 %- Also NEED an '\alias' for EACH other topic documented here.
 \title{
 Interpolate the Fritz (2006) and SKIRTOR (2016) AGN Templates
 }
 \description{
 Interpolates between non-computed values of the Fritz 2006 AGN Templates using a mixture of schemes. This is similar to the approach used for \code{\link{Dale_interp}}. The model is much too complex to summarise here, so interested readers should refer to the source papers if using this function for science.
+
+\code{Fritz_anav_interp} and \code{SKIRTOR_anav_interp} interpolate the
+angle-averaged template libraries. They have no viewing-angle argument and
+retain the luminosity normalisation of the corresponding viewing-angle models.
 }
 \usage{
 Fritz_interp(lum = 1e+44, ct = 40, rm = 60, an = 30, ta = 1, al = 4, be = -0.5,
@@ -15,6 +21,12 @@ Fritz_interp(lum = 1e+44, ct = 40, rm = 60, an = 30, ta = 1, al = 4, be = -0.5,
 
 SKIRTOR_interp(lum = 1e+44, ct = 40, rm = 60, an = 30, ta = 1, p = 1, q = 1,
   SKIRTOR = NULL)
+
+Fritz_anav_interp(lum = 1e+44, ct = 40, rm = 60, ta = 1, al = 4, be = -0.5,
+  Fritz_anav = NULL)
+
+SKIRTOR_anav_interp(lum = 1e+44, ct = 40, rm = 60, ta = 1, p = 1, q = 1,
+  SKIRTOR_anav = NULL)
 }
 \arguments{
   \item{lum}{
@@ -49,6 +61,14 @@ The \option{Fritz} data. Providing this directly reduces disk IO since it will u
 }
   \item{SKIRTOR}{
 The \option{SKIRTOR} data. Providing this directly reduces disk IO since it will use the copy already loaded into RAM. If not provided it will load the data from disk.
+}
+  \item{Fritz_anav}{
+The angle-averaged \code{Fritz_anav} template library. If NULL, it is loaded
+with \code{data("Fritz_anav")}. Pass a preloaded library to avoid repeated disk IO.
+}
+  \item{SKIRTOR_anav}{
+The angle-averaged \code{SKIRTOR_anav} template library. If NULL, it is loaded
+with \code{data("SKIRTOR_anav")}. Pass a preloaded library to avoid repeated disk IO.
 }
 }
 \details{
@@ -70,7 +90,9 @@ Previously, we have found that an advantageous property of such two-phase medium
 If \option{an} this is larger then \option{ct}/2 then by definition a significant component of the AGN accretion disk light is being attenuated by the torus (AKA type II AGN). If \option{an} is less then \option{ct}/2 then we have a direct line-of-sight to the AGN accretion disk and will observe much more UV flux in particular (AKA type I AGN, including blazars and quasars).
 }
 \value{
-Numeric Matrix; first column is observed wavelength (Angstroms) and second column is observed luminosity (erg/s/Ang).
+Data frame with columns \code{wave} (rest-frame wavelength in Angstroms) and
+\code{lum} (luminosity density in erg/s/Ang). The angle-averaged helpers return
+angle-averaged luminosity density rather than a particular viewing direction.
 }
 \references{
 Feltre et al, 2012, MNRAS, 426, 120 \cr
diff --git a/man/ProSpectSED.Rd b/man/ProSpectSED.Rd
index b25e4e9..0547fd7 100644
--- a/man/ProSpectSED.Rd
+++ b/man/ProSpectSED.Rd
@@ -19,7 +19,8 @@ ProSpectSED(SFH = SFHfunc, z = 0.1, tau_birth = 1, tau_screen = 0.3, tau_AGN = 1
   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, ...)
+  LFWHM = 470, IGMabsorb = 0, Inoue14_LAFcoef = NULL, Inoue14_DLAcoef = NULL,
+  AGN_screen_model = "viewing_angle", AGN_anav = NULL, ...)
 
 ProSpectSEDlike(parm = c(8,9,10,10,0,-0.5,0.2), Data)
 }
@@ -72,6 +73,12 @@ Object or scalar logical; optional. Pass in the Dale dust library \code{\link[Pr
 }
   \item{AGN}{
 Object; optional. Pass in the AGN template to use. Must be one of \code{\link{AGN_UnOb}}, \code{\link{AGN_UnOb_Sparse}}, \code{\link[ProSpectData:Fritz]{Fritz}}, or \code{\link[ProSpectData:SKIRTOR]{SKIRTOR}}. These objects have classes attached, when ensures the correct AGN code is executed (so users should not simply try to load their own specific AGN model).
+}
+  \item{AGN_screen_model}{
+Character scalar; \code{"viewing_angle"} (default), \code{"angle_averaged"}, or \code{"None"}. Choice of AGN SED to calculate FIR emission due to UV attenuation by the ISM. The angle-averaged mode internally uses the matching \code{Fritz_anav} or \code{SKIRTOR_anav} template. \code{None} is equivalent to zero screen optical depth for the AGN only.
+  }
+  \item{AGN_anav}{
+Object; optional. Pass in the angle-averaged AGN template to use when \option{AGN_screen_model} = \code{"angle_averaged"}. Must be one of \code{Fritz_anav} or \code{SKIRTOR_anav}. Providing this avoids reloading the averaged template from disk on each call.
 }
 
   \item{filtout}{
@@ -222,7 +229,18 @@ List, optional but recommended; pass in the spectral library directly. Must be o
 List, optional but recommended; pass in the Dale dust library directly (should be one of \option{Dale_NormTot}, \option{Dale_NormSFR} otherwise will run, but not work correctly). If set to FALSE then there will be no FIR re-emission.
 }
   \item{Data$AGN}{
-Data.frame, optional; pass in the AGN template to use. Must be one of \code{\link{AGN_UnOb}} or \code{\link{AGN_UnOb_Sparse}}.
+Object, optional; pass in the viewing-angle AGN template to use. Accepts the
+same template classes as \option{AGN}, including Fritz and SKIRTOR.
+}
+  \item{Data$AGN_screen_model}{
+Character scalar, optional; AGN emission mode passed to \code{ProSpectSED}.
+It can also be supplied as \code{Data$arglist$AGN_screen_model}; a non-NULL
+\code{Data$AGN_screen_model} takes precedence. Defaults to \code{"viewing_angle"}.
+}
+  \item{Data$AGN_anav}{
+Object, optional; matching angle-averaged AGN library passed to
+\code{ProSpectSED}. It can also be supplied as \code{Data$arglist$AGN_anav};
+a non-NULL \code{Data$AGN_anav} takes precedence. See \option{AGN_anav}.
 }
   \item{Data$filtout}{
 List; required. Pass in the photometric filters directly. If \code{filters} is a vector of filter names, this can be achieved efficiently with a command like: filtout=foreach(i = filters)\%do\%{getfilt(i)} or filtout=foreach(i = filters)\%do\%{approxfun(getfilt(i))}, where the latter will be faster (see \code{\link{getfilt}}).