#!/usr/bin/perl
#
# Compute the spectrum for thermal brehmsstrahlung
#   over a range of temperatures and densities
#   and print results in a nice table.
#
# This version integrates emissivity over frequency
#   and computes the number of photons one would
#   receive within each energy interval
#
# 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 = 50E8;
$mult_t = 3.162278;
# energy of 1 keV, in Joules
$one_kev = 1.602E-16;

# energy, keV
$start_energy = 0.10;
$end_energy = 50.0;
$mult_energy = 1.211528;

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_kev;
    $nu = $e_Joules/$h;
    $hnukt = ($h*$nu)/($k*$t);
    if ($debug > 0) { 
      printf " energy %7.1f keV  %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;
    if ($debug > 0) {
      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;
    }

    # now, repeat for the next energy 
    $next_energy = $energy * $mult_energy;
    $next_e_Joules = $next_energy*$one_kev;
    $next_nu = $next_e_Joules/$h;
    $next_hnukt = ($h*$next_nu)/($k*$t);
    if ($debug > 0) { 
      printf " next: energy %7.1f keV  %9.4e J  nu %9.4e  hv/kT %9.4e \n",
          $next_energy, $next_e_Joules, $next_nu, $next_hnukt;
    }
    $next_gff = gaunt_ff($t, $next_nu);
    if ($debug > 0) {
      printf " next_gff %9.4e  \n", $next_gff;
    }
    # compute the emissivity (erg/s-cm^3-Hz)
    $next_eff_nu = 6.8E-38*($n_cm3*$n_cm3)*sqrt($t)
                 *exp(-$next_hnukt)*$next_gff;
    if ($debug > 0) {
      printf " next_eff_nu %9.4e \n", $next_eff_nu;
    }

    # compute total energy emitted in this energy interval
    $energy_in_bin = (0.5*($eff_nu + $next_eff_nu))*($next_nu - $nu);
    # now convert from energy to number of photons,
    #     using the average energy per photon within interval
    $avg_energy_per_photon = 0.5*$h*($nu + $next_nu);
    $num_photons = $energy_in_bin / $avg_energy_per_photon;

    if ($eff_nu > 0) {
      $ratio = $num_photons / $eff_nu;
    }
    else {
      $ratio = 0;
    }
    if ($debug > 0) {
      printf " ratio %9.4e \n", $ratio;
    }

    printf " dens %9.4e  temp %8.1f nu %9.4e  keV %9.4e eff_nu %9.4e  photon %9.4e \n",
           $n_cm3, $t, $nu, $energy, $eff_nu, $num_photons;
    

  }

  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);
}