clear; clc; close all;

T1 = 5:2000;
eps1 = 0.95;
eps2 = 0.85;
lambdas = [8,14;4.5,5.5;2,2.8;1.45,1.8;0.7,1.1];
%lambdas = [8,14];

%%
TK = 273.15;
T0 = 0:2500;
kB = 1.38064852e-23;
n = 1;
h = 6.626075e-34;
c = 299792458;
F = h*c/kB/n;


f=figure(); hold on;

for iL = 1:size(lambdas,1)
    lambda = 1e-6*linspace(lambdas(iL,1),lambdas(iL,2),10000);
    clear intensity tmp int_eps1 int_eps2 int_0 T2;
    for iT = 1:length(T0)
        %tmp = Planck(lambda,T0(iT)+TK);
        Kelvin = T0(iT)+TK;
        tmp = 1./(lambda.^5.*(exp(F./(lambda.*Kelvin))-1));
        intensity(iT) = sum(tmp);
    end
    
    if false
        figure(); hold on;
        plot(T0,intensity);
        xlabel('Temperature of object, \epsilon=1  [°C]')
        ylabel(sprintf('radiated intensity [%.0f,%.0f] µm',min(lambda)*1e-3,max(lambda)*1e-3))
    end
    
    % form temperature measurement, assume what epsilon=1 intensity the meter
    % assumed
    int_eps1 = interp1(T0,intensity,T1);
    % from assumed temperature calculate measured intensity by inverting the
    % assumption epsilon=1 from epsilon=eps1
    int_0 = int_eps1*eps1;
    % from measured intensity calculate assumed temperature if a different
    % epsilon is assumed
    int_eps2 = int_0/eps2;
    % from assumed intensity for epsilon=1 calculate temperature
    T2 = interp1(intensity,T0,int_eps2);

    plot(T1,T2-T1);
    
    leg{iL} = sprintf('\\lambda = %.1f .. %.1f µm',lambdas(iL,1),lambdas(iL,2));
end
xlabel(sprintf('Temperature measurement based on \\epsilon=%.2f',eps1))
ylabel(sprintf('Temperature error T(\\epsilon=%.2f) - T(\\epsilon=%.2f)',eps2,eps1))
grid on;
legend(leg);




print(f,'IR_Thermometer_epsilon_effect1','-dpng')