Planck distribution law for black body radiation

昔描きたくて色々やって結局諦めたようなことがあったのを思い出したので書いてみたけど、特に難しいことは無かった。
何に躓いていたのか思い出せないのは少しもったいない気分。

そーす

## Constants
h <- 6.6260693 * 10^-34                 # Planck's constant
c <- 299792458                          # Speed of light
kb <- 1.3806505*10^-23                  # Boltzmann constant 

## Planck distribution law for black body radiation
planck.dist <- function(t){
  function(lambda){
    lambda <- lambda * 10^-9
    2*h*c^2 / (lambda^5 * (exp(h*c/(lambda*kb*t)) - 1))
  }
}

## Wien displacement law
wien.disp.law <- function(t){
  h*c/(4.965*kb*t) * 10^9
}

draw.planck.dists <- function(temp){
  curve(planck.dist(temp[1])(x), 0, 2000,
        xlab="λ(nm)",
        ylab="",
        axes=F)
  axis(1)
  box()
  if(length(temp) > 1){
    for(i in 2:length(temp)){
      lines(planck.dist(temp[i])(0:2000), col=i)
    }
  }
  legend("topright", inset=0.05,
         legend = paste(temp, "K"),
         col = 1:length(temp),
         lty = 1,
         box.lty = 0,
         title = "Temperature")
  ## draw lambda max
  abline(v = wien.disp.law(temp),
         col = 1:length(temp),
         lty = 3)
  legend("right", inset = 0.05,
         legend = paste(round(wien.disp.law(temp)), "nm"),
         col = 1:length(temp),
         lty = 3,
         box.lty = 0,
         title = "λmax")
}

temp <- c(6000, 5000, 4000, 3000)
png(file = "Planck's_distribution_law.png")
draw.planck.dists(temp)
dev.off()