Extend use of NIST ion. pot. to photo CS

Previously, NIST ionization potentials were not used when computing
photoionization cross sections; instead the data in phfit.dat were used.

Define method getEthresh to switch between using NIST values for valence
shells, and the Verner data for inner shell ionization thresholds.

source/atmdat_adfa.cpp

@@ -265,6 +265,26 @@ t_ADfA::t_ADfA()
 	}
 }
 
+inline long t_ADfA::set_vshell_index( long nelec, long Z ) const
+{
+	long nout = NTOT[nelec-1];
+	if( Z == nelec && Z > 18 )
+		nout = 7;
+	if( Z == (nelec + 1)
+	    && (Z == 20 || Z == 21 || Z == 22 || Z == 25 || Z == 26) )
+		nout = 7;
+	return nout;
+}
+
+double t_ADfA::getEthresh( long int nshell, long int nelec, long int Z ) const
+{
+	long nout = set_vshell_index( nelec, Z );
+	if( nshell == nout )
+		return atmdat.getIonPot(Z-1, Z-nelec) * EVRYD;
+	else
+		return double(PH1[nshell-1][nelec-1][Z-1][0]);
+}
+
 double t_ADfA::phfit(long int nz, 
 		     long int ne,
 		     long int is, 
@@ -322,12 +342,8 @@ double t_ADfA::phfit(long int nz,
 		return crs;
 	}
 
-	nout = NTOT[ne-1];
-	if( nz == ne && nz > 18 )
-		nout = 7;
-	if( nz == (ne + 1) && ((((nz == 20 || nz == 21) || nz == 22) || 
-	  nz == 25) || nz == 26) )
-		nout = 7;
+	nout = set_vshell_index( ne, nz );
+
 	if( is > nout )
 	{ 
 		return crs;
@@ -340,8 +356,10 @@ double t_ADfA::phfit(long int nz,
 
 	ASSERT( is >= 1 && is <= 7 );
 
-	if( e < PH1[is-1][ne-1][nz-1][0] )
-	{ 
+	double ethres = getEthresh( is, ne, nz );
+
+	if( e < ethres )
+	{
 		return crs;
 	}
 
@@ -358,7 +376,10 @@ double t_ADfA::phfit(long int nz,
 		}
 		else
 		{
-			einn = PH1[nint-1][ne-1][nz-1][0];
+			if( nint == nout )
+				einn = ethres;
+			else
+				einn = double(PH1[nint-1][ne-1][nz-1][0]);
 		}
 	}
 
@@ -451,7 +472,7 @@ double t_ADfA::hpfit(long int iz,
 		}
 	}
 
-	eth = ph1(0,0,iz-1,0)/POW2((double)m);
+	eth = getEthresh(1,1,iz)/POW2((double)m);
 	ex = MAX2(1. , e/eth );
 
 	/* Don't just force to be at least one...make sure e/eth is close to one or greater.	*/

source/atmdat_adfa.h

@@ -23,6 +23,18 @@ private:
 	long int NTOT[LIMELM];
 	realnum PH1[7][LIMELM][LIMELM][6];
 	realnum PH2[LIMELM][LIMELM][7];
+
+	/** set_vshell_index -- Set valence shell index for given ion
+	 * The index is used to access the photoionization cross-section
+	 * data of Verner et al. (1995) or (1996).
+	 *
+	 * \param [in] ne   number of electrons
+	 * \param [in] Z    atomic number of element
+	 *
+	 * \returns long    index of valence shell
+	 */
+	inline long set_vshell_index( long ne, long Z ) const;
+
 	/* hpfit.dat */
 	realnum PHH[NHYDRO_MAX_LEVEL][5];
 	/* rec_lines.dat */
@@ -63,6 +75,17 @@ public:
 	*/
 	realnum ph1(int i, int j, int k, int l) const { return PH1[i][j][k][l]; }
 
+	/** getEthresh -- get ionization threshold for shell
+	 * For inner shells use data of Verner et al. (1995) or (1996);
+	 * for valence shells use NIST ionization potentials.
+	 *
+	 * \param nshell [in]  shell index (1 - 1s, 2 - 2s, etc, <= 7)
+	 * \param nel    [in]  number of electrons in ion
+	 * \param Z      [in]  atomic number of element
+	 * \return ionization threshold, in eV
+	 */
+	double getEthresh( long nshell, long nel, long Z ) const;
+
 	/** sth array of cross sections for photoionization of hydrogen at threshold,
 	    0 is 1s, 1 is 2s, 2 is 2p, up to 400
 	    \param i

source/cont_createpointers.cpp

@@ -997,7 +997,7 @@ STATIC void ipShells(
 		for( nshell=0; nshell < imax; nshell++ )
 		{
 			/* ionization potential of this shell in rydbergs */
-			thresh = (double)(t_ADfA::Inst().ph1(nshell,nelec-1,nelem,0)/EVRYD);
+			thresh = t_ADfA::Inst().getEthresh(nshell+1,nelec,nelem+1)/EVRYD;
 			if( thresh <= 0.1 )
 			{
 				/* negative ip shell does not exist, set upper limit
@@ -1074,7 +1074,7 @@ STATIC void ipShells(
 	/* this statement is in ContCreatePointers but has not been done when this routine called */
 	/*iso_sp[ipH_LIKE][ipZ].fb[ipLo].ipIsoLevNIonCon = ipContEnergy(iso_sp[ipH_LIKE][ipZ].fb[ipLo].xIsoLevNIonRyd,chLab);*/
 	/*opac.ipElement[nelem][nelem][0][0] = iso_sp[ipH_LIKE][nelem].fb[ipH1s].ipIsoLevNIonCon;*/
-	opac.ipElement[nelem][nelem][0][0] = ipoint( t_ADfA::Inst().ph1(0,0,nelem,0)/EVRYD );
+	opac.ipElement[nelem][nelem][0][0] = ipoint( t_ADfA::Inst().getEthresh(1,1,nelem+1)/EVRYD );
 	ASSERT( opac.ipElement[nelem][nelem][0][0] > 0 );
 
 	/* this is the high-energy limit */

source/opacity_add1element.cpp

@@ -27,7 +27,7 @@ void OpacityAdd1Element(
 
 	/* this routine drives OpacityAdd1Subshell to put in total opacities for all shells*/
 
-	/*begin sanity check */
+	/* begin sanity check */
 	ASSERT( (nelem >=0 ) && (nelem < LIMELM) );
 
 	/* first do simple two-level systems -

source/opacity_add1subshell.cpp

@@ -101,7 +101,7 @@ void OpacityAdd1SubshellInduc(
 
 	k = ipOpac - ipLowEnergy;
 
-	/* DepartCoef is dep coef, rfield.lgInducProcess is turned off with 'no indcued' command */
+	/* DepartCoef is dep coef, rfield.lgInducProcess is turned off with 'no induced' command */
 	if( (DepartCoef > 1e-35 && rfield.lgInducProcess) && hydro.lgHInducImp )
 	{
 		iup = MIN2(ipHiEnergy,rfield.nflux);

source/opacity_createall.cpp

@@ -748,7 +748,7 @@ STATIC void OpacityCreate1Element(
 			{
 				/* photo energy MAX so that we never eval below threshold */
 				energy = MAX2(rfield.anu(ip)*EVRYD , 
-					t_ADfA::Inst().ph1(nshell,nelec-1,nelem,0));
+					t_ADfA::Inst().getEthresh(nshell+1,nelec,nelem+1));
 
 				/* the cross section in mega barns */
 				cs = t_ADfA::Inst().phfit(nelem+1,nelec,nshell+1,energy);
@@ -803,9 +803,7 @@ STATIC double Opacity_iso_photo_cs(
 	{
 		if( index==0 )
 		{
-			/* this is the ground state, use Dima's routine, which works in eV
-			 * and returns megabarns */
-			double EgammaEV = MAX2(EgammaRyd*(realnum)EVRYD , t_ADfA::Inst().ph1(0,0,nelem,0));
+			double EgammaEV = MAX2(EgammaRyd, atmdat.getIonPot(nelem, nelem)) * EVRYD;
 			crs = t_ADfA::Inst().phfit(nelem+1,1,1,EgammaEV)* 1e-18;
 			/* make sure cross section is reasonable */
 			ASSERT( crs > 0. && crs < 1e-10 );

source/pressure_total.cpp

@@ -143,7 +143,7 @@ void PresTotCurrent()
 				/* this is the sum of all the energy needed to bring the atom up
 				 * to the ion+1 level of ionization - at this point a positive number */
 				phycon.EnergyIonization += dense.xIonDense[nelem][ion] * 
-					t_ADfA::Inst().ph1(Heavy.nsShells[nelem][i-1]-1,nelec,nelem,0)/EVRYD*kadvec;
+					t_ADfA::Inst().getEthresh(Heavy.nsShells[nelem][i-1],nelec+1,nelem+1)/EVRYD*kadvec;
 			}
 		}
 	}

source/save_do.cpp

@@ -2427,7 +2427,7 @@ void SaveDo(
 										break;
 
 									/* array elements are shell, numb of electrons, element, 0 */
-									energy = t_ADfA::Inst().ph1(ips,nelem-ion,nelem,0);
+									energy = t_ADfA::Inst().getEthresh(ips+1,nelem-ion+1,nelem+1);
 
 									/* now print threshold with correct format */
 									if( energy < 10. )

source/save_opacity.cpp

@@ -203,7 +203,7 @@ void save_opacity(FILE * ioPUN,
 
 						/* ionization potential of shell */
 						fprintf(ioPUN,"\t%.2f" ,
-							t_ADfA::Inst().ph1(nshell,nelec-1,nelem,0)/EVRYD);
+							t_ADfA::Inst().getEthresh(nshell+1,nelec,nelem+1)/EVRYD);
 
 						/* set lower and upper limits to this range */
 						long ipop = opac.ipElement[nelem][ion][nshell][2];