%
%       lab2_answers.m
%
%   ....This script reproduces the data analysis and plots from Lab 2.
%
%       Written by Von P. Walden
%                  10 February 1998
%

%
% Part A
%
wv = 0.1:0.01:200; 
T = [300 1000 6000]; 
Bwv = planck_wv(T,wv); 
loglog(wv,Bwv) 
axis([.1 200 0.01 1e8]) 
xlabel('Wavelength (micrometers)'); 
ylabel('Radiance (W m-2 um-1 str-1)') 
title('Planck curves in Figure 1.3 (Part A)') 

%
% Part B
%
for i = 1:3 
    wv_max(i) = wv( find(Bwv(i,:) == max(Bwv(i,:)) )  ); 
end 
disp('Wavelengths of maximum emission (Part B): ')
wv_max

figure
plot(1./T,wv_max,'*',1./T,wv_max,'-') 
xlabel('1 / T (K -1)')
ylabel('Wavelength (micrometers)')
title('Example of Wien''s Law (Part B)')

%   ....Slope and intercept
p = polyfit(1./T,wv_max,1);
disp('Slope and intercept of Wien''s Law plot (Part B):')
slope = p(1)
intercept = p(2)

%
% Part C
%
wn = 1:1:20000; 
Bwn = planck_wn(T,wn); 

figure
loglog(wn,Bwn) 
axis([1 20000 0.01 1e8]) 
xlabel('Wavenumber (cm-1)'); 
ylabel('Radiance [W m-2 str-1 (cm-1)-1]') 
title('Planck curves vs. wavenumber (Part C)') 

for i = 1:3 
    wn_max(i) = wn( find(Bwn(i,:) == max(Bwn(i,:)) )  ); 
end 
disp('Wavenumbers of maximum emission (Part B): ')
wn_max


%
% Part D
%
disp('Intergrated flux values using TRAPZ: (Part D)')
Flux_wv = pi*trapz(wv,Bwv(1,:)) 
Flux_wn = pi*trapz(wn,Bwn(1,:))/1000 
disp('Flux from Stefan-Boltzmann Law: (Part D)')
5.67e-8*300.^4


%
% Part E
%
nBwv_6000 = Bwv(3,:)./max(Bwv(3,:)); 
nBwv_300  = Bwv(1,:)./max(Bwv(1,:)); 

figure
semilogx(wv,nBwv_6000,wv,nBwv_300)
axis([0.1 100 0 1.1])
xlabel('Wavelength (micrometers)')
ylabel('Normalized Radiance (unitless)')
title('Plots of Normalized Radiance of two Blackbodies at 6000 K and 300 K (Part E)')
text(0.3, 0.2, '6000 K')
text(7, 0.2, '300 K')
