#!/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.
#
# We compute the spectrum of a large sphere of hot gas,
#   assuming some form of density profile and
#   temperature profile, and performing the
#   integration for projection onto the sky.
#
# MWR 4/11/2011

use POSIX;

$debug = 0;

if (0 == 1) {
  for ($i = 0; $i < 10; $i += 0.1) {
    $rho = compute_density(1.0, 1.0, $i);
    printf " %6.1f  %9.4e \n", $i, $rho;
  }
  exit 0;
}

# constants
$pi = 3.14159;
# charge on ion = 1 for Hydrogen
$Z = 1;
# meters per kpc
$kpc_to_meter = 3.08E19;
# central density, atoms/m^3
$n = 1e4;
$n_cm3 = $n*1e6;
# density follows isothermal sphere with this core radius
#   (units are kpc)
$core_radius_kpc = 1.0;
$core_radius_meter = $core_radius_kpc*$kpc_to_meter;
# Boltzmann constant
$k = 1.38E-23;
# speed of light
$c = 2.9979E8;
# Planck's constant
$h = 6.626E-34;
# central temperature, Kelvin
$T0 = 100E6;
# energy of 1 keV, in Joules
$one_kev = 1.602E-16;

# column density of the ISM, atoms per cm^3
$hy_column = 0.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;


# how far off-center are we looking?  units are core radii
$mult_x_kpc = 1.211528;
for ($x_kpc = 0.1; $x_kpc < $core_radius_kpc*20.0; $x_kpc *= $mult_x_kpc) {
  $x_meter = $x_kpc*$kpc_to_meter;
  
  # we integrate over this range in the "z" coordinate
  #    units are core radii
  #    go from 0 to max, then multiply by 2 
  $z_start = 0;
  $z_end_kpc = $core_radius_kpc*20.0;
  $d_z_kpc = $core_radius_kpc*0.1;
  $z_end_meter = $z_end_kpc*$kpc_to_meter;
  $d_z_meter = $d_z_kpc*$kpc_to_meter;
  # we need these for computing volume properly
  $d_x_meter = $d_z_meter;
  $d_y_meter = $d_z_meter;
  $d_volume_meter = $d_x_meter*$d_y_meter*$d_z_meter;
  
  
  # make all entries are zero in arrays
  #    energy_array[]     energy of photons in energy bin
  #    photon_array[]     how many photons in each energy bin?
  #    close_photon_array[]     photons per energy bin at
  #                                location closest to center
  #                                for chord (i.e. z = 0)
  #
  for ($i = 0; $i < 1e6; $i++) {
    $energy_array[$i] = 0.0;
    $photon_array[$i] = 0.0;
    $close_photon_array[$i] = 0.0;
  }
  
  
  
  for ($z_meter = $z_start; $z_meter <= $z_end_meter; $z_meter += $d_z_meter) {
  
    $radius_meter = sqrt($z_meter*$z_meter + $x_meter*$x_meter);
    $radius_kpc = $radius_meter/$kpc_to_meter;
  
    # compute density at this location
    $rho = compute_density($n, $core_radius_meter, $radius_meter);
    $rho_cm3 = $rho*1e6;
  
    # compute temperature at this location
    $t = compute_temperature($T0, $core_radius_meter, $radius_meter);
  
  
    if ($debug > 0) {
      printf " z %9.4e  radius_meter %9.4e  rho %9.4e  T  %9.4e\n", 
             $z_meter, $radius_meter, $rho, $t;
    }
  
  
  
    $energy_index = 0;
    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*($rho_cm3*$rho_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*($rho_cm3*$rho_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;
  
      # add these photons to the photon_array[] for this wavelength
      $energy_array[$energy_index] = $avg_energy_per_photon;
      $photon_array[$energy_index] += $ext_photons;
  
      if ($debug > 0) {
        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;
      }
  
      # save the spectrum of the point closest to center
      #   for comparison with integrated spectrum
      if ($z_meter == $z_start) {
        $close_photon_array[$energy_index] = $photon_array[$energy_index];
      }
  
      
      $energy_index++;
  
    }
  
  }


  # now we can print out the integrated spectrum at this X position
  #   and the total number of photons observed at this position
  $sum_photons = 0;
  for ($i = 0; $i < $energy_index; $i++) {
  
    $e_Joules = $energy_array[$i];
    $nu = $e_Joules/$h;
    $avg_photon_kev = $energy_array[$i]/$one_kev;

    printf "  X %7.2f temp %9.2e  nu %9.4e keV %9.4e  photons %9.4e  close %9.4e \n",
       $x_kpc, $t, $nu, $avg_photon_kev, $photon_array[$i], $close_photon_array[$i];

    $sum_photons += $photon_array[$i];
  }
  printf "#  X %7.2f  temp %9.2e  total_photons %9.4e \n",
         $x_kpc, $t, $sum_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);

}


#############################################
# Compute the density at the given radius R
#   and core radius R_c, given a central density
#   rho0.  We assume a King model
#   of isothermal sphere
#
# Arguments:  rho0, R_c, R
# 
sub compute_density {
  my($R, $R_c, $rho0);
  my($beta);

  $beta = 1.0;

  $rho0 = $_[0];
  $R_c = $_[1];
  $R = $_[2];

  $rho = $rho0 * ((1.0 + ($R/$R_c)**2.0)**(-3.0*$beta/2.0));
  
  return($rho);
}
  

#############################################
# Compute the temperature at the given radius R
#   and core radius R_c, given a central temperature
#   T0.  
#
# Arguments:  T0, R_c, R
# 
sub compute_temperature {
  my($R, $R_c, $T0);

  $T0 = $_[0];
  $R_c = $_[1];
  $R = $_[2];

  if (0 == 1) {
    $temp = $T0 * ((1.0 + ($R/$R_c))**(-1));
  }
  else {
    $temp = $T0;
  }
  
  return($temp);
}