#--Define blackbody function:
blackbody <- function(lambda, Temp=1e3) {
  h <- 6.626068e-34 ; c <- 3e8; kb <- 1.3806503e-23   # constants
  lambda <- lambda * 1e-6   # Convert from metres to microns
  ( 2*pi*h*c^2 ) / ( lambda^5*( exp( (h*c)/(lambda*kb*Temp) ) - 1 ) )
}

#--Define title and axis labels:
main <- "Planck blackbody curves"
xlab <- expression(paste(Wavelength, " (", mu, "m)"))
ylab <- expression(paste(Intensity, " ", (W/m^3)))

#--Line styles and colours:
## (see colorbrewer2.org)
## require(RColorBrewer)
## col <- brewer.pal(3, "Dark2")
col <- c("#1B9E77", "#D95F02", "#7570B3")
lty <- 1:3	# Line type
lwd <- 2.5	# Line width

#--Set up plot device:
pdf("blackbody.pdf", height=5.5, width=7, pointsize=14)

#--Plot curves:
curve(blackbody(lambda=x), from=1, to=15, main=main, xlab=xlab, ylab=ylab,
      col=col[1], lwd=lwd)
curve(blackbody(lambda=x, T=900), add=TRUE, col=col[2], lty=lty[2], lwd=lwd)
curve(blackbody(lambda=x, T=800), add=TRUE, col=col[3], lty=lty[3], lwd=lwd)

#--Add legend:
legtext <- paste(c(1000, 900, 800), "K", sep="")
legend("topright", inset=0.05, legend=legtext, lty=lty, col=col, 
       text.col=col, lwd=lwd, bty="n")