# 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")