function [wavelr,mrlambda,wavelv,mvlambda,wavelp,mplambda,m,mv,mp...
         X,Y,Z]=...
    bbd(t,lambda1,lambda2);
%
%  bbd(t,lambda1,lambda2);
%  Calls bbody and displays results
%  for a black body at temp t, kelvin
%  over wavelengths lambda1 to lambda2, microns
%  [wavelr,mrlambda,wavelv,mvlambda,wavelp,mplambda,m,mv,mp,X,Y,Z]=...
%    bbd(t,lambda1,lambda2);
%
%  Plots and 
%  returns wavelr,mrlambda 
%          (wavelength and spectral radiant exitance, W/m^2/micron)
%
%  wavelv,mvlambda 
%          (wavelength and spectral luminous exitance, lumens/m^2/micron)
%
%  wavelp,mplambda 
%          (wavelength and spectral photon exitance, photons/sec/m^2/micron)
%
%  Also returns
%  m   radiant exitance, W/m^2/micron)
%  mv  luminous exitance, lumens/m^2/micron)
%  mp  photon exitance, photons/sec/m^2/micron)
%
% and
% X,Y,Z tristimulus values
%
%
%  By Chuck DiMarzio
%  Northeastern University
%  2006, 2009

disp(['---------']);
disp(['Temperature = ',num2str(t),' Kelvin']);
disp(['Wavelengths from ',num2str(lambda1),...
        ' to ',num2str(lambda2),' micrometer']);
disp([' ']);
[m,frac,sremax,lpeak]=bbody(t,lambda1,lambda2);
disp(['Peak Spectral Radiant Exitance, m_\lambda = ',num2str(sremax,4),...
     ' Watts/m^2/micron at \lambda = ', num2str(lpeak,4),' microns']);
disp(['Radiant Exitance, M = ',num2str(m,4),' Watts/m^2 = ',...
        num2str(frac,4),' of total in 0 to 100 microns']);
disp(['---------']);
dwavelr=min(0.005,abs(lambda2-lambda1)/20);
wavelr=[lambda1:dwavelr:lambda2];
mrlambda=bbspec(t,wavelr);  % Spectral radiant exitance, W/m^2
c=299792458; % m/s
h=6.626068E-34; % Joules/Hz =  m2 kg / s
mplambda=mrlambda.*(wavelr*1e-6)/h/c;  % Spectral photon exitance,
                                       % Photons/sec/m^2 
figure;plot(wavelr,mrlambda);
xlabel('\lambda, Wavelength, \mu m');
ylabel('Mr_\lambda, W/m^2/\mu m');
grid on;

data=load('~/Documents/data/x.ind');
wavel=data(:,1);x=data(:,2);
data=load('~/Documents/data/y.ind');
wavel=data(:,1);y=data(:,2);
data=load('~/Documents/data/z.ind');
wavel=data(:,1);z=data(:,2);
dwavel=wavel(2)-wavel(1);
test=find(wavel>lambda1 & wavel < lambda2);
if(length(test)==0);test=find(wavel<lambda1,1,'last');end;
if(length(test)==0);test=find(wavel>lambda1,1,'first');end;
wavelv=wavel(test);  % Wavelengths to plot
mlambda=bbspec(t,wavelv); 
mvlambda=mlambda.*y(test)*683;  % spectral luminous exitance lumens/m^2/micron
integrand=x(test).*mlambda;X=sum(integrand)*dwavel;
integrand=y(test).*mlambda;Y=sum(integrand)*dwavel;
integrand=z(test).*mlambda;Z=sum(integrand)*dwavel;
xx=X/(X+Y+Z);
yy=Y/(X+Y+Z);
zz=Z/(X+Y+Z);
mv=Y*683;  % Luminous exitance, l/m^2
disp(['Luminous Exitance, Mv = ',num2str(mv,4),' l/m^2']);
disp(['Chromaticity, x = ',num2str(xx)]);
disp(['              y = ',num2str(yy)]);
disp(['              z = ',num2str(zz)]);
figure;plot(wavelv,mvlambda);
xlabel('\lambda, Wavelength, \mu m');
ylabel('Mv_\lambda,  lm/m^2/\mu m');
grid on;

disp(['Luminous Efficiency, ', num2str(mv/m,4),' lumens/W']);
disp(['---------']);
mp=sum(mplambda*dwavelr);  % Photon exitance, Photons/m^2
disp(['Photon Exitance, Mp =',num2str(mp,4),' Photons/sec/m^2']);
wavelp=wavelr;
figure;plot(wavelr,mplambda);
xlabel('\lambda, Wavelength, \mu m');
ylabel('Mp_\lambda, Photons/sec/m^2/\mu m');
grid on;
disp(['Photon Efficiency, ', num2str(mp/m,4),' Photons/sec/W']);
disp(['                   ', num2str(mp/mv,4),' Photons/sec/lumen']);