昔描きたくて色々やって結局諦めたようなことがあったのを思い出したので書いてみたけど、特に難しいことは無かった。
何に躓いていたのか思い出せないのは少しもったいない気分。
そーす
## 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()