#!/usr/bin/perl
#
# Calculate a blackbody spectrum and print out intensity 
#   as a function of wavelength over a range of wavelengths.
#   Print numbers to stdout with a format like this:
#
#        nu    lambda     spec 
#
# where 
#           "nu"          is the frequency in Hz
#           "lambda"      is the wavelength in meters
#           "spec"        is the energy per square meter per second
#                              per steradian per meter
#
# usage:  bb.pl  temp  
#
# where "temp" is temperature in Kelvin
# 
# MWR 8/19/2003
# 
# Modified to compute spectrum in "fnu" units as well.
#
# MWR 3/8/2011


use POSIX;

$temp = $ARGV[0];

# wavelengths in meters
$start_lambda = 0.5e-8;
$end_lambda = 2.5e-6;
$mult_lambda = 1.05;
$d_lambda = 0.1e-6;

$pi = 3.14159;

$h = 6.626e-34;
$c = 3.0e8;
$k = 1.38e-23;
$const = 2.0*$h*($c*$c);
$exponent_const = ($h*$c)/($k*$temp);

printf "#%8s %9s %9s %9s\n", "freq(Hz)", "lambda(m)", "Blam", "Bnu";
 

for ($lambda = $start_lambda; $lambda <= $end_lambda; $lambda *= $mult_lambda) {

  $logl = log10($lambda);

  $top = $const/(pow($lambda, 5.0));
  $exponent = $exponent_const/$lambda;
  # we are careful here to avoid numerical underflow
  if ($exponent < 100) {
    $bot = exp($exponent) - 1;
    $blam = $top/$bot;
  }
  else {
    $blam = 1e-90;
  }

  # compute frequency, too
  $nu = $c/$lambda;

  # compute Bnu
  $top = 2*$h*($nu*$nu*$nu) / ($c*$c);
  $nu_exponent = ($h*$nu) / ($k*$temp);
  if ($nu_exponent < 100) {
    $bot = exp($nu_exponent);
    $bnu = $top / $bot;
  }
  else {
    $bnu = 1e-90;
  }

  printf "%9.4e %9.4e %9.4e %9.4e \n", $nu, $lambda, $blam, $bnu;

}

exit 0;

