% problem 7-1 in Bohm-Vitense

% find H+/H and Pe for the following conditions

% system constants
temp = 10080;
pg = 1e3;

% constants for Hydrogen
xion = 13.6 ;
up = 1;
u = 2;

% Boltzman constant
k = 1.38e-16;

% find theta
theta = 5040 / temp;


% calculate the constant coefficents in the Saha equation
c0 = log10(up/u) + log10(2) + 2.5 * log10(temp) - xion * theta - 0.48;

% find the number of particles per cm^3
n = pg / (k* temp);

% guess an initial ne
ne = n / 2;

% loop over the system for convergence
for i = 1:20

% find a new pe
pe = ne * k * temp;

% calculate the log ratio of H+/H
lratio = c0 - log10(pe);

% find the ratio
ratio = 10^lratio;

% calculate the abundance of H
H = n / (2* ratio + 1);

% recalculate a new ne
ne = (n - H)/2;

end

% print out the final values of interest

ratio
lratio

ne