clear; clc; close all;
h=6.626e-34;
c=3e8;
k=1.3806e-23;
c1 = 1e-7*2*pi*h*c^2;
c2 = h*c/k;

L = 4.05:0.1:17.95;
L = L*1e-6;

Ts = 273.15+[-20:20:120,200,300,400];
colors = hot(length(Ts)+4);

f = figure('Position',[50,50,600,450]);
hold on;
for ii = length(Ts):-1:1
    T = Ts(ii);
    I = c1./( L.^5 .* ( exp(c2./L/T) -1 ) );
    plot(L*1e6,I,'Color',colors(ii,:),'LineWidth',2);
    Int(ii) = 1e6*sum(I(L>8e-6&L<14e-6))*mean(diff(L));
    labels{length(Ts)+1-ii} = sprintf('% 3.0f °C    intensity % 3.0f µW/m²',T-273.15,Int(ii));
end
M = 30;
plot(8*[1,1],[0,M],'k-','Linewidth',2)
plot(14*[1,1],[0,M],'k-','Linewidth',2)
ylim([0,M])
xlim([4,18])
xlabel('Wellenlänge / µm')
ylabel('Spektrale Intensität W/m²/nm')
legend(labels)
print(f,'IR_8_14_a.png','-dpng','-r150');

f = figure('Position',[50,50,380,320]);
hold on;
plot(Ts-273.15,Int,'k-');
h = scatter(Ts-273.15,Int,20,Int);
set(h,'MarkerFaceColor',get(h,'MarkerEdgeColor'));
colormap(colors);
xlim([-25,410]);
xlabel('Temperature / °C');
ylabel('Intensity between 8 µm and 14 µm  / µW/cm²')

print(f,'IR_8_14_b.png','-dpng','-r150');