#!/usr/bin/perl
#
# Calculate a blackbody spectrum and print out intensity as
#   as a function of wavelength.  Include log(lambda) and log(Intensity),
#   too.
#
# usage:  bb.pl  temp
#
# where "temp" is temperature in Kelvin
# 
# MWR 7/23/2003
#
# fixed error in computation of the "bot" piece -- was missing "-1".
# MWR 12/9/2003

use POSIX;

$temp = $ARGV[0];

$start_lambda = 1e10;
$end_lambda = 1e-18;
$dlog_lambda = -0.05;

$h = 6.626e-34;
$c = 3.0e8;
$k = 1.38e-23;
# this should be 2*h*c^2
$const = 2.0*$h*$c*$c;
$exponent_const = ($h*$c)/($k*$temp);


$log_start = log10($start_lambda);
$log_end = log10($end_lambda);
for ($logl = $log_start; $logl > $log_end; $logl += $dlog_lambda) {

  $lambda = pow(10.0, $logl);

  $top = $const/(pow($lambda, 5.0));
  $exponent = $exponent_const/$lambda;
  if ($exponent < 100) {
    $bot = exp($exponent) - 1;
    if ($bot != 0.0) {
      $blam = $top/$bot;
    } else {
      $blam = 1e-90;
    }
  }
  else {
    $blam = 1e-90;
  }
  # printf "exponent %9.4e  bot %9.4e \n", $exponent, $bot;


  $logblam = log10($blam);

  $nu = $h*$c/$lambda;
  $lognu = log10($nu);

  printf "%9.4e %9.4e %9.4e %9.4e %9.4e\n", $logl, $lambda, $blam, $logblam, $lognu; 
}

exit 0;

