Planck’s Black-Body Radiation (1900)
R(T,ν)dν=((2π h)/(c^2))((ν^3)/((e^(h ν/kT))-1))dν
R(T,λ)dλ=((2π h)/(c^2))((λ^5)/((e^(h c/λkT))-1))dλ

rad[T_,ν_]:=(ν^3)/(Exp[ν/T]-1)

Plot[{rad[300,x],rad[400,x],rad[500,x],rad[600,x]},{x,0,6000},PlotLegends->"Expressions"]

(Legended[, (LineLegend[rad(300,x), rad(400,x), rad(500,x), rad(600,x)])])

PlanckRadiationLaw[Quantity[300,"Kelvins"],Quantity[500,"Nanometers"]]

Plot[QuantityMagnitude[PlanckRadiationLaw[Quantity[6000,"Kelvins"],Quantity[l,"Nanometers"]]],{l,150,1500},AxesLabel->{"nm","energy density"}]

Stefan-Boltzmann Law : R = ((σT)^4)

Clear[h,k,c]

radSI[T_,ν_]:=(2 Pi h /(c^2))(ν^3)/(Exp[h ν/(k T)]-1)

(∫_0^∞)radSI[T,ν]\[DifferentialD]ν

Integrate[radSI[T,ν],{ν,0,∞},Assumptions->{ T > 0, h >0, k > 0}]

radnm[T_,λ_]:=(2 Pi h (c^2)/(λ^5))/(Exp[h c/(λ k T)]-1)

h=6.62608*^-34

k=1.38065*^-23

c=2.99792458*^8

Plot[{radSI[300,x],radSI[400,x],radSI[500,x],radSI[600,x]},{x,(10^8),(10^14)},PlotLegends->"Expressions",PlotRange->{{(10^8),(10^14)},{(10^(-12)),1.5×(10^(-10))}},AxesLabel->{Hz,W/(m^2)}]

radnm[300,0.0000005]

Plot[{radnm[3000,x×(10^(-9))],radnm[4000,x×(10^(-9))],radnm[5000,x×(10^(-9))],radnm[6000,x×(10^(-9))]},{x,100,1500},PlotLegends->{"3000 K","4000 K","5000 K","6000 K"},AxesLabel->{nm,W/(m^2)},PlotLabel -> "Planck Black-Body Radiation",  PlotRange->{{100,1500},{100,1.1×(10^14)}}]

Manipulate[
Plot[radSI[temp,x],{x,(10^8),1.5×(10^14)},PlotRange->{{(10^8),1.5×(10^14)},{(10^(-12)),6×(10^(-10))}},AxesLabel->{Hz,W/(m^2)},PlotLabel -> "Planck Black-Body Radiation"],
{temp,300,1000}]

Manipulate[
Plot[radnm[Temperature,x×(10^(-9))],{x,100,1600},AxesLabel->{nm,W/(m^2)},PlotLabel -> "Planck Black-Body Radiation",  PlotRange->{{100,1600},{1000,4.5×(10^14)}}],
{Temperature,3000,8000}]

Wein’s Displacement Law : (Subscript[λ, max]) T = ((h c)/(x k))

eqn=x+5 Exp[-x] == 5

FindRoot[eqn,{x,5}]

sol=FindRoot[eqn,{x,5},WorkingPrecision->10]

y=x /. sol

tw[ λmax_]:=h c / (λmax *(10^(-9)) y k)

tw[500]

lmax[T_]:=h c / (T y k)(10^9)

lmax[5778]

lmax[12000]

Einstein Theory of Heat Capacity of Solid

Clear[cv]

cv[x_,T_]:=3 R ((((x/T)(((Exp[x/(2 T)])/(Exp[x/T]-1)))))^2)

R=8.314

cv[350,100]

cv[350,1000]

cv[350,50]

Plot[cv[1000,T],{T,10,1000},PlotRange->{{0,1000},{0,25}}]

Manipulate[
Plot[cv[θ,T],{T,0,800},AxesLabel->{"T(K)","Cv (J/mol)"},PlotLabel -> "Einstein Capacity",  PlotRange->{{0,800},{0,25}}],
{θ,50,500}]