clear all;

hbar = 1.0546e-34;

clight = 3e8;

echarge = 1.6e-19;

kb = 1.38e-23/1.6e-19; %kb in eV



L = 1; %length of cube in microns

L = L*1e-6;



Tpho = logspace(2.477, 3.477, 100);

alpha = hbar*clight*pi/L;

max = 10;



for Tindex = 1:size(Tpho, 2)

energy(Tindex) = 0;

T = Tpho(Tindex);

for l = 1:max

for m = 1:max

for n = 1:max

wavee = alpha*(l^2+m^2+n^2)^(0.5)/echarge;

energy(Tindex) = energy(Tindex) + wavee/(exp(wavee/(kb*T))-1);

end

end

end

end



energy = 2*energy;



energydensity = energy/L^3; %in eV/m^3



energydensity = energydensity*echarge;



semilogy(Tpho, energydensity, 'b')

hold

semilogy(Tpho, (4*5.67e-8*Tpho.^4/clight), 'r')