HOMEWORK 3 HELPFUL NOTES

Please be more liberal (thorough and complete) for the explanations in
your write ups.  Please attach (staple) your codes to the END of your
homework, and do not embed them within the solutions.

I have included tables of important quantities in the files

partfuncs.tab - partition functions (log10)
ionpots.tab   - ionization potentials (eV)
amu.tab       - atomic mass units of elements

The files all contain a header line with useful information.

NOTE regarding the indexing the elements in these files.  "k" is
employed to denote the species and "j" is employed to denote the
ionization stage.  Per the class notes, "j" is indexed from 1 for the
neutral.  What is different in these files is that not all species
from the periodic table are included, so the "k" index does not
correspond 1:1 with the atomic number, as it would if all elements
were tabulated.  I just wanted to bring that point out in case it
seemed confusing.  Note that the "k" index for a given species is
consistent between the tables.


PROBLEM 1.

I did not specify the number of temperature steps or increment, so you
can define these to your liking.  I would recommend steps of 100 K.
You will need the partition functions at various temperatures, thus
you will need to interpolate those data.  Note that the table is given
as a function of theta=5040/T, increasing to the right (which means
temperature is decreasing to the right).  I suggest you interpolate in
theta, even though you make your calculations in terms of T.


PROBLEM 2.

I have provded the solutions to this problem in the file
"Prob2-solutions.txt"; please format you solutions identical to those
presented in this file.

You are free to solve Eq. 1 using any method, however, I highly
recommend Zbrent's method to root solve for the electron density.  It
is extremely stable and fast.  In this directory, I have included the
Fortran functions (``zbrent.f'') and C code (``zbrent.c'').  It is
highly likely that IDL and python have the zbrent function.
Explanations for the use of the zbrent function can be found in the
book "Numerical Recipes".  Basically, you provide it the function you
wish to root solve (Eq. 1), the lower and upper bound over which to
search for the root, and the tolerance at which to find it.

You will require the partition functions at T=6000, which will require
interpolation.  You will require the ionization potentials for each
ionization stage for all species. You will also require the masses of
the elements in atomic mass units (for mass density calculations).
The data for these are given in the *.tab files.

When solving for the mean molecular weights and mass densities, I do
NOT want you to simply solve the pressure-density-temperature
relation, i.e, P=(rho/mu*m_atm)*kT.  That is basically a way to deduce
or cross-check a value given all other quantities under the assumption
of pressure equilibrium (which we assume holds).  I DO want you to
compute all quantities from the principles derived in class, including
particle number conservation, particle mass and mass fraction
relations, and charge conservation.


PROBLEM 3.

Once you have successfully solved Problem 2, this problem basically
entails putting a temperature loop in your code.  Again, I did not
specify the number of temperature steps or increment, so you can
define these to your liking.  I would recommend steps of 100 K


NOTE ABOUT INTERPOLATING PARTITION FUNCTIONS.  

Note that some partition functions are constant with temperature.
Cleary, you do not need to interpolate these.

If you use high order interpolation, such as the cubic spline methods
of SPLINT and SPLINE, you can get some "ringing".  The tabulated
partition functions (log10) vary slowly with T over certain T ranges
and then change more rapidly over other T ranges.  In the range of T
where the transition occurs, higher order interpolation can result in
"ringing" about the data points, and this can result in negative
values, which are unphysical since the tabulated data are in log10.  I
recommend that you test your interpolation as a function T carefully
before you proceed with your calculations.  You may wish to use a low
order spline or a lower order polynomial fit instead of the higher
order spline.  You can do linear interpolation if you want, but that is
always the worst case approach.