/* This file is part of Cloudy and is copyright (C) 1978-2004 by Gary J. Ferland.
 * For conditions of distribution and use, see copyright notice in license.txt */
/*he_1trans compute Aul for given line	*/
/*ritoa - converts the square of the radial integral for a transition 
 * (calculated by scqdri) to the transition probability, Aul.	*/
/*ForbiddenAuls calculates transition probabilities for forbidden transitions.	*/
/*scqdri - stands for Semi-Classical Quantum Defect Radial Integral	*/
/*Jint - used by scqdri	*/
/*AngerJ - used by scqdri */
/*DoFSMixing - applies a fine structure mixing approximation to A's.  To be replaced by 
 * method that treats the entire rate matrix.	*/

#include "cddefines.h" 
#include "physconst.h" 
#include "taulines.h"
#include "path.h"
#include "dense.h"
#include "trace.h"
#include "hydro_bauman.h"
#include "iso.h"
#include "helike.h"
#include "helike_einsta.h"

/*lint -e662 creation of  out of bound pointer */
/*lint -e661 creation of  out of bound pointer */

/* the array of transitions probabilities read from data file.  */
static double ***TransProbs;

static double ritoa( long li, long lf, long nelem, double k, double RI2 );

static double ForbiddenAuls( long ipHi, long ipLo, long nelem );

static double CalculateRadialIntegral( long nelem, long ipLo, long ipHi );

static double RadialIntegrand( double r12 );

static double WaveFunction( double r12 );

static long RI_ipHi, RI_ipLo, RI_ipLev;
static long globalZ;

/*static FILE *ofp;  */

/* used as parameters in qg32 integration */
static double vJint , zJint;

static double Jint( double theta )
{
	/*	[ cos[vx - z sin[x]] ]  */
	double 
		d0 = ( 1.0 / PI ),
		d1 = vJint * theta,
		d2 = zJint * sin(theta),
		d3 = (d1 - d2),
		d4 = cos(d3),
		d5 = (d0 * d4);

	return( d5 );
}

static double AngerJ( double vv, double zz )
{
	long int rep = 0, ddiv, divsor;

	double y = 0.0;

	/* Estimate number of peaks in integrand.  */                            
	/* Divide region of integration by number  */
	/*  peaks in region.                       */
	if( (fabs(vv)) - (int)(fabs(vv)) > 0.5 )
		ddiv = (int)(fabs(vv)) + 1;
	else 
		ddiv = (int)(fabs(vv));

	divsor  = ((ddiv == 0) ? 1 : ddiv);
	vJint = vv;
	zJint = zz;

	for( rep = 0; rep < divsor; rep++ )
	{
		double
		rl = (((double) rep)/((double) divsor)),
		ru = (((double) (rep+1))/((double) divsor)),
		x_low = (PI * rl),
		x_up  = (PI * ru);      

		y += qg32( x_low, x_up, Jint );
	}

	return( y );
}

/******************************************************************************/
/******************************************************************************/
/*                                                                            */
/*    Semi-Classical Quantum Defect Radial Integral                           */      
/*                                                                            */
/*   See for example                                                          */
/*     Atomic, Molecular & Optical Physics Handbook                           */
/*     Gordon W. F. Drake; Editor                                             */
/*     AIP Press                                                              */
/*     Woddbury, New York.                                                    */
/*     1996                                                                   */
/*                                                                            */
/* NOTE:: we do not include the Bohr Radius a_o in the                        */
/*           definition of of  R(n,L;n'L') as per Drake.                      */
/*                                                                            */
/*                                                                            */
/*                   1  (n_c)^2 | {      D_l max(L,L') }                      */
/*    R(n,L;n'L') = --- ------- | { 1 -  ------------- } J( D_n-1 ; -x )  -   */
/*                   Z   2 D_n  | {          n_c       }                      */
/*                                                                            */
/*                                                                            */
/*                    {      D_L max(L,L') }                                  */
/*                -   { 1 +  ------------- } J( D_n+1 ; -x )                  */
/*                    {          n_c       }                                  */
/*                                                                            */
/*                                                                            */
/*                     2                    |                                 */
/*                  + ---  sin(Pi D_n)(1-e) |                                 */
/*                     Pi                   |                                 */
/*                                                                            */
/*  where                                                                     */
/*        n_c = (2n*'n*)/(n*'+n*)                                             */
/*                                                                            */
/*       Here is the quantity Drake gives...                                  */
/*            n_c = ( 2.0 * nstar * npstar ) / ( nstar + npstar ) ;           */
/*                                                                            */
/*       while V.A. Davidkin uses                                             */
/*            n_c = sqrt(  nstar * npstar  );                                 */
/*                                                                            */
/*        D_n = n*' - n*                                                      */
/*                                                                            */      
/*        D_L = L*' - L*                                                      */
/*                                                                            */
/*        x = e D_n                                                           */
/*                                                                            */
/*        Lmx  = max(L',L)                                                    */
/*                                                                            */
/*        e = sqrt( 1 - {Lmx/n_c}^2 )                                         */
/*                                                                            */
/*                                                                            */
/*       Here n* = n - qd where qd is the quantum defect                      */
/*                                                                            */
/******************************************************************************/
/******************************************************************************/
double scqdri(/* upper and lower quantum numbers...n's are effective	*/
			  double nstar, long int l,
			  double npstar, long int lp,
              double iz
              )
{
	double n_c = ((2.0 * nstar * npstar ) / ( nstar + npstar ));
	/*double n_c = sqrt( nstar * npstar );*/
	
	double D_n = (nstar - npstar);
	double D_l = (double) ( l - lp );
	double lg  = (double) ( (lp > l) ? lp : l);

	double h = (lg/n_c);
	double g = h*h;
	double f = ( 1.0 - g );
	double e = (( f >= 0.0) ? sqrt( f ) : 0.0 );

	double x  = (e * D_n);
	double z  = (-1.0 * x);
	double v1 = (D_n + 1.0);
	double v2 = (D_n - 1.0); 

	double d1,d2,d7,d8,d9,d34,d56,d6_1; 

	if ( iz == 0.0 )
		iz += 1.0;

	if ( D_n == 0.0 )
	{
	  return( -1.0 );
	}

	if ( D_n < 0.0 )
	{
	  return( -1.0 );
	}

	if ( f < 0.0 )
	{
	  /* This can happen for certain  quantum defects   */
	  /* in the lower n=1:l=0 state. In which case you  */
	  /* probably should be using some other alogrithm  */
	  /* or theory to calculate the dipole moment.      */
	  return( -1.0 );
	}

	d1 = ( 1.0 / iz );

	d2 = (n_c * n_c)/(2.0 * D_n);

	d34 = (1.0 - ((D_l * lg)/n_c)) * AngerJ( v1, z );

	d56 = (1.0 + ((D_l * lg)/n_c)) * AngerJ( v2, z );

	d6_1 = PI * D_n;

	d7 = (2./PI) * sin( d6_1 ) * (1.0 - e);

	d8 = d1 * d2 * ( (d34) - (d56) + d7 );

	d9 = d8 * d8;

	ASSERT( D_n  > 0.0 );
	ASSERT( l  >= 0  );
	ASSERT( lp >= 0 );
	ASSERT( (l == lp + 1) || ( l == lp - 1) );
	ASSERT( n_c != 0.0 );
	ASSERT( f >= 0.0 );
	ASSERT( d9  > 0.0 );

	return( d9 );
}

static double ForbiddenAuls( long ipHi, long ipLo, long nelem )
{
	double A;
	/* >>refer	Helike	2pho	Derevianko, A., & Johnson, W.R. 1997, Phys. Rev. A 56, 1288
	 * numbers are not explicitly given in this paper for Z=21-24,26,27,and 29.
	 * So numbers given here are interpolated.	*/
	double As2nuFrom1S[28] = {1940.,1.82E+04,9.21E+04,3.30E+05,9.44E+05,2.31E+06,5.03E+06,1.00E+07,
		1.86E+07,3.25E+07,5.42E+07,8.69E+07,1.34E+08,2.02E+08,2.96E+08,4.23E+08,5.93E+08,8.16E+08,
		1.08E+09,1.43E+09,1.88E+09,2.43E+09,3.25E+09,3.95E+09,4.96E+09,6.52E+09,7.62E+09,9.94E+09};
	/* Important clarification, according to Derevianko & Johnson (see ref above), 2^3S can decay
	 * to ground in one of two ways: through a two-photon process, or through a single-photon M1 decay,
	 * but the M1 rates are about 10^4 greater that the two-photon decays throughout the entire
	 * sequence.  Thus these numbers, are much weaker than the effective decay rate, but should probably
	 * be treated in as a two-photon decay at some point	*/
	double As2nuFrom3S[28] = {1.25E-06,5.53E-05,8.93E-04,8.05E-03,4.95E-02,2.33E-01,8.94E-01,2.95E+00,
		8.59E+00,2.26E+01,5.49E+01,1.24E+02,2.64E+02,5.33E+02,1.03E+03,1.91E+03,3.41E+03,5.91E+03,
		9.20E+03,1.50E+04,2.39E+04,3.72E+04,6.27E+04,8.57E+04,1.27E+05,2.04E+05,2.66E+05,4.17E+05};
				
	if ( (ipLo == ipHe1s1S) && (N_(ipHi) == 2) )
	{
		if( nelem == ipHELIUM )
		{
			/* All of these but the second and third one (values 51.02 and 1E-20) are from
			 * >>refer	HeI	As	Lach, G, & Pachucki, K, 2001, Phys. Rev. A 64, 042510
			,* 1E-20 is made up
			 * 51.3 is from the Derevianko & Johnson paper cited above.	*/
			double ForbiddenHe[5] = { 1.272E-4,	51.02,	1E-20,	177.58,	0.327 };

			A = ForbiddenHe[ipHi - 1];
			putError(nelem,ipHi,ipLo,ipRad,-1);
		}
		else
		{
			switch ( (int)ipHi )
			{
			case 1: /* Parameters for 2^3S to ground transition.	*/
				/* >>refer	Helike	As	Lin, C.D., Johnson, W.R., and Dalgarno, A. 1977, 
				 * >>refercon	Phys. Rev. A 15, 1, 010015	*/
				A = (3.9061E-7) * pow( (double)nelem+1., 10.419 ) + As2nuFrom3S[nelem-2];
				break;
			case 2: /* Parameters for 2^1S to ground transition.	*/
				A = As2nuFrom1S[nelem-2];
				break;
			case 3: /* Parameters for 2^3P0 to ground transition.	*/
				A = iso.SmallA;
				break;
			case 4: /* Parameters for 2^3P1 to ground transition.	*/
				A = ( 11.431 * pow((double)nelem, 9.091) );
				break;
			case 5: /* Parameters for 2^3P2 to ground transition.	*/
				/* According to the trend in Porquet & Dubau (2000), these are about 
				 * 1E-3 times the rate of the 2^3P1 to ground transition.	*/	
				A = ( 0.011431 * pow((double)nelem, 9.091) );
				break;
			default:
				TotalInsanity();
			}
			putError(nelem,ipHi,ipLo,ipRad,-1);
		}
		return A;
	}
	
	/* The next two cases are fits to probabilities given in 
	 * >>refer	He-like	As	Johnson, W.R., Savukov, I.M., Safronova, U.I., & 
	 * >>refercon	Dalgarno, A., 2002, ApJS 141, 543J	*/
	/* originally astro.ph. 0201454 */
	/* The involve Triplet P and Singlet S.  Rates for Triplet S to Singlet P 
	 * do not seem to be available.	*/
	
	/* Triplet P to Singlet S...Delta n not equal to zero!	*/
	else if( nelem>ipHELIUM && L_(ipHi)==1 && S_(ipHi)==1 && 
		L_(ipLo)==0 && S_(ipLo)==0 && N_(ipLo) < N_(ipHi) )
	{
		A = 8.0E-3 * exp(9.283/sqrt((double)N_(ipLo))) * pow((double)nelem,9.091) /
			pow((double)N_(ipHi),2.877);
		putError(nelem,ipHi,ipLo,ipRad,-1);
	}
	
	/* Singlet S to Triplet P...Delta n not equal to zero!	*/
	else if( nelem > ipHELIUM && L_(ipHi)==0 && S_(ipHi)==0 && 
		L_(ipLo)==1 && S_(ipLo)==1 && N_(ipLo) < N_(ipHi) )
	{
		A = 2.416 * exp(-0.256*N_(ipLo)) * pow((double)nelem,9.159) / pow((double)N_(ipHi),3.111);
		
		if ( ( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P2) ) )
		{
			/* This is divided by 3 instead of 9, because value calculated is specifically for 2^3P1.
			 * Here we assume statistical population of the other two.	*/
			A *= (2.*(ipLo-3)+1.0)/3.0;
		}
		putError(nelem,ipHi,ipLo,ipRad,-1);
	}
	
	else if( ( ipLo == ipHe2s3S ) && ( ipHi == ipHe2p1P ) )
	{
		/* This transition,1.549 , given by Lach and Pachucki, 2001 for the atom */
		if( nelem == ipHELIUM )
		{
			A = 1.549;
			putError(nelem,ipHi,ipLo,ipRad,-1);
		}
		else
		{
			/* This is a fit to data given in
			 * >>refer	He-like	As	Savukov, I.M., Johnson, W.R., & Safronova, U.I. 
			 * >>refercon	astro-ph 0205163	*/
			A= 0.1834*pow((double)nelem, 6.5735);
			putError(nelem,ipHi,ipLo,ipRad,-1);
		}
	}

	else
	{
		/* Current transition is not supported.	*/
		A = iso.SmallA;
		putError(nelem,ipHi,ipLo,ipRad,-1);
	}

	/* TODO - For now just just put 1% error for forbidden lines. */
	putError(nelem,ipHi,ipLo,ipRad,.01f);

	ASSERT( A > 0.);
	return A;   
}

/* Calculates Einstein A for a given transition.	*/
double he_1trans( 
			   /* charge on c scale, Energy is wavenumbers, Einstein A	*/
			   long nelem , double Enerwn ,
			   /* quantum numbers of upper level:	*/
			   double Eff_nupper, long lHi, long sHi, long jHi,
			   /* and of lower level: */
			   double Eff_nlower, long lLo, long sLo, long jLo )
			   /* Note j is only necessary for 2 triplet P...for all other n,l,s,
			    * j is completely ignored.	*/
{
	/* this will be A returned for totally forbidden transition 
#	define SMALLA 1e-20 - 
	replaced with iso.SmallA */

	double RI2, MyRI, Aul;
	long nHi, nLo, ipHi, ipLo;
	/*long lMax = MAX2( lLo, lHi );*/

	ASSERT(nelem > ipHYDROGEN);
	
	/* Since 0.4 is bigger than any defect, adding that to the effective principle quantum number,
	 * and truncating to an integer will produce the principal quantum number.	*/
	nHi = (int)(Eff_nupper + 0.4);
	nLo = (int)(Eff_nlower + 0.4);
	
	/* Make sure this worked correctly.	*/
	ASSERT( fabs(Eff_nupper-(double)nHi) < 0.4 );
	ASSERT( fabs(Eff_nlower-(double)nLo) < 0.4 );

	ipHi = QuantumNumbers2Index[nelem][nHi][lHi][sHi];
	if( (nHi==2) && (lHi==1) && (sHi==1) )
	{
		ASSERT( (jHi>=0) && (jHi<=2) );
		ipHi -= (2 - jHi);
	}

	ipLo = QuantumNumbers2Index[nelem][nLo][lLo][sLo];
	if( (nLo==2) && (lLo==1) && (sLo==1) )
	{
		ASSERT( (jLo>=0) && (jLo<=2) );
		ipLo -= (2 - jLo);
	}

	ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n == nHi );
	if( nHi <= iso.n_HighestResolved[ipHE_LIKE][nelem] )
	{
		ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ipHi].l == lHi );
		ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ipHi].s == sHi );
	}
	ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ipLo].n == nLo );
	if( nLo <= iso.n_HighestResolved[ipHE_LIKE][nelem] )
	{
		ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ipLo].l == lLo );
		ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ipLo].s == sLo );
	}

	/* First do allowed transitions	*/
	if ( (sHi == sLo) && (abs((int)(lHi - lLo)) == 1) )
	{
		Aul = -2.;

		/* For clarity, let's do this in two separate chunks...one for helium, one for everything else.	*/
		if ( nelem == ipHELIUM )
		{
			{
				/* Experiment with calculating radial integrals myself.	*/
				/*@-redef@*/
				enum {DEBUG_LOC=FALSE};
				/*@+redef@*/
		
				if( ( DEBUG_LOC ) && (nLo<=2) && (nHi==2) )
				{
					MyRI = CalculateRadialIntegral( ipHELIUM, ipLo, ipHi );
					fprintf(ioQQQ,"ipHi\t%li\tipLo\t%li\tMyRI^2\t%.2e\n",
						ipHi,
						ipLo,
						MyRI*MyRI);
				}
			}

			/* Retrieve transition probabilities for Helium.	*/
			/* >>refer He	As	Drake, G.W.F., Atomic, Molecular, and Optical Physics Handbook */
			if( ipHi <= MAX_TP_INDEX && N_(ipHi) <= iso.n_HighestResolved[ipHE_LIKE][ipHELIUM] )
			{
				/*Must not be accessed by collapsed levels!	*/
				ASSERT( ipHi < ( iso.numLevels[ipHE_LIKE][ipHELIUM] - iso.nCollapsed[ipHE_LIKE][ipHELIUM] ) );
				ASSERT( ipLo < ( iso.numLevels[ipHE_LIKE][ipHELIUM] - iso.nCollapsed[ipHE_LIKE][ipHELIUM] ) );
				ASSERT( ipHi > 2 );

				Aul = TransProbs[nelem][ipHi][ipLo];

				putError(nelem,ipHi,ipLo,ipRad,0.0005f);
			}
			
			if( Aul < 0. )
			{
				/* Here are the Lyman transitions.	*/
				if( ipLo == ipHe1s1S )
				{
					ASSERT( (lHi == 1) && (sHi == 0) );

					/* these fits calculated from Drake A's (1996) */
					if ( nLo == 1 )
						Aul = (1.59208e10) / pow(Eff_nupper,3.0);
					ASSERT( Aul > 0.);
					putError(nelem,ipHi,ipLo,ipRad,0.005f);
				}
				
				/* last resort for transitions involving significant defects, 
				 * except that highest lLo are excluded */
				else if( lHi>=2 && lLo>=2 && nHi>nLo )
				{
					/*lint -e790 integral to float */
					Aul = H_Einstein_A(nHi ,lHi , nLo , lLo , nelem);
					/*lint +e790 integral to float */
					ASSERT( Aul > 0.);

					if( lHi + lLo >= 7 )
					{
						putError(nelem,ipHi,ipLo,ipRad,0.001f);
					}
					else
					{
						putError(nelem,ipHi,ipLo,ipRad,0.01f);
					}
				}
				else if( N_(ipHi)>10 && N_(ipLo)<=5 && lHi<=2 && lLo<=2 )
				{
					int paramSet=0;
					double emisOscStr, x, a, b, c;
					double extrapol_Params[2][4][4][3] = {
						/* these are for singlets */
						{	
							{	/* these are P to S */
								{	0.8267396	,	1.4837624	,	-0.4615955	},
								{	1.2738405	,	1.5841806	,	-0.3022984	},
								{	1.6128996	,	1.6842538	,	-0.2393057	},
								{	1.8855491	,	1.7709125	,	-0.2115213	}},
							{	/* these are S to P */
								{	-1.4293664	,	2.3294080	,	-0.0890470	},
								{	-0.3608082	,	2.3337636	,	-0.0712380	},
								{	0.3027974	,	2.3326252	,	-0.0579008	},
								{	0.7841193	,	2.3320138	,	-0.0497094	}},
							{	/* these are D to P */
								{	1.1341403	,	3.1702435	,	-0.2085843	},
								{	1.7915926	,	2.4942946	,	-0.2266493	},
								{	2.1979400	,	2.2785377	,	-0.1518743	},
								{	2.5018229	,	2.1925720	,	-0.1081966	}},
							{	/* these are P to D */
								{	0.0000000	,	0.0000000	,	0.0000000	},
								{	-2.6737396	,	2.9379143	,	-0.3805367	},
								{	-1.4380124	,	2.7756396	,	-0.2754625	},
								{	-0.6630196	,	2.6887253	,	-0.2216493	}},
						},
						/* these are for triplets */
						{	
							{	/* these are P to S */
								{	0.3075287	,	0.9087130	,	-1.0387207	},
								{	0.687069	,	1.1485864	,	-0.6627317	},
								{	0.9776064	,	1.3382024	,	-0.5331906	},
								{	1.2107725	,	1.4943721	,	-0.4779232	}},
							{	/* these are S to P */
								{	-1.3659605	,	2.3262253	,	-0.0306439	},
								{	-0.2899490	,	2.3279391	,	-0.0298695	},
								{	0.3678878	,	2.3266603	,	-0.0240021	},
								{	0.8427457	,	2.3249540	,	-0.0194091	}},
							{	/* these are D to P */
								{	1.3108281	,	2.8446367	,	-0.1649923	},
								{	1.8437692	,	2.2399326	,	-0.2583398	},
								{	2.1820792	,	2.0693762	,	-0.1864091	},
								{	2.4414052	,	2.0168255	,	-0.1426083	}},
							{	/* these are P to D */
								{	0.0000000	,	0.0000000	,	0.0000000	},
								{	-1.9219877	,	2.7689624	,	-0.2536072	},
								{	-0.7818065	,	2.6595150	,	-0.1895313	},
								{	-0.0665624	,	2.5955623	,	-0.1522616	}},
						}
					};

					if( lLo==0 )
					{
						paramSet = 0;
					}
					else if( lLo==1 && lHi==0 )
					{
						paramSet = 1;
					}
					else if( lLo==1 && lHi==2 )
					{
						paramSet = 2;
					}
					else if( lLo==2 )
					{
						paramSet = 3;
						ASSERT( lHi==1 );
					}

					a = extrapol_Params[sHi][paramSet][nLo-2][0];
					b = extrapol_Params[sHi][paramSet][nLo-2][1];
					c = extrapol_Params[sHi][paramSet][nLo-2][2];
					x = log( iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipLo]*RYD_INF/Enerwn );

					emisOscStr = exp(a+b*x+c*x*x)/pow(Eff_nupper,3.)*
						(2.*lLo+1)/(2.*lHi+1);

					Aul = TRANS_PROB_CONST*Enerwn*Enerwn*emisOscStr;

					if ( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P1) || (ipLo == ipHe2p3P2) )
					{
						Aul *= (2.*(ipLo-3)+1.0)/9.0;
					}

					ASSERT( Aul > 0. );
					putError(nelem,ipHi,ipLo,ipRad,0.01f);
				}
				else
				{
					/* Calculate the radial integral from the quantum defects.	*/
					RI2 = scqdri(Eff_nupper,lHi,Eff_nlower,lLo,(double)(ipHELIUM));
					ASSERT( RI2 > 0. );
					/* Convert radial integral to Aul.	*/
					Aul = ritoa(lHi,lLo,ipHELIUM,Enerwn,RI2);
					/* radial integral routine does not recognize fine structure.
					 * Here we split 2^3P.	*/
					if ( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P1) || (ipLo == ipHe2p3P2) )
					{
						Aul *= (2.*(ipLo-3)+1.0)/9.0;
					}

					ASSERT( Aul > 0. );
					putError(nelem,ipHi,ipLo,ipRad,0.03f);
				}
			}
		}

		/* Heavier species	*/
		else
		{
			/* There is probably no need to go to this detail for the ions.
			 * These are remnants of the original version in which Helium was done
			 * in the same loops as the ions.  But of course, since the ions are 
			 * more hydrogenic for any given level, (not to mention far less abundant)
			 * this detail is not necessary, and the algorithm for ions should be simplified.	*/
			
			/* Retrieve transition probabilities for Helike ions.	*/
			/* >>refer He-like	As	Johnson, W.R., Savukov, I.M., Safronova, U.I., & 
			 * >>refercon	Dalgarno, A., 2002, ApJS 141, 543J, originally astro.ph. 0201454 */
			if( ipHi <= MAX_TP_INDEX && N_(ipHi) <= iso.n_HighestResolved[ipHE_LIKE][nelem] )
			{
				/*Must not be accessed by collapsed levels!	*/
				ASSERT( ipHi < ( iso.numLevels[ipHE_LIKE][nelem] - iso.nCollapsed[ipHE_LIKE][nelem] ) );
				ASSERT( ipLo < ( iso.numLevels[ipHE_LIKE][nelem] - iso.nCollapsed[ipHE_LIKE][nelem] ) );
				ASSERT( ipHi > 2 );

				Aul = TransProbs[nelem][ipHi][ipLo];
				putError(nelem,ipHi,ipLo,ipRad,-1);
			}
			
			if( Aul < 0. )
			{
				/* Do same-n transitions. */
				if( nLo == nHi )
				{
					/* These are 2p3Pj to 2s3S fits to (low Z) Porquet & Dubau (2000) & 
					 * (high Z) NIST Atomic Spectra Database.	*/
					if ( ipLo == ipHe2s3S )
					{
						if (ipHi == ipHe2p3P0)
							Aul = 3.31E7 + 1.13E6 * pow((double)nelem+1.,1.76);
						else if (ipHi == ipHe2p3P1)
							Aul = 2.73E7 + 1.31E6 * pow((double)nelem+1.,1.76);
						else if (ipHi == ipHe2p3P2)
							Aul = 3.68E7 + 1.04E7 * exp(((double)nelem+1.)/5.29);
						else
						{
							fprintf(ioQQQ,"WHOA!  Stop in Helike.c");
							cdEXIT(EXIT_FAILURE);
						}
					}
					
					/* These are 2p1P to 2s1S fits to data from TOPbase.	*/
					else if ( ( ipLo == ipHe2s1S ) && ( ipHi == ipHe2p1P) )
					{
						Aul = 5.53E6 * exp( 0.171*(nelem+1.) );
					}
					
					else 
					{
						/* This case should only be entered if n > 2.  Those cases were done above.	*/
						ASSERT( nLo > 2 );

						/* Triplet P to triplet S, delta n = 0	*/
						if ( (lHi == 1) && (sHi == 1) && (lLo == 0) && (sLo == 1))
						{
							Aul = 0.4 * 3.85E8 * pow((double)nelem,1.6)/pow((double)nHi,5.328);
						}
						/* Singlet P to singlet D, delta n = 0	*/
						else if ( (lHi == 1) && (sHi == 0) && (lLo == 2) && (sLo == 0))
						{
							Aul = 1.95E4 * pow((double)nelem,1.6) / pow((double)nHi, 4.269);
						}
						/* Singlet P to singlet S, delta n = 0	*/
						else if ( (lHi == 1) && (sHi == 0) && (lLo == 0) )
						{
							Aul = 6.646E7 * pow((double)nelem,1.5) / pow((double)nHi, 5.077);
						}
						else /*if ( (lHi == 2) && (sHi == 1)  && (lLo == 1) ) */
						{
							Aul = 3.9E6 * pow((double)nelem,1.6) / pow((double)nHi, 4.9);
							if ( (lHi >2) || (lLo > 2) )
								Aul *= (lHi/2.);
							if (lLo > 2)
								Aul *= (1./9.);
						}
					}
					ASSERT( Aul > 0.);
				}

				/* assume transitions involving F and higher orbitals are hydrogenic.	*/
				else if ( ((lHi > 2) || (lLo > 2)) )
				{
					ASSERT( nHi > nLo );
					/*lint -e790 integral to float */
					Aul = H_Einstein_A(nHi ,lHi , nLo , lLo , nelem);
					/*lint +e790 integral to float */
					ASSERT( Aul > 0.);
				}
				
				/* These transitions are of great importance, but the below radial integral 
				 * routine fails to achieve desirable accuracy, so these are fits as produced 
				 * from He A's for nupper through 9.  They are transitions to ground and 
				 * 2, 3, and 4 triplet S.	*/
				else if ( ( ipLo == 0 ) || ( ipLo == ipHe2s1S ) || ( ipLo == ipHe2s3S ) 
					|| ( ipLo == ipHe3s3S ) || ( ipLo == ipHe4s3S ) )
				{
					/* Here are the Lyman transitions.	*/
					if ( ipLo == 0 )
					{
						ASSERT( (lHi == 1) && (sHi == 0) );

						/* In theory, this Z dependence should be Z^4, but values from TOPbase 
						 * suggest 3.9 is a more accurate exponent.	Values from 
						 * >>refer	He-like	As	Johnson, W.R., Savukov, I.M., Safronova, U.I., & 
						 * >>refercon	Dalgarno, A., 2002, ApJS 141, 543J	*/
						/* originally astro.ph. 0201454  */
						Aul = 1.375E10 * pow((double)nelem, 3.9) / pow((double)nHi,3.1);
					}

					/* Here are the Balmer transitions.	*/
					else if ( ipLo == ipHe2s1S )
					{
						ASSERT( (lHi == 1) && (sHi == 0) );

						Aul = 5.0e8 * pow((double)nelem,4.) / pow((double)nHi, 2.889);
					}

					/* Here are transitions down to triplet S	*/
					else
					{
						ASSERT( (lHi == 1) && (sHi == 1) );

						if ( nLo == 2 )
							Aul = 1.5 * 3.405E8 * pow((double)nelem,4.) / pow((double)nHi, 2.883);
						else if ( nLo == 3 )
							Aul = 2.5 * 4.613E7 * pow((double)nelem,4.) / pow((double)nHi, 2.672);
						else 
							Aul = 3.0 * 1.436E7 * pow((double)nelem,4.) / pow((double)nHi, 2.617);
					}

					ASSERT( Aul > 0.);
				}
				
				/* Every other allowed transition is calculated as follows.	*/
				else
				{
					/* Calculate the radial integral from the quantum defects.	*/
					RI2 = scqdri(Eff_nupper,lHi,Eff_nlower,lLo,(double)(nelem));
					/* Convert radial integral to Aul.	*/
					Aul = ritoa(lHi,lLo,nelem,Enerwn,RI2);
					/* radial integral routine does not recognize fine structure.
					 * Here we split 2^3P.	*/
					if ( ( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P1) || (ipLo == ipHe2p3P2) ) && (Aul > iso.SmallA) )
					{
						Aul *= (2.*(ipLo-3)+1.0)/9.0;
					}

				}
				putError(nelem,ipHi,ipLo,ipRad,-1);
			}

			/* TODO - for now just give ions some a 5% error across the board */
			putError(nelem,ipHi,ipLo,ipRad,0.05f);
		}
	}

	/* Now do forbidden transitions from 2-1 ... */
	/* and those given by  
	 * >>refer	He-like	As	Johnson, W.R., Savukov, I.M., Safronova, U.I., & 
	 * >>refercon	Dalgarno, A., 2002, ApJS 141, 543J	*/
	/* originally astro.ph. 0201454  
	 * for heavy elements. These are triplet P to singlet S, 
	 * going either up or down...Triplet S to Singlet P are not included, as they are far weaker.	*/
	else 
	{
		ASSERT( (sHi != sLo) || (abs((int)(lHi - lLo)) != 1) );
		Aul = ForbiddenAuls(ipHi, ipLo, nelem);
		ASSERT( Aul > 0. );
	}
	
	Aul = MAX2( Aul, iso.SmallA );
	ASSERT( Aul >= iso.SmallA );

	/* negative energy for a transition with substantial transition probability
	 * would be major logical error - but is ok for same-n l transitions */
	if( Enerwn < 0. && Aul > iso.SmallA )
	{
		fprintf( ioQQQ," he_1trans hit negative energy, nelem=%li, val was %f \n", nelem ,Enerwn );
	}

	return Aul;
}

void DoFSMixing( long nelem, long ipLoSing, long ipHiSing )
{
	long int nHi, lHi, sHi, nLo, lLo, sLo, ipHiTrip, ipLoTrip;
	float Ass, Att, Ast, Ats;
	float SinHi, SinLo, CosHi, CosLo;
	double HiMixingAngle, LoMixingAngle , error;
	float Kss, Ktt, Kts, Kst, fss, ftt, fssNew, fttNew, ftsNew, fstNew, temp;
	
	nHi = iso.quant_desig[ipHE_LIKE][nelem][ipHiSing].n;
	lHi = iso.quant_desig[ipHE_LIKE][nelem][ipHiSing].l;
	sHi = iso.quant_desig[ipHE_LIKE][nelem][ipHiSing].s;
	nLo = iso.quant_desig[ipHE_LIKE][nelem][ipLoSing].n;
	lLo = iso.quant_desig[ipHE_LIKE][nelem][ipLoSing].l;
	sLo = iso.quant_desig[ipHE_LIKE][nelem][ipLoSing].s;

	if( sHi==1 || sLo ==1 )
		return;

	if( abs(lHi - lLo)!=1 )
		return;

	if( nLo < 2 )
		return;

	if( (lHi<=1) || (lLo<=1) )
		return;

	if( (nHi==nLo) && (lHi==1) && (lLo==2) )
		return; 
	
	if( (nHi > nLo ) && ( lHi != 1 ) && ( lLo != 1) )
		return;
								
	ASSERT( lHi > 0 );
	/*ASSERT( (lHi > 1) && (lLo > 1) );*/
	
	ipHiTrip = QuantumNumbers2Index[nelem][nHi][lHi][1];
	ipLoTrip = QuantumNumbers2Index[nelem][nLo][lLo][1];
	
	if( lHi == 2 )
	{
		HiMixingAngle = 0.01;
	}
	else if( lHi == 3 )
	{
		HiMixingAngle = 0.5;
	}
	else
	{
		HiMixingAngle = PI/4.;
	}
	
	if( lLo == 2 )
	{
		LoMixingAngle = 0.01;
	}
	else if( lLo == 3 )
	{
		LoMixingAngle = 0.5;
	}
	else
	{
		LoMixingAngle = PI/4.;
	}

	/* These would not work correctly if l<=1 were included in this treatment!	*/
	ASSERT( ipHiTrip > ipLoTrip );
	ASSERT( ipHiTrip > ipLoSing );
	ASSERT( ipHiSing > ipLoTrip );
	ASSERT( ipHiSing > ipLoSing );
	
	SinHi = (float)sin( HiMixingAngle );
	SinLo = (float)sin( LoMixingAngle );
	CosHi = (float)cos( HiMixingAngle );
	CosLo = (float)cos( LoMixingAngle );
	
	Kss = EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoSing].EnergyWN;
	Ktt = EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].EnergyWN;
	Kst = EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoTrip].EnergyWN;
	Kts = EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoSing].EnergyWN;

	fss = EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Aul/(float)TRANS_PROB_CONST/POW2(Kss)/
		EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoSing].gLo*
		EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoSing].gHi;

	ftt = EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Aul/(float)TRANS_PROB_CONST/POW2(Ktt)/
		EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].gLo*
		EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].gHi;

	temp = (float)( sqrt(fss/Kss)*CosHi*CosLo + sqrt(ftt/Ktt)*SinHi*SinLo );
	fssNew = Kss*(float)POW2( temp );
	temp = (float)( sqrt(fss/Kss)*SinHi*SinLo + sqrt(ftt/Ktt)*CosHi*CosLo );
	fttNew = Ktt*(float)POW2( temp );
	temp = (float)( sqrt(fss/Kss)*CosHi*SinLo - sqrt(ftt/Ktt)*SinHi*CosLo );
	fstNew = Kst*(float)POW2( temp );
	temp = (float)( sqrt(fss/Kss)*SinHi*CosLo - sqrt(ftt/Ktt)*CosHi*SinLo );
	ftsNew = Kts*(float)POW2( temp );

	Ass = (float)TRANS_PROB_CONST*POW2(Kss)*fssNew*EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoSing].gLo/
		EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoSing].gHi;
	
	Att = (float)TRANS_PROB_CONST*POW2(Ktt)*fttNew*EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].gLo/
		EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].gHi;

	Ast = (float)TRANS_PROB_CONST*POW2(Kst)*fstNew*EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoTrip].gLo/
		EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoTrip].gHi;

	Ats = (float)TRANS_PROB_CONST*POW2(Kts)*ftsNew*EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoSing].gLo/
		EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoSing].gHi;

	error = fabs( ( EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Aul+
		EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Aul /*+
		EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoTrip].Aul+
		EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoSing].Aul */)/
		(Ass+Ast+Ats+Att) - 1.f );

	if( error > 0.001 )
	{
		fprintf( ioQQQ, "FSM error %e LS %li HS %li LT %li HT %li Ratios Ass %e Att %e Ast %e Ats %e\n", error,
			ipLoSing, ipHiSing, ipLoTrip, ipHiTrip,
			Ass/EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Aul,
			Att/EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Aul,
			Ast/EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoTrip].Aul,
			Ats/EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoSing].Aul );
	}

	EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Aul = Ass;
	EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Aul = Att;
	EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoTrip].Aul = Ast;
	EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoSing].Aul = Ats;

	return;
}

/*ritoa converts the square of the radial integral for a transition 
 * (calculated by scqdri) to the transition probability, Aul.	*/
static double ritoa(long li, long lf, long nelem, double k, double RI2)
{  

#define BohrRadiusCM	5.29177249e-9
	
	/*	Variables are as follows:				*/
	/*	lg = larger of li and lf				*/
	/*	fmean = mean oscillator strength		*/
	/*		for a given level.					*/
	/*	mu = reduced mass of optical electron.	*/
	/*	EinsteinA = Einstein emission coef.		*/
	/*	w = angular frequency of transition.	*/
	/*	RI2_cm = square of rad. int. in cm^2.	*/
	long lg;
	double fmean,mu,EinsteinA,w,RI2_cm;
	
	mu = ELECTRON_MASS/(1+ELECTRON_MASS/(dense.AtomicWeight[nelem]*ATOMIC_MASS_UNIT));

	w = 2. * PI * k * SPEEDLIGHT;
	
	RI2_cm = RI2 * BohrRadiusCM * BohrRadiusCM;

	lg = (lf > li) ? lf : li;

	fmean = 2.0*mu*w*lg*RI2_cm/((3.0*H_BAR) * (2.0*li+1.0));

	EinsteinA = TRANS_PROB_CONST*k*k*fmean;

	/* ASSERT( EinsteinA > SMALLFLOAT ); */

	return EinsteinA;
}

#undef DEBUG_LOC 

void HelikeTransProbSetup( void )
{

#	define chLine_LENGTH 1000
	char chLine[chLine_LENGTH] , 
		/* this must be longer than chDataPath, set in path.h */
		chFilename[FILENAME_PATH_LENGTH_2];

	FILE *ioDATA;
	int lgEOL;

	long nelem, ipLo, ipHi, i, i1, i2, i3;
	
	if( (TransProbs = (double ***)MALLOC(sizeof(double **)*(unsigned)LIMELM ) )==NULL )
		BadMalloc();

	for( nelem=ipHELIUM; nelem < LIMELM; ++nelem )
	{

		if( (TransProbs[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(MAX_TP_INDEX+1) ))==NULL )
			BadMalloc();
		
		for( ipLo=ipHe1s1S; ipLo <= MAX_TP_INDEX ;++ipLo )
		{
			if( (TransProbs[nelem][ipLo] = (double*)MALLOC(sizeof(double)*(unsigned)MAX_TP_INDEX ))==NULL )
				BadMalloc();
		}
	}

	/********************************************************************/
	/*************** Read in data from he_transprob.dat	*****************/
	
	/* check on path if file not here and path set */
	/* path was parsed in getset */
	if( lgDataPathSet == TRUE )
	{
		/*path set, so look only there */
		strcpy( chFilename , chDataPath );
		strcat( chFilename , "he_transprob.dat" );
	}
	else
	{
		/* path not set, check local space only */
		strcpy( chFilename , "he_transprob.dat" );
	}

	if( trace.lgTrace )
		fprintf( ioQQQ," HeCreate opening he_transprob.dat:");

	if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
	{
		fprintf( ioQQQ, " HeCreate could not open he_transprob.dat\n" );
		if( lgDataPathSet == TRUE )
			fprintf( ioQQQ, " even tried path\n" );

		if( lgDataPathSet == TRUE )
		{
			fprintf( ioQQQ, " HeCreate could not open he_transprob.dat\n");
			fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
			fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
			fprintf( ioQQQ, " Computation can not continue without he_transprob.dat.\n");
			cdEXIT(EXIT_FAILURE);
		}
	}
	else 
	{
		/* check that magic number is ok */
		if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
		{
			fprintf( ioQQQ, " HeCreate could not read first line of he_transprob.dat.\n");
			puts( "[Stop in HeCreate]" );
			cdEXIT(EXIT_FAILURE);
		}
		i = 1;
		i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
		i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( i1 !=TRANSPROBMAGIC || i2 != N_HE1_TRANS_PROB )
		{
			fprintf( ioQQQ, 
				" HeCreate: the version of he_transprob.dat is not the current version.\n" );
			fprintf( ioQQQ, 
				" HeCreate: I expected to find the number %i %i and got %li %li instead.\n" ,
				TRANSPROBMAGIC, N_HE1_TRANS_PROB, i1, i2 );
			fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
			puts( "[Stop in HeCreate]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* Initialize TransProbs[nelem][ipLo][ipHi] before filling it in.	*/
		for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
		{	
			for( ipHi=0; ipHi <= MAX_TP_INDEX; ipHi++ )
			{
				for( ipLo=0; ipLo < MAX_TP_INDEX; ipLo++ )
				{
					TransProbs[nelem][ipHi][ipLo] = -1.;
				}
			}
		}

		for( ipLo=1; ipLo <= N_HE1_TRANS_PROB; ipLo++ )
		{
			char *chTemp;
			
			/* get next line image */
			if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
				BadRead();
			
			while( chLine[0]=='#' )
			{
				if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
					BadRead();
			}

			i3 = 1;
			i1 = (long)FFmtRead(chLine,&i3,INPUT_LINE_LENGTH,&lgEOL);
			i2 = (long)FFmtRead(chLine,&i3,INPUT_LINE_LENGTH,&lgEOL);
			/* check that these numbers are correct */
			if( i1<0 || i2<=i1 )
			{
				fprintf( ioQQQ, " HeCreate detected insanity in he_transprob.dat.\n");
				puts( "[Stop in HeCreate]" );
				cdEXIT(EXIT_FAILURE);
			}
				
			chTemp = chLine;
			
			/* skip over 2 tabs to start of data */
			for( i=0; i<1; ++i )
			{
				if( (chTemp = strchr( chTemp, '\t' )) == NULL )
				{
					fprintf( ioQQQ, " HeCreate could not init he_transprob\n" );
					puts( "[Stop in HeCreate]" );
					cdEXIT(EXIT_FAILURE);
				}
				++chTemp;
			}

			/* now read in the data */
			for( nelem = ipHELIUM; nelem < LIMELM; nelem++ )
			{
				float a;
				
				if( (chTemp = strchr( chTemp, '\t' )) == NULL )
				{
					fprintf( ioQQQ, " HeCreate could not scan he_transprob\n" );
					puts( "[Stop in HeCreate]" );
					cdEXIT(EXIT_FAILURE);
				}
				++chTemp;

				sscanf( chTemp , "%e" , &a );
				TransProbs[nelem][i2][i1] = a;
				/*TransProbs[nelem][i2][i1] = FFmtRead(chLine,&i3,chLine_LENGTH,&lgEOL);*/
				
				if( lgEOL )
				{
					fprintf( ioQQQ, " HeCreate detected insanity in he_transprob.dat.\n");
					puts( "[Stop in HeCreate]" );
					cdEXIT(EXIT_FAILURE);
				}
			}
		}

		/* check that ending magic number is ok */
		if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
		{
			fprintf( ioQQQ, " HeCreate could not read last line of he_transprob.dat.\n");
			puts( "[Stop in HeCreate]" );
			cdEXIT(EXIT_FAILURE);
		}
		i = 1;
		i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
		i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( i1 !=TRANSPROBMAGIC || i2 != N_HE1_TRANS_PROB )
		{
			fprintf( ioQQQ, 
				" HeCreate: the version of he_transprob.dat is not the current version.\n" );
			fprintf( ioQQQ, 
				" HeCreate: I expected to find the number %i %i and got %li %li instead.\n" ,
				TRANSPROBMAGIC, N_HE1_TRANS_PROB, i1, i2 );
			fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
			puts( "[Stop in HeCreate]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* close the data file */
		fclose( ioDATA );
	}

	return;
}

static double CalculateRadialIntegral( long nelem, long ipLo, long ipHi )
{
	double RadInt = 0.;

	RI_ipHi = ipHi;
	RI_ipLo = ipLo;

	/* Z is needed also, but cannot be passed to integrand directly.	*/
	globalZ = nelem;
	
	if( ( nelem == ipHELIUM) )
	{
		/* Evaluate radial integral	out to 10 Bohr radii.	*/
		RadInt = qg32( 0., 0.2, RadialIntegrand );
		RadInt += qg32( 0.2, 0.4, RadialIntegrand );
		RadInt += qg32( 0.4, 0.6, RadialIntegrand );
		RadInt += qg32( 0.8, 1.0, RadialIntegrand );
		RadInt += qg32( 1.0, 2.0, RadialIntegrand );
		RadInt += qg32( 2.0, 4.0, RadialIntegrand );
		RadInt += qg32( 4.0, 10., RadialIntegrand );
	}
	else
		fprintf(ioQQQ, "No can do in CalculateRadialIntegral, nelem %li, iplo %li, ipHi %li\n",
			nelem, ipLo, ipHi);
	
	/* mulitply by normalization constants */
	RadInt *= 1.;

	return RadInt;
}

static double RadialIntegrand( double r12 )
{
	double result;

	RI_ipLev = RI_ipHi;

	result = r12 * WaveFunction( r12 ); 
	
	RI_ipLev = RI_ipLo;

	result *= WaveFunction( r12 );

	return result;
}

static double WaveFunction( double r12 )
{
	double value,PsiSlow2,alphaOfR2, PsiFastOneTwo;
	long Z = globalZ;
	/* use classical radius of electron for R1.	*/
	/* Should it be multiplied by Z^2?	*/
	double R1 = POW2((double)Z) * 5.3E-5;
	double R2 = r12 + R1;
	double a=0.,b=0.,mu=0.,nu=0.,Beta=0.;
	long l = iso.quant_desig[ipHE_LIKE][Z][RI_ipLev].l;
	/* long s = iso.quant_desig[ipHE_LIKE][Z][RI_ipLev].s;	*/


	/*  This wavefunction from Koga and Matsuhashi 89.
	 *  It is a product of two wavefunctions...one for the 1s ground electron,
	 *  and the other for the excited state.	*/

	/* Excited function	*/
	if( l == 0 )
	{
		Beta = 0.411;
		PsiSlow2 = sqrt( POW3(Beta) / PI ) * exp( -1. * Beta * R2 );
		/* parameters for alphaOfR2, for S orbitals.	*/
		a = -0.1;
		b = 0.14;
		mu = 1.34;
		nu = 2.86;
	}
	else if( l == 1 )
	{
		Beta = 0.5;
		/* The leading number is sqrt( POW3(2/sqrt(3)) / PI )	*/
		PsiSlow2 = 1.1547 * pow( Beta, 2.5 ) * R2 * exp( -1. * Beta * R2 );
		/* parameters for alphaOfR2, for P orbitals.	*/
		a = -1;
		b = 0.2;
		mu = 0.1;
		nu = 4.0;
	}
	else
	{
		/* for now only do S and P	*/
		PsiSlow2 = 0.;
	}

	alphaOfR2 = 2 - (1. + a * pow( R2, mu ) )/ ( 1. + b * pow( R2, nu ) );
	
	/* Ground function	*/
	PsiFastOneTwo = sqrt( POW3(alphaOfR2) / PI ) * exp( -1. * alphaOfR2 * R1 );

	value = PsiFastOneTwo * PsiSlow2;

	return value;
}
/*lint +e662 creation of  out of bound pointer */

#undef PARALLEL
/*lint +e662 creation of  out of bound pointer */
/*lint +e661 creation of  out of bound pointer */









