Learn R Programming

cosmoFns (version 1.1-1)

sedFitThin: Optically-thin SED fit

Description

Function takes Herschel-SPIRE photometry and fits optically-thin greybody function for a single-component temperature and galaxy luminosity. Function generates nsamp realizations of observed flux densities with standard deviations for error analysis.

Usage

sedFitThin(s, e = s*0.2, z = 2.5, nsamp = 100, alpha = 2, beta = 1.5,
wl= c(250, 350, 500), sc.start = 1.e-6, T.start = 50)

Arguments

s

Vector of observed-frame flux densities [Jy]

e

Vector of standard deviation of observed-frame flux density [Jy]

z

Galaxy redshift

nsamp

Number of realizations for Monte-Carlo calculation

alpha

Index of power-law for short-wavelength extension

beta

Dust emissivity power law

wl

Vector of observed-frame wavelengths corresponding to s and e [microns]

sc.start

Initial guess for fit luminosity density scaling factor

T.start

Initial guess for dust temperature [K]

Value

List of class sedfit with elements:

td

Mean of dust temperature distribution

e.td

Standard deviation of dust temperature distribution

lum.gb

Mean of greybody luminosity distribution

e.lum.gb

Standard deviation of greybody luminosity distribution

lum.gbpl

Mean of greybody-power law luminosity distribution

e.lum.gbpl

Standard deviation of greybody-power law luminosity distribution

scaleFactor

Conversion between observed frame flux density and rest frame luminosity density

success

Fraction of fit attempts that converged

results

Matrix with nsamp rows and 5 columns: dust temperature in K, greybody luminosity, luminosity for greybody with smoothly-joined power law to short wavelengths, luminosity density scaling, and transition frequency in GHz for power law extension. The first row contains results for the center-of-error input flux densities s.

Details

Conversion from observed to rest frame is from equation (24) in Hogg 2000. Dust temperature and 8-1000 micron luminosity derivation is described in Blain, Barnard & Chapman 2003. Galaxy SEDs typically fall off more slowly than greybody on the Wien side; see plot generated by examples below to visualize power-law extension suggested by Blain et al. 2003.

References

Hogg 2000, astro-ph 9905116v4; Blain, Barnard & Chapman 2003, MNRAS 338, 733.

Examples

Run this code
# NOT RUN {
s <- c(0.242, 0.293, 0.231)
e <- c(0.037, 0.045, 0.036)
z <- 2.41
beta <- 1.5
alpha <- 2
X <- sedFitThin(s=s, e=e, z=z, alpha=alpha, beta=beta, nsamp=100)
str(X)

## Make a plot
# greybody in blue, power-law extension in red dashed line
# functions
# optically thin greybody
otGreybody <- function(nu, T, beta, sc=1) {
                       # nu in GHz, T in K, beta and sc unitless
                       sc*nu^(3+beta)/(exp(0.04801449*nu/T) - 1)
              }
# high frequency tail
hfTail <- function(nu, alpha) nu^-alpha
#
# setups for 8-1000 microns:
nu.low <- 3e5/1000
nu.high <- 3e5/8
l.nue <- s*X$scaleFactor
#
# greybody
nue.sweep <- seq(nu.low, nu.high, len=350)
pred <- otGreybody(nue.sweep, X$results[1,1], beta=beta,
                   X$results[1,4])
ylim <- range(pred, l.nue)
par(fig=c(0,1,0.2,1), mgp=c(1.8, 0.6, 0))
plot(3e5/nue.sweep, pred, t='l', ylim=ylim, log='xy', col=4,
     xlab='Rest frame wavelength [microns]',
     ylab=expression(paste('Luminosity density [ ', L[sun],
                           ' ', Hz^-1, ']')))
# power law
nue.sweep <- seq(X$results[1,5], nu.high, len=100)
val.t <- otGreybody(nu=X$results[1,5], T=X$results[1,1], beta=beta,
     sc=X$results[1,4])
lines(3e5/nue.sweep, val.t*hfTail(nue.sweep/X$results[1,5], alpha=alpha),
      col=2, lwd=1, lty=2)
# data
wl <- c(250, 350, 500)
nue <- 3e5/wl*(1+z)
points(3e5/nue, l.nue, pch=16, col=3)
# }

Run the code above in your browser using DataLab