echo on
clc

% -------------------- black-body.m -------------------------
%
% Lecture 1:   The Planck Blackbody Radiation Law
% 			   Numerical Calculations with MATLAB
%
%				Chemistry 455
%
%				James B. Callis
%				University of Washington
%				Summer 1997
%
%
%	Press any key to continue
pause
clc

% 				What is a Blackbody?
%
% Answer:
%
% It is an idealized source of electromagnetic radiation that 
% is produced by heating matter.
%
% Blackbodies are generally hollow spheres with a perfectly 
% absorbing interior, maintained at temperature T.
%
% At thermal equilibrium, as much energy that is emitted 
% by the body is absorbed by the body.
%
% A small hole is drilled in the blackbody so that a sample 
% of energy can be obtained.
%
%	Press any key to continue
pause
clc

%		What is the Nature of the Distribution of Radiation
%		Emitted by a Blackbody?
%
% The distribution of frequencies is measured in terms of energy 
% density (which has the mks units of Joules/meter^3) as 
% a function of frequency (units of sec-1).
%
% The spectrum of radiation emitted by the body is solely a 
% function of the temperature of the body. 
%
% As the temperature of the blackbody is raised, the spectrum of 
% light shifts toward the blue (higher frequency). (See Figure
% 1-1 of McQuarrie. 
%
%	Press any key to continue
pause
clc

% Classical Physics Cannot Explain Blackbody Radiation
%
% Rayleigh-Jeans Law:
%
% p(v,T)*dv = [(8*pi*k.*T/c^3).*v.^2]*dv
%				(note the use of MATLAB notation)
%
% where p(v,T) is the density of radiation energy with
% frequencies between v and v + dv. k is the Boltzmann 
% constant, T is the temperature in Kelvin and c is the
% speed of light. 
%
% Note that the spectrum increases constantly as v increases.
%
% Note also that the total amount of energy emitted by such a 
% body is infinite.
%	Press any key to continue
pause
clc

% MATLAB can be extremely useful in understanding functional
% relationships like the Rayleigh-Jeans Law. We can visualize
% how the dependent variable, p, is influenced by variation in
% independent variable v by making a plot like that of 
% figure 1-1 of McQuarrie. 
%
% To make a plot in MATLAB, we carry out three steps:
%
%			1 - Define the constants.
%			2 - Calculate the values of the dependent and
%				independent variables.
%			3 - Make the plot with titles and labels.
%
%
%	Press any key to continue
pause
clc

 
% Step 1 - Define the constants:
%
%
h = 6.626e-34; % Planck's constant in J-sec
c = 3.00e8; % Speed of light in a vacuum (m/sec)
T= 5000; %(Temperature in degrees K)
k = 1.381e-23; % Boltzmann constant in J-K-1
%
%
%
%
%	Press any key to continue
pause
clc

% Step 2 - Calculate the arrays of values of independent
% and dependent values:

% First calculate the values of the independent variable:
v = .015e14:.015e14:15e14;% Do the entire range of frequencies
% Now calculate the values of the radiation density according
% to the R-J law. Note the use of dot multiplication notation.
p = (8*pi*k.*T/c^3).*v.^2; % The Rayleigh-Jeans law
%
%
%
%	Press any key to continue
pause
clc
% Step 3 - Make a plot:


plot(v,p)
xlabel('Frequency, sec-1')
ylabel('Energy density, J-m-3')
title('Rayleigh-Jeans Radiation Law')
pause
clg
clc
% Planck's Law
%
% In 1900, Planck proposed the following empirical law, which
% satisfactorily duplicated the observed blackbody radiation:
%
% p_v*dv = [(8*pi*h*v.^3/c^3).*(1./(exp(h*v./(k*T))-1))]*dv;
%
% Where h is an empirical constant whose value, Planck adjusted
% to agree with experiment. 
%
% Much to Planck's delight, h was a constant for all blackbodies
% at all temperatures.
%
% h = 6.626e(-34) joule-sec (note MATLAB notation)
%
%	Press any key to continue
pause
clc

% Again, MATLAB can be very helpful in understanding functions
% like Planck's Blackbody Radiation Law.
%
% One strategy to explore complex functions is to break the 
% equation into its parts. We see that the Planck law consists 
% of two parts:
%
% Part a: (8*pi*h*v.^3)/c^3
%
% Part b: (1./(exp(h*v./(k*T))-1))
%
% We could plot these two parts together with the complete
% expression to see how the function gets its unique form.
%
%	Press any key to continue
pause
clc

% Here is the MATLAB code:
%
p_a = (8*pi*h*v.^3)/c^3.; % This generates the first part
p_a_norm = p_a/max(p_a); % This normalizes the first part
p_b = (1./(exp(h*v./(k*T))-1)); % This generates the second
p_b_norm = p_b/max(p_b); % This normalizes the second part
p_p = p_a.*p_b; %Planck's law is just the product of the two
p_p_norm = p_p/max(p_p); % This normalizes the second part
%
plot(v,p_p_norm, v, p_a_norm, v,p_b_norm)
xlabel('Frequency, sec-1')
ylabel('Energy density, J-m-3')
title('Planck Radiation Law and Its Components')
pause
clg
clc

% Now we would like to compare the Planck Radiation law
% and the Rayleigh-Jeans law as in Figure 1-1 of McQuarrie
% It is highly desirable to plot the energy density for
% four temperatures on one plot. MATLAB can do this conveniently
% by use of the function meshgrid. [X,Y] = MESHGRID(x,y) 
% transforms the domain specified by vectors x and y into 
% arrays X and Y that can be used for the evaluation
% of functions of two variables and 3-D surface plots.
%
%	Press any key to continue
pause
clc
% Here we generate the R-J Law curve as per Figure 1-1


v_rj = .015e14:.015e14:1.5e14; % R-J law calculation is only over
% a limited range of frequencies.
p_rj = (8*pi*k.*6000/c^3).*v_rj.^2; % The Rayleigh-Jeans law for 
% T = 6000 K. 
%
%	Press any key to continue
pause
clc
 	  
% Now Do the Planck Law for 4 Temperatures
% Here is the temperature vector:
%
temp_vec = 3000:1000:6000;
%
% Now generate a pair of matrices representing a rectangular
% grid of points in v and T (wavelength and temperature).

v = .015e14:.015e14:15e14;% Do the entire range of frequencies
[V,T] = meshgrid(v,temp_vec);
p_nu = (8*pi*h*V.^3/c^3).*(1./(exp(h*V./(k*T))-1));


%	Press any key to continue
pause
clc

plot(v_rj,p_rj,v,p_nu)
title('Comparison of Planck and Rayleigh-Jeans Laws')
ylabel('Energy Density, J/m^3')
xlabel('Frequency, Hz')
legend('R-J Law (6000 K)', 'T = 3000 K', 'T = 4000 K', 'T = 5000 K', 'T= 6000 K')
pause
clg

clc
% Having verified the excellent agreement of his equation with
% experiment, Planck set out to find an explanation. He found 
% that he could derive his expression from classical physics
% (statistical thermodynamics) by making the assumption that 
% radiation was being emitted and absorbed by oscillators in the 
% wall, each with a characteristic frequency, v. 
% 
% However, he made the revolutionary assumption that the energy
% stored in the oscillators could only take on discrete values:
%
% e = n*h*v,  n = 0,1,2,3....
%
%	Press any key to continue
pause
clc
% Triumphs of the Planck Law
%
% 1. The discovery of h, one of the universal constants of nature 
%
% 2. The interpretation of the constant of the Wein displacement
%    Law in terms of fundamental physical constants (h, c and k).
%
% 3. The interpretation of the constant, sigma of the Stefan-
%    Boltzmann law in terms of the fundamental physical constants
%    (h,c and k).
%
% End of MATLAB script for Lecture 1.