    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %    Calculate the Radiation       %
    %    Curves for a blackbody        %
    %                                  %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%   Define the wavelength
%%%%%%%%%   Range and increment.

T = 300        %  sample value of temperature 
lambdamin = 0.1;    %  Minimum wavelength (microns)
lambdamax = 100.;   %  Maximum wavelength (microns)

%  T = 6000        %  sample value of temperature 
%  lambdamin = 0.1;    %  Minimum wavelength (microns)
%  lambdamax = 5.;   %  Maximum wavelength (microns)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NN = 251;  %  Number of values to be calculated.

lambdamin = lambdamin/10^6  %  Change to metres
lambdamax = lambdamax/10^6  %  Change to metres
lambdadiff = (lambdamax-lambdamin)/(NN-1); % Wavelength increment (m).

%%%%%%%%  Define a few constants
c = 3*10^8;             %   Speed of light
h = 6.626*10^(-34);     %   Planck's constant
k = 1.38*10^(-23);      %   Boltzmann's constant

c1 = 2*pi*h*c^2;        %   First constant in Planck formula 
c2 = h*c/k;        %  Second constant in Planck formula 


%%%%%%%%   Loop through the wavelengths.

for N = 1:NN
  lambda = lambdamin + (N-1)*lambdadiff;
  LAMBDA(N) = lambda;
  B(N)  = c1*lambda^(-5) / ( exp(c2/(lambda*T)) -1) ; % exact
  Bx(N) = c1*lambda^(-5) / ( exp(c2/(lambda*T))   ) ; % approx
  error(N) = (B(N)-Bx(N));
end

%%%%%%%%   Plot the results
plot(LAMBDA,B,'b-',LAMBDA,Bx,'r:',LAMBDA,error,'m.')

 title('Balckbody Radiation')
 xlabel('Wavelength')
 ylabel('B')
 
 legend('Exact','Approx','Error')
 

return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
