Contents

RADIOMETRY plots various black-body calculations

by Chuck DiMarzio Northeastern University December 2008

!! This file may be copied, used, or modified for educational and
!! research purposes provided that this header information is not
!! removed or altered, and provided that the book is cited in
!! publications, as DiMarzio, Charles A., Optics for Engineers,
!! CRC Press, Boca Raton, FL, 2011.
!! http://www.crcpress.com
!! Other distribution is prohibited without permission.
Requires BBODY, BBSPEC, and data files sunspecexo.dat sunspecsl.dat
[t,lambda]=meshgrid(fliplr([100,300,500,1000,3000,5000,10000]),...
                    (0.01:0.01:100));
m_lambda=bbspec(t,lambda);
%

Plot black body curves

[mpoint,test]=max(m_lambda,[],1);
lpoint=lambda(test,:);
symb=['o';'x';'+';'^';'>';'<';'v'];
set(0,'DefaultAxesColorOrder',[0 0 0]);
fig1=figure;
for n=1:length(lpoint);
  loglog(lpoint(n),mpoint(n),symb(n));hold on;
end;
legend(num2str(t(1,:)'),'Location','NorthEast');
loglog(lambda,m_lambda,lpoint,mpoint);axis([0,100,1e-2,1e10]);
xlabel('\lambda, Wavelength, \mu m');
ylabel('m_\lambda, Spect. Radiant Exitance, W/m^2/\mu m');
grid on;
hold off;

Outdoor radiant exitance levels

tsun=5000;dlambda=0.01;      % effective solar temperature
lambda1=(dlambda:dlambda:100);  % Set wavelength for sun spectra
[msun,fracsun,sremax,lpeak]=bbody(tsun,0.1,100);
% msun is M, radiant exitance, fracsun is fraction in the band,
% sremax is the mzximum spectral radiant exitance, at wavelength lpeak
msun_lambda=bbspec(tsun,lambda1);
mcloud_lambda=1000*msun_lambda/msun;  % Nominal solar irradiance is 1000w/m^2
%                             Assume 100% diffuse reflectance
moose=msun_lambda./lambda1.^4;
msky_lambda=27*pi*moose/sum(moose)/dlambda;
mnight_lambda=bbspec(300,lambda1);

fig2=figure;
testsun=find(msun_lambda==max(msun_lambda));
testcloud=find(mcloud_lambda==max(mcloud_lambda));
testsky=find(msky_lambda==max(msky_lambda));
testnight=find(mnight_lambda==max(mnight_lambda));
loglog(lambda1(testsun),msun_lambda(testsun),'o',...
     lambda1(testcloud),mcloud_lambda(testcloud),'x',...
     lambda1(testsky),msky_lambda(testsky),'+',...
     lambda1(testnight),mnight_lambda(testnight),'*');hold on;
legend('Sun','Cloud','Blue Sky','Night',2);
loglog(lambda1,msun_lambda,lambda1,mcloud_lambda,...
       lambda1,msky_lambda,lambda1,mnight_lambda)
xlabel('\lambda, Wavelength, \mu m');
ylabel('m_\lambda, Spect. Radiant Exitance, W/m^2/\mu m');
axis([0,100,1e-2,1e10]);
grid on;
hold off;

Plot actual solar spectra

load astm;
lambdaastm=data(:,1);
etr=data(:,2);
gt=data(:,3);
dpcc=data(:,4);
lambda2=lambdaastm/1000;  % Raw data is in nm, plot in microns
e_sun_exo_lambda=etr*1000;     % Spectral Irradiance convert from
                               % w/m^2/sr/nm to w/m^2/sr/\mu m
e_sun_sl_lambda=gt*1000;     % spectral Irradiance at surface
[m6k,fracsun,sremax,lpeak]=bbody(6000,0.1,100);
% m6k is M, radiant exitance, fracsun is fraction in the band,
% sremax is the mzximum spectral radiant exitance, at wavelength lpeak
e6k_lambda=bbspec(6000,lambda2)/m6k*1480;
[m5k,fracsun,sremax,lpeak]=bbody(5000,0.1,100);
% m5k is M, radiant exitance, fracsun is fraction in the band,
% sremax is the mzximum spectral radiant exitance, at wavelength lpeak
e5k_lambda=bbspec(5000,lambda2)/m5k*1250;
fig3=figure;plot(lambda2,e_sun_exo_lambda,'--',...
            lambda2,e_sun_sl_lambda,'--',...
            lambda2,e6k_lambda,'-',...
            lambda2,e5k_lambda,'-');
grid on;
xlabel('\lambda, Wavelength, \mu m');
ylabel('e_\lambda, Spect. Irradiance, W/m^2/\mu m');

solar radiance = integral of spectral radiance d lambda

intetr=trapni(etr,lambdaastm)
intgt=trapni(gt,lambdaastm)
intetr =
   1.3480e+03
intgt =
   1.0005e+03

Total radiant and luminous exitance vs. Temperature

dlambda=0.005;lambda_axis=(dlambda:dlambda:50);
t_axis=(600:50:10000);
[t,lambda]=meshgrid(t_axis,...
                    lambda_axis);
L_lambda=bbspec(t,lambda)/pi;
v=ones(size(lambda,1),1)*dlambda;
L=L_lambda'*v;
[x,y,z]=xyz(lambda_axis'*1000);  % xyz wants nanometers for lambda
LL=L_lambda'*y*683*dlambda;

fig4=figure;
axes('position',[0.2,0.2,0.6,0.6]);
[ax,h1,h2]=plotyy(t_axis,L,t_axis,LL,'semilogy');
grid on;
xlabel('T, Temperature, K');
ylabel(ax(1),'M, Radiance, W/m^2/sr');
ylabel(ax(2),'M, Luminance, L/m^2/sr');
axis(ax(1),[0,10000,1e-10,1e12]);
axis(ax(2),[0,10000,1e-10,1e12]);

hold on;plot([100,1500],1.1e7*[1,1],[8500,9900],7e8*[1,1]);
text(100,3e7,'Solar Disk');text(8200,2e9,'Solar Disk');

plot([1100,2300],7e-10*[1,1],[8600,9900],5e-7*[1,1]);
text(1100,7e-9,'Min. Vis.');text(8200,5e-6,'Min. Vis.');

plot([1200,2700],13*[1,1],[8300,9900],2500*[1,1]);
text(1500,130,'Lunar Disk');text(8100,25000,'Lunar Disk');

plot([5000,6000],LL(find(t_axis==5000))*[1,1]);
plot(6000,LL(find(t_axis==5000)),'>');
plot([5000,4000],L(find(t_axis==5000))*[1,1]);
plot(4000,L(find(t_axis==5000)),'<');

hold off;

fig5=figure;plot(t_axis,LL./L);
grid on;
xlabel('T, Temperature, K');
ylabel('Efficiency, lm/W');

% To get the colors right for low temperatures, the x,y,z fit isn't good
% enough.  Load the actual data, and use that
load xyzdata;% a(1) is wavelength, a(2) is xbar, b(2) is ybar, c(2) is
                % zbar
lambdaa_axis=a(:,1);dlambdaa=5;
% t_axis=[600:50:10000];  % Already defined
[t,lambdaa]=meshgrid(t_axis,...
                    lambdaa_axis);
L_lambdaa=bbspec(t,lambdaa)/pi;


XX=L_lambdaa'*a(:,2)*683*dlambdaa;
YY=L_lambdaa'*b(:,2)*683*dlambdaa;
ZZ=L_lambdaa'*c(:,2)*683*dlambdaa;


xhot=XX./(XX+YY+ZZ);
yhot=YY./(XX+YY+ZZ);

fig6=figure;plot(t_axis,xhot,'-',t_axis,yhot,'-.');
legend('x','y');
grid on;
xlabel('T, Temperature, K');
ylabel('x,y, Color Coordinates');
save hotcolors t_axis xhot yhot

fivek=find(t_axis==5000)
radianceofsolardisk=L(fivek)  % Radiance of solar disk
luminanceofsolardisk=YY(fivek) % Luminance of solar disk
fivek =
    89
radianceofsolardisk =
   1.1292e+07
luminanceofsolardisk =
   9.1929e+11

Replot of the solar spectra on separate plots

fig31=figure;plot(lambda2,e_sun_exo_lambda,':',...
            lambda2,e6k_lambda,'-');
grid on;
xlabel('\lambda, Wavelength, \mu m');
ylabel('e_\lambda, W/m^2/\mu m');

fig32=figure;plot(lambda2,e_sun_sl_lambda,':',...
            lambda2,e5k_lambda,'-');
grid on;
xlabel('\lambda, Wavelength, \mu m');
ylabel('e_\lambda, W/m^2/\mu m');