#!/usr/bin/perl
#
# Calculate a blackbody spectrum and print out intensity in fnu
#   as a function of frequency.  Include log(lambda) and log(Intensity),
#   too.
#
# usage:  bb_fnu.pl  temp
#
# where "temp" is temperature in Kelvin
# 
# MWR 12/9/2003
#

use POSIX;

$debug = 0;

$temp = $ARGV[0];

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

# all units here are MKS
$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/($k*$temp);

# calculate the starting and ending boundaries in frequency space
$start_freq = $c/$start_lambda;
$end_freq = $c/$end_lambda;
$dlog_freq = 0.1;

$energy_sum = 0.0;

$log_start = log10($start_freq);
$log_end = log10($end_freq);
for ($logf = $log_start; $logf < $log_end; $logf += $dlog_freq) {

  $freq = pow(10.0, $logf);

  $top = $const*($freq*$freq*$freq);
  $exponent = $exponent_const*$freq;
  if ($exponent < 100) {
    $bot = exp($exponent) - 1;
    $bnu = $top/$bot;
  }
  else {
    $bnu = 1e-90;
  }
  # printf "exponent %9.4e  bot %9.4e \n", $exponent, $bot;


  # compute the energy contributed by a bin of frequencies at this 
  #   frequency
  $next_freq = pow(10.0, $logf + $dlog_freq);
  $d_freq = $next_freq - $freq;
  $energy = $bnu*$d_freq;
  $energy_sum += $energy;
  if ($debug > 0) {
    printf "  freq %9.4e  d_freq %9.4e  bnu %9.4e  energy %9.4e  sum %9.4e \n",
            $freq, $d_freq, $bnu, $energy, $energy_sum;
  }


  $logbnu = log10($bnu);

  $lam = $c/$freq;
  $loglam = log10($lam);

  printf "%9.4e %9.4e %9.4e %9.4e %9.4e %9.4e\n", $logf, $freq, $bnu, $logbnu, $lam, $loglam; 
}

exit 0;