% Matlab script to...
% Calculate and plot Planck's law for given temperature (T).
%
% Originator: G. J. Hakim, University of Washington
%
% To run: open a matlab session and type: planck
%
%
disp('====================================================')
disp(' ')
disp('Compute Planck function for monochromatic irradiance.')
disp('Works well for temperatures less than ~10^5 K')
disp(' ')
disp('====================================================')
disp(' ')

% constants:
c1 = 3.7400e-16; c2 = 0.0144;

% read temperature
T = 300;
while (T >= 0)
  T = input('Enter the absolute temperature (-1 to exit)')

  % Use Wien's law to center curve
  lmax = 2.88e-3/T;
  lambda = lmax/10:1e-8:lmax*10;

  E = c1.*lambda.^(-5)./( exp((c2/T)./lambda) -1 );

  if (T > 0) 
	 figure(1); clf
	 plot(lambda,E,'r','linewidth',2);
	 xlabel ('Wavelength (m)','fontsize',14,'fontweight','bold')
	 ylabel ('Monochromatic Irradiance (W/m^3)','fontsize',14,'fontweight','bold')
	 h = gca; 
	 set(h,'fontweight','bold')
	 set(h,'linewidth',2)
	 cT = num2str(T);
	 ttl = strcat('Planck function for T = ',cT,' K');
	 title (ttl,'fontsize',16,'fontweight','bold');
  end
end

disp(' Goodbye!')


