# 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