The Planck Law is really easily derived. In the following I use natural units with [itex]\hbar=c=k_{\text{B}}=1[/itex].
You start from the quantization of the free electromagnetic field (most simply in the radiation gauge, [itex]A^0=0[/itex], [itex]\vec{\nabla} \cdot \vec{A}=0[/itex]). The Fock basis is given by
[tex]|\{N(\vec{p},\lambda) \} \rangle,[/tex]
where I use a cube of length [itex]L[/itex] as a quantization volume with periodic boundary conditions. Then [itex]\vec{p} \in \frac{2 \pi}{L} \mathbb{Z}^3[/itex] and [itex]\lambda \in \{-1,1\}[/itex] (helicities).
Since photons have spin 1, they are bosons and thus the occupation numbers all run from 0 to [itex]\infty[/itex]. The canonical partition sum for photons of fixed energy [itex]E=|\vec{p}|[/itex] thus is
[tex]Z(\beta,E)=\sum_{N=0}^{\infty} \exp(-2 \beta E N)=\left (\frac{1}{1-\exp(-\beta E)} \right )^2,[/tex]
where the square comes from the two helicities (polarizations) of each photon mode.
The average total energy of the photons of energy [itex]E[/itex] is thus
[tex]\langle E_{\text{tot}} \rangle_{E}=-\partial_{\beta} \ln Z(\beta,E)=\frac{2 E}{\exp(\beta E)-1}.[/tex]
To get the energy spectrum, we have to count the states. In a momentum-volume element [itex]\Delta^3 \vec{p}[/itex] we have [itex]\mathrm{\Delta}^3 \vec{p} L^3/(2 \pi)^3[/itex] states. Thus the energy-density spectrum is in the limit [itex]L \rightarrow \infty[/itex]
[tex]\mathrm{d} u(E)=\frac{\mathrm{d}^3 \vec{p}}{(2 \pi)^3} \frac{2 |\vec{p}|}{\exp(\beta |\vec{p}|)-1}.[/tex]
Because of [itex]E^2=\vec{p}^2[/itex] we have
[tex]\mathrm{d}^3 \vec{p} = \vec{p}^2 \mathrm{d} |\vec{p}| \mathrm{d} \Omega= E^2 \mathrm{d} E \mathrm{d} \Omega.[/tex]
Integrating over the full solid angle gives
[tex]\mathrm{d} u=\mathrm{d} E \frac{8 \pi}{(2 \pi)^3} \frac{E^3}{\exp(\beta E)-1}.[/tex]
This is Planck's Law:
[tex]u(E)=frac{8 \pi}{(2 \pi)^3} \frac{E^3}{\exp(\beta E)-1}.[/tex]
Usually it's written in terms of the frequency [itex]\nu=E/(2 \pi)[/itex] which gives
[tex]u(\nu)=2 \pi u(E)=\frac{16 \pi^2 \nu^3}{\exp(2 \pi \beta \nu)-1}[/tex]
or reinstalling all factors [itex]h=2 \pi \hbar[/itex], [itex]c[/itex] and [itex]k_{\text{B}}[/itex]
[tex]u(\nu)=\frac{8 \pi h \nu^3}{c^3} \frac{1}{\exp[h \nu/(k_{\text{B}} T)]-1}.[/tex]
See Wikipedia for a thorough further discussion:
https://en.wikipedia.org/wiki/Planck's_law