# Optically-thin SED fit 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. ```r 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. ## Returns 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`. ## References Hogg 2000, astro-ph 9905116v4; Blain, Barnard & Chapman 2003, MNRAS 338, 733. ## Author(s) A. Harris ## 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. ## Examples ```r 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) ```

sedFitThin function