function bfract(L,T)
% blackbody fractional function
%
% syntax : f = bfract(L,T);
% L : wavelength (m)
% T : temperature (K)
%
% f is the portion of radiation emitted by a blackbody 
% at temperature T in the wavelength interval [0,L]
%

disp('---------------------------------------------------')
disp([' wavelength (m)  : ',sprintf('%0.5e',L)])
disp([' temperature (K) : ',num2str(T)])
disp('---------------------------------------------------')
disp([' blackbody fractional function (-) : ',sprintf('%0.5e',fract(L*T))])
disp('---------------------------------------------------')


function c = fract(a)
%
% Function fract(a) calculates the blackbody fractional function
% F_{0-a}, where a is the wavelength-temperature product (m*K).
%
C2 = 1.438775225e-2;
zeta = C2/a;
n = 10;
b = 0.;
for i = 1:n
    b = b + exp(-i*zeta)/i*(zeta^3+3.*zeta^2/i+6.*zeta/i^2+6./i^3);
end
c = b*15.0/pi^4;

