Cloudy & Associates

Commit 6bc75456 authored by Chatzikos, Marios's avatar Chatzikos, Marios
Browse files

Use NIST ionization potentials, not phfit.dat

Previously, the code used function ph1() to access the valence
threshold energies (ionization potentials) stored in file phfit.dat.

Define method atmdat.getIonPot() to access the NIST IPs (read from file
ionization_potentials.dat), and use it in lieu of the call to ph1(),
where possible.
parent 69c3fa58
......@@ -281,6 +281,17 @@ struct t_atmdat : public module {
/** accurate ionization potentials in Ryd */
double EIonPot[LIMELM][LIMELM];
/** getIonPot -- get ionization potential, in Ryd
*
* \param nelem [in] element index (0 for H, 29 for Zn)
* \param ion [in] ion charge (0 for neutral, nelem for H-like)
* \return double ionization potential, in Ry
*/
inline double getIonPot( long int nelem, long int ion )
{
return EIonPot[nelem][ion];
}
/** CharExcIon is ionization, */
/** [0] is Atom^0 + H+ => Atom+1 + H0
* [n] is Atom^+n + H+ => Atom^+n-1 + H0 */
......
......@@ -335,7 +335,7 @@ void ContCreatePointers(void)
/* check that ionization potential of neutral carbon valence shell is
* positive */
ASSERT( t_ADfA::Inst().ph1(2,5,5,0) > 0. );
ASSERT( atmdat.getIonPot(ipCARBON, 0) > 0. );
/* now fill in all sub-shell ionization array indices for elements heavier than He,
* this must be done after previous loop on iso.ipIsoLevNIonCon[ipH_LIKE] since hydrogenic species use
......@@ -348,13 +348,13 @@ void ContCreatePointers(void)
}
/* most of these are set in ipShells, but not h-like or he-like, so do these here */
Heavy.Valence_IP_Ryd[ipHYDROGEN][0] = t_ADfA::Inst().ph1(0,0,ipHYDROGEN,0)/EVRYD;
Heavy.Valence_IP_Ryd[ipHELIUM][0] = t_ADfA::Inst().ph1(0,1,ipHELIUM,0)/EVRYD;
Heavy.Valence_IP_Ryd[ipHELIUM][1] = t_ADfA::Inst().ph1(0,0,ipHELIUM,0)/EVRYD;
Heavy.Valence_IP_Ryd[ipHYDROGEN][0] = atmdat.getIonPot(ipHYDROGEN, 0);
Heavy.Valence_IP_Ryd[ipHELIUM][0] = atmdat.getIonPot(ipHELIUM, 0);
Heavy.Valence_IP_Ryd[ipHELIUM][1] = atmdat.getIonPot(ipHELIUM, 1);
for( long nelem=2; nelem<LIMELM; ++nelem )
{
Heavy.Valence_IP_Ryd[nelem][nelem-1] = t_ADfA::Inst().ph1(0,1,nelem,0)/EVRYD;
Heavy.Valence_IP_Ryd[nelem][nelem] = t_ADfA::Inst().ph1(0,0,nelem,0)/EVRYD;
Heavy.Valence_IP_Ryd[nelem][nelem-1] = atmdat.getIonPot(nelem, nelem-1);
Heavy.Valence_IP_Ryd[nelem][nelem] = atmdat.getIonPot(nelem, nelem);
if( dense.lgElmtOn[nelem])
{
/* now confirm that all are properly set */
......
......@@ -142,7 +142,7 @@ double hydro_energy(long nelem, long n, long /* l */, long /* s */, long /* g *
else
/* Dima's data in phfit.dat have ionization potentials in eV
* with four significant figures*/
HIonPoten = t_ADfA::Inst().ph1(0,0,nelem,0)/EVRYD;
HIonPoten = atmdat.getIonPot(nelem, nelem)/EVRYD;
ASSERT(HIonPoten > 0.);
return HIonPoten/POW2((double)n)*RYD_INF;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment