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

L = linspace(0,10000e-9,500);

Ts = linspace(100,6000,10);
colors = hot(length(Ts)+4);

f = figure('Position',[50,50,700,500]);
hold on;
for ii = 1:length(Ts)
    T = Ts(ii);
    I = c1./( L.^5 .* ( exp(c2./L/T) -1 ) );
    plot(L*1e6,I,'Color',colors(ii,:),'LineWidth',2);
    labels(ii,:) = char(sprintf('%04.0f °C',T-273.15));
end
xlabel('Wellenlänge / µm')
ylabel('Spektrale Intensität W/m²/nm')

hc = colorbar();
set(hc,'Position',[0.83 0.3 0.03 0.6]);
colormap(colors(1:end-2,:))
caxis([min(Ts)-mean(diff(Ts))/2,max(Ts)+mean(diff(Ts))/2])
set(hc,'ytick',Ts,'yticklabel',labels);

ax2 = axes('Position',[0.28,0.65,0.5,0.3]); hold on;
hold on;
for ii = 1:length(Ts)
    T = Ts(ii);
    I = c1./( L.^5 .* ( exp(c2./L/T) -1 ) );
    plot(L*1e6,I/max(I),'Color',colors(ii,:),'LineWidth',2);
end
xlabel('Wellenlänge / µm')
ylabel('normierte Intensität')

ax3 = axes('Position',[0.28,0.25,0.5,0.3]); hold on;
hold on;
for ii = 1:length(Ts)
    T = Ts(ii);
    I = c1./( L.^5 .* ( exp(c2./L/T) -1 ) );
    plot(L*1e6,I,'Color',colors(ii,:),'LineWidth',2);
end
xlabel('Wellenlänge / µm')
ylabel('Spektrale Intensität W/m²/nm')
ylim([0,60])



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