% absorption program for Bohm-Vitense problem 7-3

% find the hydrogen bound absorption coefficent
% for levels n <=4

temp = 5800;
theta = 5040/ temp;
% temp = 5040 / theta;

% ionization energy
xion = 13.6;
ergperev = 1.6022e-12;
xion = xion * ergperev;

% physical constants
c = 3e10;
z = 1;
g = 1;
h = 6.57e-27;
k = 1.38e-16;
sigma = 5.67e-5;
m = 9.1094e-28;
e = 4.8032e-10;
pi = 3.14159265;

% calculate the coefficient for the bound free coefficent
mult =  (64.0 * pi**4 / (3 * sqrt(3.0))) * z**4 * m * e**10/ (c * h**6)  * g;

% loop through an index
for j = 1:130

% calculate wavelength in angstroms and cm
i = 200 + j *100;
lambda = i * 1e-8;

% find the frequency and energy
energy = h * c /lambda;
freq = c / lambda;

% determine the number of states in the system
x = (energy/ xion) **(-0.5);
n = floor(x) + 1;

% sum over all excited states from the lowest to the eighth
atmp = 0;
for nn = n:4

% calculate the boltmann occupation numbers
gn = 2.0 * nn**2.0;
g1 = 2.0;

% calculate the excitation energy
xn = xion * (1.0 - 1.0 / nn**2)  ;

% find the number in the state using the boltzman equation
numbn = gn/ g1 * exp(-xn/(k*temp));

% find the absorption coefficent for this level
an = mult / (nn**5 * freq**3);

% sum to find the total absorption coefficent
atmp = atmp + an * numbn;
end

% save the results into an array
ll(j) = log10(i)     ;
kappa(j) = log10(atmp);
ni(j) = n;

end


gset output 'bh7p3.ps'
gset term post color
% plot the results - log10(kappa/atom) vs log10(wavelength in angstroms)
plot(ll,kappa)