Diatomic partition functions in IDL

Recently I posted a set of IDL routines for the atomic partition functions (PF) based on various polynomial fits that are available in the literature. Here I do the same for the diatomic molecules.

The routines for the diatomic PF use an internal catalog of 325 diatomic molecules. The ionized molecules are separate entries. The catalog is stored as an IDL structure in catalogue_of_diatomics.sav. The structure has the following tags for every molecule:

Name	Chemical formula, e.g. 'H2', 'CO', 'H2+'
Nc	Number of constituents, always 2 for diatomics
Nuc	Number of uniq constituents, e.g. 2 for CO, 1 for H2
Constituents	Constituent atoms, e.g. ['H', 'H'], ['C', 'O']
Charge	Charge, e.g. -1, 0, +1
Code	Kurucz's code for molecules, e.g. CO '0608.00'
D0	Dissociation energy in eV
Data	List of internally available data sets for this molecule
Type	Type of data: pf for partition function, kp for chemical equilibrium constant, eint for internal energy (one entry for each dataset)
References	List of ADS references (for each dataset)

Input of the routine are the temperature (in Kelvin) and the chemical formula of the molecule (e.g. 'H2', 'CO'). The output is specified by the keywords DATA and TYPE. Depending on TYPE, the output can be the partition function ('pf'), the chemical equilibrium constant for pressure ('kp') or the internal energy due to the rotational-vibrational bound states ('eint'). The chemical equilibrium constant is given for the pressures (to get its value for the number density divide it by kT). The internal energy in erg per particle. (All the units in the output and used internally in the code are in CGS).

This code uses various polynomial fits available in the literature and evaluates them for the given temperature and the specified molecule. The keyword DATA specifies the source from which the data is taken and the particular fit. The following sources are available:

Reference	Function	Molecules	DATA keyword value
Vardya 1961	Kp	H2 and H2+	'vardya'
Vardya 1965	Eint	H2 and H2+	'vardya'
Tsuji 1973	U	95 diatomics	'tsuji1973'
Sauval and Tatum 1984	U, Kp, D0	294 diatomics	'sauvalandtatum1984'
Irwin 1981	U	H2 and CO	'irwin1981'
Bohn and Wolf 1984	U	H2 and CO	'bohn1984'

None of the original datasets contain all three functions. The missing functions are computed from the present ones. The chemical equilibrium constant:

$$K_p = \frac{Ua\,Ub}{Uab} \left(\frac{2\pi m_{\mathrm{ab}} kT}{h^2}\right)^\frac{3}{2} \mathrm{e}^{-\frac{D_0}{kT}} ,$$
and the internal energy:
$$ E_{\mathrm{int}} = kT \frac{\partial \ln U}{\partial \ln T},$$
where $U_{\mathrm{a}}$ and $U_{\mathrm{b}}$ are the partition functions of the atoms A and B, $U_{\mathrm{ab}}$ is the partition function of the diatomic, $m_{\mathrm{ab}}$ is the reduced mass ($m_{\mathrm{ab}} = m_{\mathrm{a}}*m_{\mathrm{b}}/(m_{\mathrm{a}}+m_{\mathrm{b}})*m_{\mathrm{h}}$) and $D_0$ is the dissociation energy.

The default data is 'sauvalandtatum1984' for all the molecules that are included in their list and 'tsuji1973' for the rest.

The dissociation constants are available only for the molecules of Sauval and Tatum. Therefore the chemical equilibrium constant cannot be computed for the molecules that are exclusiv for Tsuji if the user does not supply $D_0$ manually.

Note that each of these data sources strictly applies only to a limited range of temperatures: Gray (2520 - 25200 K), Irwin (1000 - 16000 K), Sauval and Tatum (1000 - 9000 K). Out of these limits the values returned by this function may not be reliable and it is up to the user to decide on the extrapolation procedure in that case. If keyword RANGE is set in a call to a partition function, then the temperatures out of the range are clipped to it.

Notes:

•  The set of Irwin (1981) for the diatomics replicates an earlier set of Tatum (1966) that is already included and updated in Sauval and Tatum (1984). 
•  There is also a set of Kurucz that can be easily added to the catalogue, but it is not supported yet.
•  The data of Irwin had been updated since 1987, but I do not include them here as I miss the full reference and the data description.

If you use any of the data in a publication, please acknowledge it by citing the corresponding paper.


Routines

The tarball is here, diatomics.tar. The individual functions and the data files are:

diatomics.pro
diatomics_irwin.pro
diatomics_sauval_and_tatum.pro
diatomics_vardya.pro
diatomics_tsuji.pro
diatomics_bohn.pro
catalogue_of_diatomics.sav
sauvalandtatum1984_pf.sav
sauvalandtatum1984_eqc.sav
tsuji1973_eqc.sav


Installation

Before running the code for the first time, copy all the functions into a directory within your IDL path. Open diatomics_path.pro and set the path to that exact directory. That's all.


Example

Let's compare PFs from different sources for the neutral H2:

IDL> name = 'H2'
IDL> t = 2000. + DINDGEN(6)*2000.
IDL> pfi = DIATOMICS(t, z, data = 'i')
IDL> pfs = DIATOMICS(t, z, data = 's')
IDL> pfb = DIATOMICS(t, z, data = 'b') 

And let's do the same for the chemical equilibrium constants:

IDL> name = 'H2'
IDL> t = 2000. + DINDGEN(6)*2000.
IDL> kpi = DIATOMICS(t, z, data = 'i')
IDL> kps = DIATOMICS(t, z, data = 's')
IDL> kpv = DIATOMICS(t, z, data = 'v')
IDL> kpt = DIATOMICS(t, z, data = 't')
IDL> kpb = DIATOMICS(t, z, data = 'b') 

UPDATE (Feb.2017)
I found a mistake in catalogue_of_diatomics.sav. Two molecules were labeled as CoF (168 and 179, starting from 0). The first one is actually CuF, the second CdF. The corrected version of the catalogue is uploaded.