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