#!/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
#
# This version includes absorption by the ISM.
#
# 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;

# column density of the ISM, atoms per cm^3
$hy_column = 1.0E20;
# coefficients of fit to compute absorption cross-section
#    Taken from Morrison and McCammon,  ApJ 270, 119 (1983)
#    energy values are in keV
$cross_min_energy = 0.030;
$cross_max_energy = 10.0;
$cross_array_size = 14;
$energy_bot_array[0] = 0.030; $c0_array[0] =  17.3; $c1_array[0] =  608.1; $c2_array[0] = -2150;
$energy_bot_array[1] = 0.100; $c0_array[1] =  34.6; $c1_array[1] =  267.9; $c2_array[1] = -476.1;
$energy_bot_array[2] = 0.284; $c0_array[2] =  78.1; $c1_array[2] =   18.8; $c2_array[2] =    4.3;
$energy_bot_array[3] = 0.400; $c0_array[3] =  71.4; $c1_array[3] =   66.8; $c2_array[3] =  -51.4;
$energy_bot_array[4] = 0.532; $c0_array[4] =  95.5; $c1_array[4] =  145.8; $c2_array[4] =  -61.1;
$energy_bot_array[5] = 0.707; $c0_array[5] = 308.9; $c1_array[5] = -308.6; $c2_array[5] =  294.0;
$energy_bot_array[6] = 0.867; $c0_array[6] = 120.6; $c1_array[6] =  169.3; $c2_array[6] =  -47.7;
$energy_bot_array[7] = 1.303; $c0_array[7] = 141.3; $c1_array[7] =  146.8; $c2_array[7] =  -31.5;
$energy_bot_array[8] = 1.840; $c0_array[8] = 202.7; $c1_array[8] =  104.7; $c2_array[8] =  -17.0;
$energy_bot_array[9] = 2.471; $c0_array[9] = 342.7; $c1_array[9] =   18.7; $c2_array[9] =    0.0;
$energy_bot_array[10] = 3.210; $c0_array[10] = 352.2; $c1_array[10] =   18.7; $c2_array[10] =    0.0;
$energy_bot_array[11] = 4.038; $c0_array[11] = 433.9; $c1_array[11] =   -2.4; $c2_array[11] =    0.75;
$energy_bot_array[12] = 7.111; $c0_array[12] = 629.0; $c1_array[12] =   30.9; $c2_array[12] =    0.0;
$energy_bot_array[13] = 8.331; $c0_array[13] = 701.2; $c1_array[13] =   25.2; $c2_array[13] =    0.0;
$energy_bot_array[14] = 10.00; $c0_array[14] =   0.0; $c1_array[14] =    0.0; $c2_array[14] =    0.0;


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


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) {

    if (0 == 1) {
      $cross = calc_cross_section($energy);
      printf " keV %7.2f  cross %9.4e \n", $energy, $cross;
      next;
    }


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

    # now, compute the optical depth due to absorption by hydrogen
    #    along the line of sight
    $cross = calc_cross_section($energy);
    $opt_depth = $cross*$hy_column;
    $ext_factor = exp(0.0 - $opt_depth);
    $ext_photons = $num_photons*$ext_factor;

    printf " dens %9.4e  temp %8.1f nu %9.4e  keV %9.4e eff_nu %9.4e  photon %9.4e  ext_phot %9.4e \n",
           $n_cm3, $t, $nu, $energy, $eff_nu, $num_photons, $ext_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);
}


############################################
# Compute the cross-section of absorption by the ISM,
#   given an energy in keV.  
#   Taken from Morrison and McCammon, 
#      ApJ 270, 119 (1983)
#
sub calc_cross_section {
  my($energy_keV);
  my($index, $i, $cross);

  $energy_keV = $_[0];

  if ($energy_keV < $cross_min_energy) {
    return(0);
  }
  if ($energy_keV > $cross_max_energy) {
    return(0);
  }
  
  # figure out the appropriate index in arrays
  $index = -1;
  for ($i = 0; $i < $cross_array_size; $i++) {
    if (($energy_keV >= $energy_bot_array[$i]) &&
        ($energy_keV < $energy_bot_array[$i+1])) {
      $index = $i;
      last;
    }
  }
  if ($index == -1) {
    printf STDERR "calc_cross_section: invalid index for energy $energy_keV \n";
    exit(1);
  }
  
  # compute the cross-section per hydrogen atom
  $cross = $c0_array[$index] + 
           $c1_array[$index]*$energy_keV +
           $c2_array[$index]*$energy_keV*$energy_keV;
  if ($debug > 1) {
    printf " A: cross %12.5f \n", $cross;
  }
  $cross /= ($energy_keV*$energy_keV*$energy_keV);
  if ($debug > 1) {
    printf " B: cross %12.5f \n", $cross;
  }
  $cross *= 1E-24;

  return($cross);

}