sedFitThin {cosmoFns}

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]

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.

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.

Note

Fit will sometimes crash on numerical derivative and throw an error. In this case the routine will halt without producing results. The more usual lack of convergence is reported as a warning, and the corresponding results will be NA in the output matrix.

Author(s)

A. Harris

References

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

Examples

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)