
(* :Title: BlackBodyRadiation *)

(* :Author: Barbara Ercolano *)

(* :Summary:
This package provides functions giving the basic properties of black body 
radiation at a specified temperature, 
including the total power emitted, 
the peak wavelength and and the power emitted by a narrow band about the 
peak. The package also provides a function for plotting a blackbody 
profile at a specified temperature.
*)

(* :Context: Miscellaneous`BlackBodyRadiation` *)

(* :Copyright: Copyright 1997-1999, Wolfram Research, Inc. *)

(* :Source: *)

(* :Mathematica Version: 3.0 *)

(* :History:
   Version 1.0 by Barbara Ercolano, 1997
   Version 1.1 changes by ECM, 1998
*)

(* :Package Version: 1.1 *)

BeginPackage["Miscellaneous`BlackBodyRadiation`", "Miscellaneous`Units`",
		"Miscellaneous`SIUnits`", "Miscellaneous`PhysicalConstants`"];

PeakWavelength::usage = 
  "PeakWavelength[temperature] gives the wavelength of the maximum emission \
of a blackbody at the specified temperature."

BlackBodyProfile::usage=   
  "BlackBodyProfile[temperature1, temperature2, ...] plots blackbody \
spectral distribution profiles at the specified temperatures.";  

TotalPower::usage=   
"TotalPower[temperature] gives the total radiative power emitted by a
black body at the specified temperature.";

MaxPower::usage=
"MaxPower[temperature, wavelengthband] gives the radiative power emitted by
a blackbody in the specified wavelength band about the peak
wavelength at the specified temperature.";

Begin["`Private`"]


MaxPower::wvl = "Value `1` is not an appropriate wavelength band.";

PeakWavelength::temp =
BlackBodyProfile::temp =
TotalPower::temp =
MaxPower::temp =
  "Value `1` is not an appropriate temperature.";

BlackBodyProfile::temps =
	 "Values `1` are not appropriate temperatures.";

convertToKelvin[T_] := Module[{a},
		Replace[T, {a_.*( Celsius|Centigrade) :> N[a] + 273.15,
		 a_. * Fahrenheit :>  (N[a]-32)*5/9 + 273.15,
		a_. * Rankine :> N[a] *5/9, a_. * Kelvin :> N[a]}]]

planckConstant = PlanckConstant/(Joule Second);
speedOfLight = SpeedOfLight/(Meter/Second);
boltzmannConstant =BoltzmannConstant/(Joule/Kelvin);

pf1 = (2 Pi planckConstant  speedOfLight^2);
pf2 = planckConstant speedOfLight/boltzmannConstant;

planckFunction[t_, x_] := (pf1/ x^5)/(Exp[pf2/(x*t)]-1)

peakWavelength[t_] := (2.90 10^-3)/t;

PeakWavelength[T_] := Module[{t},peakWavelength[t] Meter
	  /; (t=convertToKelvin[T];
				 If[MatchQ[t, _Real | _Integer], True,
			   Message[PeakWavelength::temp, T]; False, False])]

maxline[t_]:=  Module[{p=peakWavelength[t],q},
		q = planckFunction[t,p];
		Graphics[{Dashing[{0.01, 0.01}], Line[{{p, 0},  {p, q}}],
			Text[t K, {p, q}, {0, -1}]}]]

blackBodyProfile[t_, opt___?OptionQ]:=
   Module[{x},
	{maxline[t], Plot[planckFunction[t,x], {x, 0,  6 peakWavelength[t]},
		opt, DisplayFunction -> Identity]}
   ];

BlackBodyProfile[Ts__, opt___?OptionQ] := Module[{ts},
    Show[Map[blackBodyProfile, ts],
        opt,
        DisplayFunction->$DisplayFunction,
        PlotRange->All,
        Axes -> True,
        AxesLabel -> {"Wavelength, m",
                      "Power, W/m^2"}] /;
           (ts=Map[convertToKelvin, {Ts}];
	    If[MatchQ[ts, {(_Real | _Integer)...}], True,
	       If[Length[{Ts}]>1,
			   Message[BlackBodyProfile::temps, {Ts}],
			   Message[BlackBodyProfile::temp, Ts]]; False, False])]

(* original code, see below.
totalPower[t_]:=
  Module[{x, slwcon, ovfl, ncvb,  nint},
	slwcon = (Head[NIntegrate::slwcon] === $Off);  Off[NIntegrate::slwcon];
 	ovfl = (Head[General::ovfl] === $Off); Off[General::ovfl];
	ncvb = (Head[NIntegrate::ncvb] === $Off); Off[NIntegrate::ncvb];
	nint = NIntegrate[planckFunction[t, x], {x, 0, Infinity}] Watt/Meter^2;
	If[!slwcon, On[NIntegrate::slwcon]];
	If[!ovfl, On[General::ovfl]];
	If[!ncvb, On[NIntegrate::ncvb]];
	nint
  ]
*)
(* The total power is given directly by the Stephan-Boltzmann law;
   no need to compute the integral numerically. We can verify the
   result by symbolic integration. Note that because we want to
   express the units in Watt/Meter^2, we strip the units from the
   computation of \[Sigma] by the use of First[Cases[..., _?NumberQ]],
   which is marginally cleaner than multiplying by the inverse of the
   units involved in this instance. The units have already been
   stripped from t elsewhere in the code. *)
totalPower[t_] :=
    First[Cases[
      2 Pi^5 BoltzmannConstant^4/(15 SpeedOfLight^2 PlanckConstant^3),
     _?NumberQ]] * t^4 (Watt/Meter^2)

TotalPower[T_]:=
 Module[{t},
   (
   totalPower[t]
   ) /; (t=convertToKelvin[T];
	 If[MatchQ[t, _Real | _Integer], True,
	    Message[TotalPower::temp, T]; False, False])
 ]

MaxPower[T_, delamda_] :=
 Module[{t,delam},
  (
  0.657548 (delam/peakWavelength[t]) totalPower[t]
  ) /;
   (t=convertToKelvin[T];
    delam = N[Convert[delamda, Meter]/Meter];
    If[MatchQ[t, _Real | _Integer],
        If[MatchQ[delam, _Real | _Integer], True,
            Message[MaxPower::wvl, delamda]; False,
             False],
       Message[MaxPower::temp, T]; False,
       False])
 ]

End[];

EndPackage[];
