Simon Jeffery's Software Store
TAP: equation of state
Quick Reference
| SUBROUTINE TAP_SAHA ( ZELEM, TEMP,
NEL, ION_MIN, ION_MAX, U, CHI, IONFRN, STATUS ) |
Calculate relative abundances of ions of species K. |
| FUNCTION TAP_PARFN ( ZELEM, ZION,
TEMP, PEL ) |
Calculate the partition function for a given ion,
at a given temperature and electron pressure.
|
| FUNCTION TAP_IONPT ( ZELEM, ZION
) |
Provides ionisation potential for a given ion.
|
Long Reference
DOUBLE PRECISION FUNCTION TAP_IONPT ( ZELEM, ZION )
Purpose:
Provides ionisation potential for a given ion
Invocation:
RESULT = TAP_IONPT( atom_z, ion_z, status )
Description:
This function provides ionisation potentials expressed in electron
volts. Recommended values are derived from analyses of optical
spectra. The following sources have been used, in descending priority
1: Moore 1970, table 1
2: Allen 1973, table 17 (electron affinities)
3. Allen 1973, table 16
The ionisation potentials given by Moore were derived by multiplying
the observed series limits by 0.000123981, i.e., by using the
conversion factor
1eV = 806573. m**(-1)
The 1986 CODATA (Cohen & Taylor 1987) recommended value for this
quantity is
1eV = 806544.10 m**(-1)
+/- .24
In order to preserve the form of the data tabulated by Moore,
and thus to to allow easy verification, but also to be achieve the
highest possible standard of accuracy, this version applies
the correction at execution time. With 1986 CODATA values, the
difference amounts to about 36 parts per million, as anticipated
by Moore.
Arguments:
ZELEM = INTEGER (Given)
The atomic number of the ion ( eg C: = 6)
ZION = INTEGER (Given)
The atomic charge on the ion ( eg CI: = 0, OVI: = 5 )
Returned Value:
TAP_IONPT = DOUBLE PRECISION
Ionisation potential in electron volts (eV)
References:
- Moore C.E., 1970. NSRDS - NBS 34 "Ionisation Potentials and
Ionization Limits Derived from the Analyses of Optical Spectra"
- Cohen E.R. & Taylor B.N., 1987. Journal of Research of the
National Bureau of Standards, 92, 85.
- Allen C.W.. 1973. Astrophysical Quantities, 3rd Edition
atlas9
- Lester version of Kurucz code in CCP7 library
Accuracy:
The accuracy of the data stored reflects the number of
significant figures given by the sources, such that the final
significant figure is believed to be meaningful.
Implementation Deficiencies:
The uncertainty of 30 ppm in the conversion factor used by Moore
has been eliminated by applying a correction factor. This process
arbitrarily increases the number of "significant figures" in the
returned value cf. the stored value, and has the effect of
slightly increasing the execution time. These effects should be
removed in a future version.
In order to allow ionizations to be calculated for rare-earths,
a third and sometimes even a second ionization potential has been
adopted as in atlas9. Cases are marked - #1 -
Pitfalls:
If unusual results are obtained, check that the ion identifier has
been entered correctly: i.e., 0 for neutrals, 1 for singly
ionized, etc.. Do not use the frequently used "spectrum" identifier,
i.e. I for neutrals, II for singly ionized, etc.
Error conditions:
The status flag will be set if any of the following conditions are
encountered:
an illegal value for the atomic number ( zelem )
an illegal value for the ion charge ( zion )
a datum is unavailable
In all cases, the result of the function call will be 0 eV.
Bugs:
Data from Allen for higher ionization levels than given by Moore
has not yet been entered.
SUBROUTINE TAP_SAHA ( ZELEM, TEMP, NEL, ION_MIN, ION_MAX, U, CHI, IONFRN, STATUS )
Purpose:
Calculate relative abundances of ions of species K.
Language:
Fortan 90
Invocation:
CALL TAP_SAHA ( zelem, temp, nel, ion_min, ion_max, u, chi,
: ion_frn, status )
Description:
Calculate relative ionisation fractions of ions j = 1,.... for
atomic species (Z), given temperature and electron density,
using the Saha ionisation equation. By recursive application of
the usual Saha equation relating relative abundances of successive
ionisation stages (Njk/Nj+1k), a general equation giving the
fraction Njk/Nk is obtained (Mihalas eqn. 5-17). This is
evaluated using scaled logarithmic variables.
The programmer must supply the atomic number of the species to be
treated, the gas temperature (K) and electron density (m**-3).
He is also required to provide the partition functions and
ionisation potentials for each ionisation stage to be considered.
These requirements allow for some flexibility in the way the
subroutine may be used, and for the inclusion of additional
phsyics (e.g. pressure ionisation, electron shielding, etc.)
By specifying ION_MIN and ION_MAX, the programmer may treat only
selected ionisation stages (saving unnecessary computations).
By using an internal OFFSET, the locations 1,2,... in the partition
function, ionisation potential and returned ionisation fraction
arrays map onto the ionisation stages ION_MIN ... ION_MAX.
Consider two examples:
i) calculate ion fractions for neutrals through triply ionised
ii) calculate ion fractions for -ve ions through doubly ionized
Arrays U, CHI and ION should have data in the locations 1 onwards
example (i) example (ii)
U 0 1 2 3 -1 0 1 2
CHI 0 1 2 -1 0 1
IONFRN 0 1 2 3 -1 0 1 2
Arguments:
ZELEM = INTEGER (Given)
Atomic number (ie hydrogen = 1, helium = 2, etc)
TEMP = DOUBLE PRECISION (Given)
Temperature (degrees K)
NEL = DOUBLE PRECISION (Given)
Electron density (m**-3)
ION_MIN = INTEGER (Given)
Minimum ionic charge to be treated (ie negative = -1, ...)
ION_MAX = INTEGER (Given)
Maximum ionic charge to be treated (ie neutral = 0, ...)
U = VARIABLE LENGTH DOUBLE PRECISION ARRAY (Given)
Partition function for each ionisation stage
CHI = VARIABLE LENGTH DOUBLE PRECISION ARRAY (Given)
Ionisation potential for each ionisation stage
IONFRN = VARIABLE LENGTH DOUBLE PRECISION ARRAY (Returned)
Array containing fraction in each ionisation stage
STATUS = INTEGER (Returned)
The global status.
References:
Mihalas D., 1978, Stellar Atmospheres (2nd ed.) Freeman
Bugs:
Does not treat negative ions.
{note_any_bugs_here}
DOUBLE PRECISION FUNCTION TAP_PARFN ( ZELEM, ZION, TEMP, PEL )
Purpose:
Calculate the partition function for a given ion,
at a given temperature and electron pressure
Language:
Fortran 90
Invocation:
U = TAP_PARFN ( zelem, zion, temp, pel )
Description:
Calculate partition functions using approximation formulae
due to Traving et al. (1966a,b) and Unsoeld (1948).
If these formulae have not been included (coded) then
a ground-state partition function is returned.
Arguments:
ZELEM = INTEGER (Given)
Atomic number (ie hydrogen = 1, helium = 2, etc )
ZION = INTEGER (Given)
Atomic charge (ie neutral = 0, singly ionised = 1, etc )
TEMP = DOUBLE PRECISION (Given)
Temperature (degrees K)
PEL = DOUBLE PRECISION (Given)
Electron pressure (N m**-2)
Returned Value:
TAP_PARFN = DOUBLE PRECISION
Partition function U(T,Pe)
Algorithm:
U = g0 + Sum(i=1,ns)[U'(i)]
U' = Sum(j=1,nk)[g(j).10^(-Chi(j).theta)]
+ 2.gpr.Qas(l,z).10^(-Chi(ion).theta)
theta = 5040 / T
Chebyshev approximation
f(theta) = Sum(j=1,nk)[g(j).10^(-Chi(j).theta)]
~ psi(theta) = Sum(k=1,nm)[alpha(k).10^(-gamma(k).theta)]
Accuracy:
Chebyshev coefficients are valid for 0 <= Chi.theta <= 800
In most cases, the maximum error at max(theta) is limited to
|f(theta) - psi(theta)| < 0.01
Limitations:
If Chi.theta > 800, the ground state statistical weight is returned.
Where Chebyshev polynomials are not yet entered,
the ground state statistical weight is returned.
If ground state statistical weights are unavailable
a default value U = 1.0 is returned.
Species currently included as a Chbyshev poynomial
(indicated by ion identifier):
ZELEM ZION + I
1 H I
2 He I II
3 Li
4 Be
5 B
6 C I II III IV V
7 N I II III IV V
8 O I II III IV
9 F
10 Ne I II III IV
11 Na
12 Mg I II III IV
13 Al I II III IV
14 Si I II III IV V
15 P I II III IV
16 S I II III IV
17 Cl
18 A I II III IV
19 K
20 Ca
21 Sc
22 Ti
23 V
24 Cr
25 Mn
26 Fe I II III IV
27 Co
28 Ni
29 Cu
30 Zn
31 Ga
32 Ge
33 As
34 Se
35 Br
36 Kr
37 Rb
38 Sr I II III IV
39 Y I II III IV
40 Zr I II III IV
Note:
An accurate expression for THETA is given by
THETA = EV_K / TEMP
but in this subroutine it is more appropriate to use the
expression used by Traving et al. (1966a)
THETA = 5040. / TEMP
References:
Traving et al. 1966a. Abh.der.Hamburg.Sternw. 8,3
Traving et al. 1966b. Abh.der.Hamburg.Sternw. 8,26
Unsold, 1948. Zeit.f.Astrophysik. 24,356.
External Routines Used:
TAP:
TAP_PAR, [routine_used]...
[optional_subroutine_items]...
Authors:
UH: U. Heber (Kiel)
CSJ: C.S.Jeffery (St Andrews/CCP7)
RED: R.E.Dudley (St Andrews) : Chebychev polynomial data
RWJR: R.W.J.Rolleston (QUB) : Chebychev polynomial data
CSJ: C.S.Jeffery (St Andrews/CCP7) : Version for TAP
{enter_new_authors_here}
History:
12-MAR-1992 (CSJ):
TAP version.
29-FEB-1993 (CSJ):
Sr, Y and Zr added
9-MAY-1994 (CSJ):
Converted to double precision function for
consistency with other TAP subprograms.
29-SEP-2000 (CSJ)
gsw's added for Z > 40.
{enter_changes_here}
Bugs:
Ground state statistical weights provided for first three ionisation
states only, for Z >= Sc.
No ground state statistical weights for As and Br, bad values returned.
No ground state statistical weights for negative ions.
For As, Br and Z>40, gsw's are taken from Kurucz code
(PFSAHA).
There are probably thousands of app. formulae for
pf's in the literature, the OP are known to have done
their own thing as well. This code really need cleaning
up and making completely general - for given conditions.
{note_any_bugs_here}