#!/usr/bin/perl
#
# Compute the spectrum for thermal brehmsstrahlung
#   over a range of temperatures and densities
#   and print results in a nice table.
#
# MWR 4/11/2011

use POSIX;

$debug = 0;

# constants
$pi = 3.14159;
# charge on ion = 1 for Hydrogen
$Z = 1;
# average density, atoms/m^3
$n = 1e4;
$n_cm3 = $n*1e6;
# Boltzmann constant
$k = 1.38E-23;
# speed of light
$c = 2.9979E8;
# Planck's constant
$h = 6.626E-34;
# temperature, Kelvin
$start_t = 0.1E6;
$end_t = 50E6;
$mult_t = 3.162278;
# energy of 1 eV, in Joules
$one_ev = 1.602E-19;

# energy, eV
$start_energy = 0.01;
$end_energy = 150;
$mult_energy = 1.467799;

for ($t = $start_t; $t <= $end_t; $t *= $mult_t) {

  if ($debug > -1) { 
    printf "#  temp %9.4e  \n", $t;
  }

  for ($energy = $start_energy; $energy <= $end_energy; $energy *= $mult_energy) {

    # compute frequency in Hz
    $e_Joules = $energy*$one_ev;
    $nu = $e_Joules/$h;
    $hnukt = ($h*$nu)/($k*$t);
    if ($debug > 0) { 
      printf " energy %7.1f eV  %9.4e J  nu %9.4e  hv/kT %9.4e \n",
          $energy, $e_Joules, $nu, $hnukt;
    }

    $gff = gaunt_ff($t, $nu);
    if ($debug > 0) {
      printf " gff %9.4e  \n", $gff;
    }

    # compute the emissivity (erg/s-cm^3-Hz)
    $eff_nu = 6.8E-38*($n_cm3*$n_cm3)*sqrt($t)
                 *exp(-$hnukt)*$gff;
    printf " dens %9.4e  temp %8.1f nu %9.4e  eV %9.4e eff_nu %9.4e \n",
           $n_cm3, $t, $nu, $energy, $eff_nu;
    

  }

  printf "\n\n\n";

}





exit 0;


#############################################
# Compute Gaunt factor, given temperature (K)
#   and frequency (Hz)
#
# 
sub gaunt_ff {
  my($temp, $nu);
  my($gff);

  $temp = $_[0];
  $nu = $_[1];

  $gff = (3.0/sqrt($pi))*log((9*$k*$temp)/(4*$h*$nu));

  # if the Gaunt factor is less than 0,  just return 0
  if ($gff < 0) {
    $gff = 0;
  }
 
  return($gff);
}