clear; clc; close all;

T1 = -10:150;
eps1 = 0.95;
lambdas = [8,14];
eps2 = 0.8:0.01:1;

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

colors = jet(length(eps2));

f=figure(); hold on;
for iE = 1:length(eps2)
    lambda = 1e-6*linspace(lambdas(1),lambdas(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
    
    % 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(iE);
    % from assumed intensity for epsilon=1 calculate temperature
    T2 = interp1(intensity,T0,int_eps2);

    plot(T1,T2-T1,'k');
    plot(T1,T2-T1,'--','Color',colors(iE,:));
    h=text(T1(end)+1,T2(end)-T1(end),sprintf('\\epsilon = %.2f',eps2(iE)));
    set(h,'Color',colors(iE,:));
end
xlabel(sprintf('Temperature measurement based on \\epsilon=%.2f',eps1))
ylabel('Temperature error')
grid on;
xlim([-10,172])
%legend(leg);




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