#!/usr/bin/perl

# Ay 126 problem set 1.
# Ionization fraction in the Stromgen sphere
# refer to the problem set at http://www.astro.caltech.edu/~varun/ay126/hw1.pdf
# Varun Bhalerao

# after going through the algebra, we get:
# (1-y)/y^2 = r^2/p * 0.00062739892782623223
# here r = R/R_s, p = 1-r^2
# lets solve this and generate a plot

open (OUT, ">ionization_fraction");

for ($r=0.001; $r<1.0; $r+=0.001){
    $p = 1 - $r**2;
    $k = $r**2 / $p * 0.000627;
    $y = (-1.0 + sqrt(1 + 4*$k)) / 2.0 / $k;
    print OUT "$r $y\n";
}
print OUT "1.0 0.0\n";      #by definition !

# plot the file "ionization_fraction" in gnuplot !