% function handles to integrate the blackbody emittance with respect to
% wavenumber, and then compare to the stefan-boltzmann result.
bbemit = @(T,nu) 1.191e-8*pi*nu.^3 ./ (exp(1.439*nu./T)-1);
stefbolt = @(T) 5.67e-8 * T.^4;

T = 255;

stefbolt(T)

bbemit_T = @(nu) bbemit(T,nu);
quad(bbemit_T, 0, 8000)