#!/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
#
# Add a piece to integrate flux over some range of wavelengths.
#
# MWR 2/27/2015
# 
# This version computes the number of photons emitted short of the
#   13.6-eV, 912-Angstrom Lyman limit.
#
# MWR 10/27/2025

use POSIX;
$debug = 0;

$MYPI = 3.1415926;
$TINY = 1.0e-12;

$temp = $ARGV[0];

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

$filter = "none";
if ($filter eq "J") {
  $central_lambda = 1.26E-6;
  $width_lambda = 0.16E-6;
}
if ($filter eq "H") {
  $central_lambda = 1.60E-6;
  $width_lambda = 0.23E-6;
}
if ($filter eq "K") {
  $central_lambda = 2.15E-6;
  $width_lambda = 0.26E-6;
}
if ($filter eq "L") {
  $central_lambda = 3.80E-6;
  $width_lambda = 0.35E-6;
}
if ($filter eq "M") {
  $central_lambda = 4.70E-6;
  $width_lambda = 0.30E-6;
}
if ($filter ne "none") {
  $start_lambda = $central_lambda + 0.5*$width_lambda;
  $end_lambda = $central_lambda - 0.5*$width_lambda;
}
else {
  $start_lambda = 1e10;
  $end_lambda = 1e-18;
}
# $dlog_lambda = -0.0001;
$dlog_lambda = -0.001;



$sum_flux = 0.0;

$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);
# wavelength of Lyman limit photon (Angstroms)
$lambda_lyman_A = 912.0;
$lambda_lyman_m = ($lambda_lyman_A)*1.0e-10;
if ($debug > 0) {
  printf " lambda_lyman_A %12.5e  lambda_lyman_m %12.5e \n",
     $lambda_lyman_A, $lambda_lyman_m;
}

# number of photons with energies above Lyman limit
$sum_lyman_photons = 0;
# number of all photons emitted at all energies
$sum_all_photons = 0;


$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 < 50) {
    $bot = exp($exponent) - 1;
    # printf " exp %10.3f bot is %12.5e \n", $exponent, $bot;
    if ($bot < $TINY) {
      $blam = 1e-90;
    }
    else {
      $blam = $top/$bot;
    }
  }
  else {
    $blam = 1e-90;
  }
  # printf "exponent %9.4e  bot %9.4e \n", $exponent, $bot;


  $logblam = log10($blam);

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

  # add up the energy in this little piece
  $next_lambda = pow(10.0, ($logl + $dlog_lambda));
  $d_lambda = $lambda - $next_lambda;
  $energy = $MYPI*$blam*$d_lambda;
  # add up the number of photons in this piece,
  #   and the number with energies above Lyman limit
  $energy_per_photon = $h*$nu;
  $num_photons = $energy / $energy_per_photon;
  $sum_all_photons += $num_photons;
  if ($lambda < $lambda_lyman_m) {
    $sum_lyman_photons += $num_photons;
  }
  if ($debug > 0) {
    printf " energy_per_photon %9.4e  num_photon %9.4e sum_lyman_photon %9.4e \n",
           $energy_per_photon, $num_photons, $sum_lyman_photons;
  }
  $sum_energy += $energy;

  printf "%9.4e %9.4e %9.4e %9.4e %9.4e %9.4e  %10.5e  %9.4e %9.4e \n", 
       $logl, $lambda, $blam, $logblam, 
                  $nu, $lognu, $sum_energy,
                  $sum_all_photons, $sum_lyman_photons;
}

exit 0;