% problem 7-2 in Bohm-Vitense

% find He+/He, He++/He+, and Pe for the following conditions

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

% constants for Helium
xion1 = 24.58 ;
xion2 = 54.4;
upp = 1;
up = 2;
u = 1;


% Boltzman constant
k = 1.38e-16;

% find theta
theta = 5040 / temp;


% calculate the constant coefficents in the Saha equation
c1 = log10(up/u) + log10(2) + 2.5 * log10(temp) - xion1 * theta - 0.48;
c2 = log10(upp/up) + log10(2) + 2.5 * log10(temp) - xion2 * 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 He+/He and He++/He+
lratio1 = c1 - log10(pe);
lratio2 = c2 - log10(pe);

% find the ratio
ratio1 = 10^lratio1;
ratio2 = 10^lratio2;

% calculate the abundance of H
Hep = n / (3*ratio2 + 2 + 1/ ratio1);
Hepp = ratio2 * Hep;
He = n - 3 * Hepp - 2* Hep;

% recalculate a new ne
ne = 2*Hepp + Hep;

end

% print out the final values of interest

ratio1
lratio1

ratio2
lratio2

  
ne

Pe = ne * k * temp

Hepp

Hep

He