# This is an -script
# by R.E. Benestad, 14.06.2010
# R can be downloaded from http://cran.r-project.org 
# Fleagle & Businger (1980) 'An introduction to Atmospheric Physics'  eq. 5.18, p. 215
# to execute, write 'source("blackbody.R")' at the R-prompt.

rm(list=ls())


h<-6.64e-34
k<-1.38e-23
c<-3.00e8
alpha<-2*h*c^2

Ll <- function(Temp,l,alpha) alpha/( l^5 * (exp(h*c/(k*l*Temp))-1) )

i<-seq(log(100),log(20000),by=0.01)
l<-exp(i) * 1.00e-9

sun.bb<-Ll(6000,l,alpha)
earth.bb<-Ll(288,l,alpha)

T1 <- Ll(250,l,alpha)
T2 <- Ll(200,l,alpha)
T3 <- Ll(300,l,alpha)
T4 <- Ll(350,l,alpha)
T5 <- Ll(150,l,alpha)

map <- matrix(rep(NA,length(i)*10),length(i),10)
for (ii in 1:length(i)) map[ii,] <- ii

x <- l*1.00e9; y <- seq(1,10,length=10)
map[(x < 380) | (x > 750),] <- NA

spectrum <- rgb(seq(0,1,length=100),
                sin(pi*seq(0,1,length=100)),
                seq(1,0,length=100))

par(ps=12,bg="black",col.axis="white",col.main="white",col.lab="white")
plot(l*1.00e9,sun.bb,type="l",log="xy",lwd=3,
     main="Solar & Terestrial heatflux",
     xlab="wave length [nm]",col="yellow",
     ylab="Flux density",ylim=c(1,1e13))
image(x,y,map,col=spectrum,add=TRUE)
text(15,1.0e12,"Planck's law",cex=2,col="white")

lines(l*1.00e9,T1,lwd=1,col="grey")
lines(l*1.00e9,T2,lwd=1,col="grey")
lines(l*1.00e9,T3,lwd=1,col="grey")
lines(l*1.00e9,T4,lwd=1,col="grey")
lines(l*1.00e9,T5,lwd=1,col="grey")

lines(l*1.00e9,earth.bb,lwd=3,col="lightblue")

grid()
text(500,8.0e11,"6000K",col="yellow")
text(1e4,5.0e7,"288K",col="lightblue")

dev2bitmap(file="blackbody.png",res=150)
dev2bitmap("blackbody.pdf",type="pdfwrite")