#!/usr/bin/perl
#
# Calculate several quantities as a function of temperature
#
#    a. the ratio of hydrogen in the excited n=2 state to hydrogen
#             in the ground n=1 state
#
#    b. the ratio of ionized to neutral hydrogen
#
#    c. the ratio of hydrogen in the n=2 state to all hydrogen
#
# MWR 1/31/2004
#
# Modified to use number density of electrons ne which is appropriate
#   for photosphere of expanding type II SN.
#
# MWR 1/31/2004


$debug = 0;

# perform all calculations in MKS
$c = 2.9979e8;
$h = 6.626e-34;
$k = 1.38e-23;
$mh = 1.67e-27;
$me = 9.1094e-31;
$pi = 3.14159;
$qe = 1.6022e-19;

# partition functions of hydrogen in ground and ionized stages
$z1 = 1.0;
$z2 = 1.0;

# energy difference between hydrogen in first and second energy levels
#   first, in eV
$e2_e1 = 10.2;
#   now, in Joules
$e2_e1 *= $qe;

# ionization potential of hydrogen from the ground state
#   first, in eV
$chi = 13.6;
#   now, in Joules
$chi *= $qe;

# number density of free electrons 
# we start with mass density of grams per cubic cm
$rho = 3.0e-12;
# convert to MKS (kg per cubic meter)
$rho *= (0.001)/(1e-6);
# now convert to atoms per cubic meter, assuming all hydrogen
$rho /= $mh;
# now assume all atoms are ionized, and calc number density
#   of electrons per cubic meter
$ne = $rho;
printf "ne is %10.4e \n", $ne;


$start_temp = 3000;
$end_temp = 300000;
$dtemp = 500;

for ($temp = $start_temp; $temp <= $end_temp; $temp += $dtemp) {

  # calculate ratio of ionized to neutral hydrogen
  $x1 = 2.0*($z2/$z1)*(1.0/$ne);
  $x2 = (2.0*$pi*$me*$k*$temp)/($h*$h);
  $x2 = $x2**1.5;
  $x3 = ($chi)/($k*$temp);
  $x3 = exp(-$x3);
  $ion_to_neut = $x1*$x2*$x3;
  $ion_to_tot = ($ion_to_neut)/($ion_to_neut + 1.0);

  if ($debug > 0) {
    printf "Q %7.0f   %10.4e %10.4e %10.4e \n", $temp, $x1, $x2, $x3;
  }

  printf "%7.0f   %10.4e %10.4e ", $temp, $ion_to_neut, $ion_to_tot;


  # calculate ratio of hydrogen in n=2 to n=1
  $x1 = (8.0/2.0);
  $x2 = ($e2_e1)/($k*$temp);
  $x2 = exp(-$x2);
  $n2_to_n1 = $x1*$x2;

  printf "%7.0f   %10.4e ", $temp, $n2_to_n1;


  # now, calculate the ratio of neutral hydrogen in n=2 to all hydrogen
  $n2_to_all = (1.0 - $ion_to_tot)*($n2_to_n1);
  printf "   %10.4e ", $n2_to_all;

  printf "\n";
}


exit 0;
