Как построить формулу с заданным диапазоном?

Глядя на страницу википедии уравнение спектральной яркости, кажется, что ваша формула немного неверна. Ваша формула умножает дополнительный pi (не уверен, что это задумано), и -1 находится внутри exp, а не снаружи:

L <- ((2*pi*h*c^2)/l^5)*((1/(exp((h*c)/(l*k*T)-1))))

Ниже приведена исправленная формула. Также обратите внимание, что я преобразовал его в функцию с параметром l, поскольку это переменная:

T <- 6000                # Absolute temperature
h <- 6.626070040*10^-34  # Plank's constant
c <- 2.99792458*10^8     # Speed of light in the medium
k <- 1.38064852*10^-23   # Boltzmann constant

L <- function(l){((2*h*c^2)/l^5)*((1/(exp((h*c)/(l*k*T))-1)))}

# Plotting
plot(L, xlim = c(20*10^-9,2000*10^-9),
     xlab = "Wavelength (nm)",
     ylab = bquote("Spectral Radiance" ~(KW*sr^-1*m^-2*nm^-1)),
     main = "Plank's Law",
     xaxt = "n", yaxt = "n")
xtick <- seq(20*10^-9, 2000*10^-9,by=220*10^-9)
ytick <- seq(0, 4*10^13,by=5*10^12)
axis(side=1, at=xtick, labels = (1*10^9)*seq(20*10^-9,2000*10^-9,by=220*10^-9))
axis(side=2, at=ytick, labels = (1*10^-12)*seq(0, 4*10^13,by=5*10^12))

Сюжет выше неплох, но я думаю, что с ggplot2 мы можем добиться большего:

h <- 6.626070040*10^-34  # Plank's constant
c <- 2.99792458*10^8     # Speed of light in the medium
k <- 1.38064852*10^-23   # Boltzmann constant

L2 <- function(l, T){((2*h*c^2)/l^5)*((1/(exp((h*c)/(l*k*T))-1)))} # Plank's Law
classical_L <- function(l, T){(2*c*k*T)/l^4} # Rayleigh-Jeans Law

library(ggplot2)
ggplot(data.frame(l = c(20*10^-9,2000*10^-9)), aes(l)) +
  geom_rect(aes(xmin=390*10^-9, xmax=700*10^-9, ymin=0, ymax=Inf), 
            alpha = 0.3, fill = "lightblue") +
  stat_function(fun=L2, color = "red", size = 1, args = list(T = 3000)) +
  stat_function(fun=L2, color = "green", size = 1, args = list(T = 4000)) +
  stat_function(fun=L2, color = "blue", size = 1, args = list(T = 5000)) +
  stat_function(fun=L2, color = "purple", size = 1, args = list(T = 6000)) +
  stat_function(fun=classical_L, color = "black", size = 1, args = list(T = 5000)) +
  theme_bw() +
  scale_x_continuous(breaks = seq(20*10^-9, 2000*10^-9,by=220*10^-9),
                     labels = (1*10^9)*seq(20*10^-9,2000*10^-9,by=220*10^-9),
                     sec.axis = dup_axis(labels = (1*10^6)*seq(20*10^-9,2000*10^-9,by=220*10^-9),
                                         name = "Wavelength (\U003BCm)")) +
  scale_y_continuous(breaks = seq(0, 4*10^13,by=5*10^12),
                     labels = (1*10^-12)*seq(0, 4*10^13,by=5*10^12),
                     limits = c(0, 3.5*10^13)) +
  labs(title = "Black Body Radiation described by Plank's Law", 
       x = "Wavelength (nm)",
       y = expression("Spectral Radiance" ~(kWsr^-1*m^-2*nm^-1)),
       caption = expression(''^'\U02020' ~'Spectral Radiance described by Rayleigh-Jeans Law, which demonstrates the ultraviolet catastrophe.')) +
  annotate("text", 
           x = c(640*10^-9, 640*10^-9, 640*10^-9, 640*10^-9, 
                 150*10^-9, (((700-390)/2)+390)*10^-9, 1340*10^-9), 
           y = c(2*10^12, 5*10^12, 14*10^12, 31*10^12, 
                 35*10^12, 35*10^12, 35*10^12), 
           label = c("3000 K", "4000 K", "5000 K", "6000 K", 
                    "UV", "VISIBLE", "INFRARED"),
           color = c(rep("black", 4), "purple", "blue", "red"),
           alpha = c(rep(1, 4), rep(0.6, 3)),
           size = 4.5) +
  annotate("text", x = 1350*10^-9, y = 23*10^12, 
           label = deparse(bquote("Classical theory (5000 K)"^"\U02020")),
           color = "black", parse = TRUE)

Заметки:

1. Я создал L2, также сделав абсолютную температуру T переменной
2. Для каждого T я строю график функции L2, используя разные цвета для представления. Я также добавил функцию classical_L, чтобы продемонстрировать классическую теорию спектральной яркости.
3. geom_rect создает светло-голубую заштрихованную область для "ВИДИМОГО" диапазона длин волн света.
4. scale_x_continuous устанавливает разрывы оси x, а labels устанавливает метки деления оси. Обратите внимание, я умножил seq на (1*10^9), чтобы преобразовать единицы в нанометры (нм). Вторая ось X добавлена ​​для отображения шкалы микрометра.
5. Аналогичным образом scale_y_continuous устанавливает разрывы и метки для оси y. Здесь я умножил на (1*10^-12) или (1*10^(-3-9)), чтобы преобразовать из ватт (Вт) в киловатты (кВт) и из обратного метра (м ^ -1) в обратный нанометр (нм ^ -1).
6. bquote правильно отображает верхние индексы в метке оси Y
7. annotate устанавливает координаты и текст меток кривых. Я также добавил метки для длин волн света "УФ", "ВИДИМОЙ" и "ИНФРАКРАСНЫЙ".

ggplot2

Сюжет из википедии:

Источник изображения: https://upload.wikimedia.org/wikipedia/commons/thumb/1/19/Black_body.svg/600px-Black_body.svg.png