/* 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 */
/*ContCreatePointers set up pointers for lines and continua called by cloudy after input read in 
 * and continuum mesh has been set */
/*fiddle adjust energy bounds of certain cells so that match ionization edges exactly */
/*ipShells assign continuum energy pointers to shells for all atoms,
 * called by ContCreatePointers */
/*LimitSh sets upper energy limit to subshell integrations */
#include "cddefines.h"
#include "physconst.h"
#include "lines_service.h"
#include "iso.h"
#include "secondaries.h"
#include "taulines.h"
#include "elementnames.h"
#include "ionbal.h"
#include "helike.h" 
#include "rt.h"
#include "opacity.h"
#include "yield.h"
#include "dense.h"
#include "he.h"
#include "fe.h"
#include "rfield.h"
#include "oxy.h"
#include "atomfeii.h"
#include "atoms.h"
#include "trace.h"
#include "kshllenr.h"
#include "hmi.h"
#include "heavy.h"
#include "predcont.h"
#include "atmdat.h"
#include "ipoint.h"
#include "h2.h"
#include "continuum.h"

/*LimitSh sets upper energy limit to subshell integrations */
static long LimitSh(long int ion, 
  long int nshell, 
  long int nelem);

static void ipShells(
	/* nelem is the atomic number on the C scale, Li is 2 */
	long int nelem);

#define TwoS	(1+ipISO)

/*fiddle adjust energy bounds of certain cells so that match ionization edges exactly */
static void fiddle(long int ipnt, 
  double exact);

void ContCreatePointers(void)
{
	char chLab[5];
	long int 
	  i, 
	  ion, 
	  ipHi, 
	  ipLo, 
	  ipISO,
	  iWL_Ang,
	  j, 
	  nelem,
	  nshells;
	double energy,
		xnew;
	/* counter to say whether pointers have ever been evaluated */
	static int nCalled = 0;

#	ifdef DEBUG_FUN
	fputs( "<+>ContCreatePointers()\n", debug_fp );
#	endif

	/* create the hydrogen atom for this coreload, routine creates space then zeros it out
	 * on first call, on second and later calls it only zeros things out */
	iso_create();

	/* create internal static variables needed to do the he-like series */
	HeCreate();

	/* create internal static variables needed to do the H2 molecule */
	H2_Create();

	/* nCalled is local static variable defined 0 when defined. 
	 * it is incremented below, so that space only allocated one time per coreload. */
	if( nCalled > 0 )
	{
		if( trace.lgTrace )
		{
			fprintf( ioQQQ, " ContCreatePointers called, not evaluating.\n" );
		}
		
#		ifdef DEBUG_FUN
		fputs( " <->ContCreatePointers()\n", debug_fp );
#		endif
		return;
	}
	else
	{
		if( trace.lgTrace )
		{
			fprintf( ioQQQ, " ContCreatePointers called first time.\n" );
		}
		++nCalled;
	}

	for( i=0; i < rfield.nupper; i++ )
	{
		/* this is array of labels for lines and continua, set to blanks at first */
		strcpy( rfield.chContLabel[i], "    ");
		strcpy( rfield.chLineLabel[i], "    ");
	}

	/* we will generate a set of array indices to ionization edges for
	 * the first thirty elements.  First set all array indices to
	 * totally bogus values so we will crash if misused */
	for( nelem=0; nelem<LIMELM; ++nelem )
	{
		if( dense.lgElmtOn[nelem] )
		{
			for( ion=0; ion<LIMELM; ++ion )
			{
				for( nshells=0; nshells<7; ++nshells )
				{
					for( j=0; j<3; ++j )
					{
						opac.ipElement[nelem][ion][nshells][j] = INT_MIN;
					}
				}
			}
		}
	}

	/* pointer to excited state of O+2 */
	opac.ipo3exc[0] = ipContEnergy(3.855,"O3ex");

	/* main hydrogenic arrays - THIS OCCURS TWICE!! H and He here, then the
	 * remaining hydrogenic species near the bottom.  This is so that H and HeII get
	 * their labels stuffed into the arrays, and the rest of the hydrogenic series 
	 * get whatever is left over after the level 1 lines.
	 * to find second block, search on "ipZ=2" */
	/* NB note that no test for H or He being on exists here - we will always
	 * define the H and He arrays even when he is off, since we need to
	 * know where the he edges are for the bookkeeping that occurs in continuum
	 * binning routines */
	/* this loop is over H, He-like only */
	for( ipISO=ipH_LIKE; ipISO<=ipHE_LIKE; ++ipISO )
	{
		/* this will be over HI, HeII, then HeI only */
		for( nelem=ipISO; nelem < 2; nelem++ )
		{
			/* generate label for this ion */
			sprintf( chLab, "%2s%2ld",elementnames.chElementSym[nelem], nelem+1-ipISO );

			/* array index for continuum edges for ground */
			iso.ipIsoLevNIonCon[ipISO][nelem][0] = ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][0],chLab);
			for( ipHi=1; ipHi < iso.numLevels[ipISO][nelem]; ipHi++ )
			{
				/* array index for continuum edges for excited levels */
				iso.ipIsoLevNIonCon[ipISO][nelem][ipHi] = ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][ipHi],chLab);

				/* define all line array indices */
				for( ipLo=0; ipLo < ipHi; ipLo++ )
				{
					/* some lines have negative or zero energy */
					/* >>chng 03 apr 22, this check was if less than or equal to zero,
					 * changed to lowest energy point so that ultra low energy transitions are
					 * not considered.	*/
					if( EmisLines[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD < continuum.filbnd[0] )
					{
						EmisLines[ipISO][nelem][ipHi][ipLo].ipCont = -1;
						EmisLines[ipISO][nelem][ipHi][ipLo].ipFine = -1;
					}
					else
					{
						/* some energies are negative for inverted levels */
						EmisLines[ipISO][nelem][ipHi][ipLo].ipCont = 
							(int)ipLineEnergy(EmisLines[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD , chLab ,
							iso.ipIsoLevNIonCon[ipISO][nelem][ipLo]);
						EmisLines[ipISO][nelem][ipHi][ipLo].ipFine = 
							ipFineCont(EmisLines[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD );
						/* check that energy scales are the same, to within energy resolution of coa */
						ASSERT( EmisLines[ipISO][nelem][ipHi][ipLo].ipFine < 0 ||
							fabs(rfield.anu[EmisLines[ipISO][nelem][ipHi][ipLo].ipCont-1] -
							rfield.fine_anu[EmisLines[ipISO][nelem][ipHi][ipLo].ipFine-1]) /
							rfield.anu[EmisLines[ipISO][nelem][ipHi][ipLo].ipCont-1] < 
							/* can only assert more than width of coarse cell since close to ionization limit,
							 * course index is kept below next ionization edge */
							1.5*rfield.widflx[EmisLines[ipISO][nelem][ipHi][ipLo].ipCont-1]/rfield.anu[EmisLines[ipISO][nelem][ipHi][ipLo].ipCont-1]);
					}
				}
			}/* ipHi loop */
		}/* nelem loop */
	}/* ipISO */

	/* need to increase the cell for the HeII balmer continuum by one, so that
	 * hydrogen lyman continuum will pick it up. */
	nelem = 1;
	/* copy label over to next higher cell */
	if( strcmp( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]-1] , "He 2" ) == 0)
	{
		strcpy( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]], 
				 rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]-1] );
		/* set previous spot to blank */
		strcpy( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]-1] , "    ");
	}
	else if( strcmp( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]-1] , "H  1" ) == 0)
	{
		/* set previous spot to blank */
		strcpy( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]] , "He 2");
	}
	else
	{
		fprintf(ioQQQ," insanity heii pointer fix contcreatepointers\n");
	}
	/* finally increment the two HeII pointers so that they are above the Lyman continuum */
	++iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s];
	++iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2p];

	/* array of either 1 or 0 for energy of elec to be counted
	 * as able to produce secondary ionizations
	 * below 100eV no secondary ionization */
	secondaries.ipSecIon = ipoint(7.353);

	/* this is highest energy where k-shell opacities are counted
	 * can be adjusted with "set kshell" command */
	KshllEnr.KshellLimit = ipoint(KshllEnr.EnergyKshell);

	/* pointers for molecules
	 * H2+ dissociation energy 2.647 eV but cs small below 0.638 Ryd */
	opac.ih2pnt[0] = ipContEnergy(0.502,"H2+ ");
	opac.ih2pnt[1] = ipoint(1.03);

	/* pointers for most prominent PAH features
	 * energies given to ipContEnergy are only to put lave in the right place
	 * wavelenghts are rough observed values of blends
	 * 7.6 microns */
	i = ipContEnergy(0.0117, "PAH " );

	/* feature near 6.2 microns */
	i = ipContEnergy(0.0147, "PAH " );

	/* 3.3 microns */
	i = ipContEnergy(0.028, "PAH " );

	/* 11.2 microns */
	i = ipContEnergy(0.0080, "PAH " );

	/* 12.3 microns */
	i = ipContEnergy(0.0077, "PAH " );

	/* 13.5 microns */
	i = ipContEnergy(0.0069, "PAH " );


	/* fix pointers for hydrogen and helium */

	/* pointer to Bowen 374A resonance line */
	he.ip374 = ipLineEnergy(1.92,"He 2",0);
	he.ip660 = ipLineEnergy(1.38,"He 2",0);

	/* set up series of continuum pointers to be used in routine lines to
	 * enter continnum points into line array */
	for( i=0; i < NPREDCONT; i++ )
	{
		/* EnrPredCont contains array of energy points, as set in zerologic.c */
		ipPredCont[i] = ipoint(EnrPredCont[i]) - 1;
	}

	/* pointer to energy defining effect x-ray column */
	rt.ipxry = ipoint(73.5);

	/* pointer to Hminus edge at 0.754eV */
	hmi.iphmin = ipContEnergy(0.05544,"H-  ");

	hmi.iheh1 = ipoint(1.6);
	hmi.iheh2 = ipoint(2.3);

	/* pointer to carbon k-shell ionization */
	opac.ipCKshell = ipoint(20.6);

	/* pointer to threshold for pair production */
	opac.ippr = ipoint(7.51155e4) + 1;

	/* pointer to x-ray - gamma ray bound; 100 kev */
	rfield.ipEnerGammaRay = ipoint(rfield.EnerGammaRay);

	Fe2_ovr_DataInit();
	for( i=0; i < NFEII; i++ )
	{
		fe2ovr_la.ipfe2[i] = ipoint(fe2ovr_la.fe2enr[i]);
	}

	/* make low energy edges of some cells exact,
	 * index is on c scale
	 * 0.99946 is correct energy of hydrogen in inf mass ryd */
	fiddle(iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1,0.99946);
	fiddle(iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2p]-1,0.24986);
	/* confirm that labels are in correct location */
	ASSERT( strcmp( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1], "H  1" ) ==0 );
	ASSERT( strcmp( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2p]-1], "H  1" ) ==0 );

	fiddle(iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s]-1,4.00000);
	ASSERT( strcmp( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s]-1], "He 2" ) ==0 );

	/* pointer to excited state of O+2
	 * ipo3exc(1) = ipContEnergy(3.855 ,'O3ex' ) */
	fiddle(opac.ipo3exc[0]-1,3.855);

	/* now fix widflx array so that it is correct */
	for( i=1; i<rfield.nupper-1; ++i )
	{
		rfield.widflx[i] = ((rfield.anu[i+1] - rfield.anu[i]) + (rfield.anu[i] - rfield.anu[i-1]))/2.f;
	}

	/* these are indices for centers of B and V filters,
	 * taken from table on page 202 of Allen, AQ, 3rd ed */
	/* the B filter array offset */
	rfield.ipB_filter = ipoint( RYDLAM / WL_B_FILT );
	/* the V filter array offset */
	rfield.ipV_filter = ipoint( RYDLAM / WL_V_FILT );

	/* these are the lower and upper bounds for the G0 radiation field
	 * used by Tielens & Hollenbach in their PDR work */
	rfield.ipG0_TH85_lo =  ipoint(  6.0 / EVRYD );
	rfield.ipG0_TH85_hi =  ipoint( 13.6 / EVRYD );

	/* this is the limits for Draine & Bertoldi Habing field */
	rfield.ipG0_DB96_lo =  ipoint(  RYDLAM / 1110. );
	rfield.ipG0_DB96_hi =  ipoint( RYDLAM / 912. );

	/* this index is to 1000A to obtain the extinction at 1000A */
	rfield.ip1000A = ipoint( RYDLAM / 1000. );

	/* now save current form of array */
	for( i=0; i < rfield.nupper; i++ )
	{
		rfield.AnuOrg[i] = rfield.anu[i];
		rfield.anusqr[i] = (float)sqrt(rfield.AnuOrg[i]);
	}

	/* following order of elements is in roughly decreasing abundance
	 * when ipShells gets a cell for the valence IP it does so through
	 * ipContEnergy, which makes sure that no two ions get the same cell
	 * earliest elements have most precise ip mapping */

	/* set up shell pointers for hydrogen */
	nelem = ipHYDROGEN;
	ion = 0;

	/* the number of shells */
	Heavy.nsShells[nelem][0] = 1;

	/*pointer to ionization threshold in energy array*/
	Heavy.ipHeavy[nelem][ion] = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s];
	opac.ipElement[nelem][ion][0][0] = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s];

	/* upper limit to energy integration */
	opac.ipElement[nelem][ion][0][1] = rfield.nupper;

	if( dense.lgElmtOn[ipHELIUM] )
	{
		/* set up shell pointers for helium */
		nelem = ipHELIUM;
		ion = 0;

		/* the number of shells */
		Heavy.nsShells[nelem][0] = 1;

		/*pointer to ionization threshold in energy array*/
		Heavy.ipHeavy[nelem][ion] = iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0];
		opac.ipElement[nelem][ion][0][0] = iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0];

		/* upper limit to energy integration */
		opac.ipElement[nelem][ion][0][1] = rfield.nupper;

		/* (hydrogenic) helium ion */
		ion = 1;
		/* the number of shells */
		Heavy.nsShells[nelem][1] = 1;

		/*pointer to ionization threshold in energy array*/
		Heavy.ipHeavy[nelem][ion] = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s];
		opac.ipElement[nelem][ion][0][0] = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s];

		/* upper limit to energy integration */
		opac.ipElement[nelem][ion][0][1] = rfield.nupper;
	}

	/* check that ionization potential of neutral carbon valence shell is
	 * positive */
	ASSERT( PH1COM.PH1[2][5][5][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
	 * iso.ipIsoLevNIonCon[ipH_LIKE] rather than ipoint in getting array index within continuum array */
	for( i=2; i<LIMELM; ++i )
	{
		if( dense.lgElmtOn[i])
		{
			/* i is the atomic number on the c scale, 2 for Li */
			ipShells(i);
		}
	}

	/* most of these are set in ipshells, but not h-like or he-like, so do these here */
	Heavy.Valence_IP_Ryd[ipHYDROGEN][0] = PH1COM.PH1[0][0][ipHYDROGEN][0]/EVRYD* 0.9998787 ;
	Heavy.Valence_IP_Ryd[ipHELIUM][0] = PH1COM.PH1[0][1][ipHELIUM][0]/EVRYD* 0.9998787 ;
	Heavy.Valence_IP_Ryd[ipHELIUM][1] = PH1COM.PH1[0][0][ipHELIUM][0]/EVRYD* 0.9998787 ;
	for( nelem=2; nelem<LIMELM; ++nelem )
	{
		Heavy.Valence_IP_Ryd[nelem][nelem-1] = PH1COM.PH1[0][1][nelem][0]/EVRYD* 0.9998787 ;
		Heavy.Valence_IP_Ryd[nelem][nelem] = PH1COM.PH1[0][0][nelem][0]/EVRYD* 0.9998787 ;
		if( dense.lgElmtOn[nelem])
		{
			/* now confirm that all are properly set */
			for( j=0; j<=nelem; ++j )
			{
				ASSERT( Heavy.Valence_IP_Ryd[nelem][j] > 0.05 );
			}
			for( j=0; j<nelem; ++j )
			{
				ASSERT( Heavy.Valence_IP_Ryd[nelem][j] < Heavy.Valence_IP_Ryd[nelem][j+1]);
			}
		}
	}

	/* array indices to bound compton electron recoil ionization of all elements */
	for( nelem=0; nelem<LIMELM; ++nelem )
	{
		if( dense.lgElmtOn[nelem])
		{
			for( ion=0; ion<nelem+1; ++ion )
			{
				/* this is the threshold energy to compton ionize valence shell electrons */
				energy = sqrt( Heavy.Valence_IP_Ryd[nelem][ion] * EN1RYD * ELECTRON_MASS * SPEEDLIGHT * SPEEDLIGHT ) / EN1RYD;
				/* the array index for this energy */
				ionbal.ipCompRecoil[nelem][ion] = ipoint( energy );
			}
		}
	}

	/* oxygen
	 * pointers for excited states
	 * IO3 is pointer to O++ exc state, is set above */
	oxy.i2d = ipoint(1.242);
	oxy.i2p = ipoint(1.367);
	opac.ipo1exc[0] = ipContEnergy(0.856,"O1ex");
	opac.ipo1exc[1] = ipoint(2.0);

	/* atomic oxygen valence shell 
	Heavy.nsShells[ipOXYGEN][0]*/
	/* this is the valence shell index */
	/* finish off excited state photoionization
	 * do not use ipContEnergy since only upper limit */
	opac.ipo3exc[1] = ipoint(5.0);

	/* upper level of 4363 */
	opac.ipo3exc3[0] = ipContEnergy(3.646,"O3ex");
	opac.ipo3exc3[1] = ipoint(5.0);

	/* following are various pointers for OI - Lybeta pump problem
	 * these are delta energies for boltzmann factors
	 * TODO	2	this is redundant with contents of oxygen line arrays
	 * use them instead
	 * when removing this, make sure all line intensity predictions
	 * also go into oi line arrays */
	atoms.ipoiex[0] = ipoint(0.7005);
	atoms.ipoiex[1] = ipoint(0.10791);
	atoms.ipoiex[2] = ipoint(0.06925);
	atoms.ipoiex[3] = ipoint(0.01151);
	atoms.ipoiex[4] = ipoint(0.01999);

	/* >>chng 97 jan 27, move nitrogen after oxygen so that O gets the
	 * most accurate pointers
	 * Nitrogen
	 * in1(1) is thresh for photo from excited state */
	opac.in1[0] = ipContEnergy(0.893,"N1ex");

	/* upper limit */
	opac.in1[1] = ipoint(2.);

	if( (trace.lgTrace && trace.lgConBug) || (trace.lgTrace && trace.lgPointBug) )
	{
		fprintf( ioQQQ, "   ContCreatePointers:%ld energy cells used. N(1R):%4ld N(1.8):%4ld  N(4Ryd):%4ld N(O3)%4ld  N(x-ray):%5ld N(rcoil)%5ld\n", 
		  rfield.nupper, 
		  iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s], 
		  iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][ipH1s], 
		  iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s], 
		  opac.ipo3exc[0], 
		  opac.ipCKshell, 
		  ionbal.ipCompRecoil[ipHYDROGEN][0] );


		fprintf( ioQQQ, "   ContCreatePointers: ipEnerGammaRay: %5ld IPPRpari produc%5ld\n", 
		  rfield.ipEnerGammaRay, opac.ippr );

		fprintf( ioQQQ, "   ContCreatePointers: H pointers;" );
		for( i=0; i <= 6; i++ )
		{
			fprintf( ioQQQ, "%4ld%4ld", i, iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][i] );
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "   ContCreatePointers: Oxy pnters;" );

		for( i=1; i <= 8; i++ )
		{
			fprintf( ioQQQ, "%4ld%4ld", i, Heavy.ipHeavy[ipOXYGEN][i-1] );
		}
		fprintf( ioQQQ, "\n" );
	}

	/* Magnesium
	 * following is energy for phot of MG+ from exc state producing 2798 */
	opac.ipmgex = ipoint(0.779);

	/* pointers to ir triplet, kh
	 * icaxyz = ipoint( 0.11 )
	 * icakh = ipoint( 0.235 )
	 * lower, upper edges of Ca+ excited term photoionizaiton
	 * ica2ex(1) = ipoint(0.72) */
	opac.ica2ex[0] = ipContEnergy(0.72,"Ca2x");
	opac.ica2ex[1] = ipoint(1.);

	/* set all pointers to ionizaiton edges with labels */
	/* set up factors and pointers for Fe continuum fluorescence */
	fe.ipfe10 = ipoint(2.605);

	/* following is WL(CM)**2/(8PI) * conv fac for RYD to NU *A21 */
	fe.pfe10 = (float)(2.00e-18/rfield.widflx[fe.ipfe10-1]);

	/* this is 353 pump, f=0.032 */
	fe.pfe11a = (float)(4.54e-19/rfield.widflx[fe.ipfe10-1]);

	/* this is 341.1 f=0.012 */
	fe.pfe11b = (float)(2.75e-19/rfield.widflx[fe.ipfe10-1]);
	fe.pfe14 = (float)(1.15e-18/rfield.widflx[fe.ipfe10-1]);

	/* set up energy pointers for line optical depth arrays
	 * this also increments flux, sets other parameters for lines */

	/* level 1 heavy elements line array */
	for( i=1; i <= nLevel1; i++ )
	{
		/* put null terminated line label into chLab using line array*/
		chIonLbl(chLab, &TauLines[i]);

		TauLines[i].ipCont = 
			(int)ipLineEnergy(TauLines[i].EnergyWN * WAVNRYD, chLab ,0);
		TauLines[i].ipFine = 
			ipFineCont(TauLines[i].EnergyWN * WAVNRYD );
		/* for debugging pointer - index on f not c scale, 
		 * this will find all lines that entered a given cell 
		   if( TauLines[i].ipCont==603 )
			fprintf(ioQQQ,"( level1 %s\n", chLab);*/

		if( TauLines[i].gf > 0. )
		{
			/* the gf value was entered
			 * derive the A, call to function is gf, wl (A), g_up */
			TauLines[i].Aul = 
				(float)(eina(TauLines[i].gf,
			  TauLines[i].EnergyWN,TauLines[i].gHi));
		}
		else if( TauLines[i].Aul > 0. )
		{
			/* the Einstion A value was entered
			 * derive the gf, call to function is A, wl (A), g_up */
			TauLines[i].gf = 
				(float)(GetGF(TauLines[i].Aul,
			  TauLines[i].EnergyWN,TauLines[i].gHi));
		}
		else
		{
			fprintf( ioQQQ, " level 1 line does not have valid gf or A\n" );
			fprintf( ioQQQ, " This is ContCreatePointers\n" );
			puts( "[Stop in ContCreatePointers]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* used to get damping constant */
		/* >>chng 02 aug 13, in terms of energy in wavenumbers */
		TauLines[i].dampXvel = 
			/*(float)(TauLines[i].Aul*
		  TauLines[i].WLAng*1e-8/PI4);*/
			(float)(TauLines[i].Aul/
		  TauLines[i].EnergyWN/PI4);

		/* derive the abs coef, call to function is gf, wl (A), g_low */
		TauLines[i].opacity = 
			(float)(abscf(TauLines[i].gf,
		  TauLines[i].EnergyWN,TauLines[i].gLo));
		/*fprintf(ioQQQ,"TauLinesss\t%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
			i,TauLines[i].opacity , TauLines[i].gf , TauLines[i].Aul ,TauLines[i].EnergyWN,TauLines[i].gLo);*/

		/* excitation energy of transition in degrees kelvin */
		TauLines[i].EnergyK = 
			(float)(T1CM)*TauLines[i].EnergyWN;

		/* energy of photon in ergs */
		TauLines[i].EnergyErg = 
			(float)(ERG1CM)*TauLines[i].EnergyWN;

		/*line wavelength must be gt 0 */
		ASSERT( TauLines[i].WLAng > 0 );
	}

	/* set the ipCont struc element for the H2 molecule */
	H2_ContPoint();

	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	{
		/* do remaining part of the he iso sequence */
		for( nelem=2; nelem < LIMELM; nelem++ )
		{
			if( dense.lgElmtOn[nelem])
			{
				/* generate label for this ion */
				sprintf( chLab, "%2s%2ld",elementnames.chElementSym[nelem], nelem+1-ipISO );
				/* Lya itself is the only transition below n=3 - we explicitly do not
				* want to generate pointers for 2s-1s or 2p-2s */
				/* array index for continuum edges for levels in He-like ions  */
				iso.ipIsoLevNIonCon[ipISO][nelem][0] = 
					ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][0],chLab);

				for( ipHi=1; ipHi < iso.numLevels[ipISO][nelem] ; ipHi++ )
				{
					/* array index for continuum edges for levels in He-like ions  */
					iso.ipIsoLevNIonCon[ipISO][nelem][ipHi] = ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][ipHi],chLab);

					/* define all he-like line pointers */
					for( ipLo=0; ipLo < ipHi; ipLo++ )
					{
						/* some lines have negative or zero energy */
						/* >>chng 03 apr 22, this check was if less than or equal to zero,
						 * changed to lowest energy point so that ultra low energy transitions are
						 * not considered.	*/
						if( EmisLines[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD < continuum.filbnd[0] )
						{
							/* some energies are negative for inverted levels */
							EmisLines[ipISO][nelem][ipHi][ipLo].ipCont = -1;
							EmisLines[ipISO][nelem][ipHi][ipLo].ipFine = -1;
						}
						else
						{
							EmisLines[ipISO][nelem][ipHi][ipLo].ipCont = 
								(int)ipLineEnergy(EmisLines[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD , chLab ,
								iso.ipIsoLevNIonCon[ipISO][nelem][ipLo]);
							EmisLines[ipISO][nelem][ipHi][ipLo].ipFine = 
								ipFineCont(EmisLines[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD );
						}
					}
				}
				iso.ipIsoLevNIonCon[ipISO][nelem][0] = ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][0],chLab);
			}
		}
	}
	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	{
		/* this will be over HI, HeII, then HeI only */
		for( nelem=ipISO; nelem < LIMELM; nelem++ )
		{
			if( dense.lgElmtOn[nelem])
			{
				ipLo = 0;
				/* these are the extra lyman lines */
				for( ipHi=2; ipHi < iso.nLyman[ipISO]; ipHi++ )
				{
					/* some energies are negative for inverted levels */
					iso.ExtraLymanLines[ipISO][nelem][ipHi].ipCont = 
						(int)ipLineEnergy(iso.ExtraLymanLines[ipISO][nelem][ipHi].EnergyWN * WAVNRYD , chLab ,
						iso.ipIsoLevNIonCon[ipISO][nelem][ipLo]);

					iso.ExtraLymanLines[ipISO][nelem][ipHi].ipFine = 
						ipFineCont(iso.ExtraLymanLines[ipISO][nelem][ipHi].EnergyWN * WAVNRYD );
				}
			}
		}
	}

	/* some lines do not exist, the flag indicating this is ipCont == -1 */
	/* for H-like sequence, these are 2p-2s (energies degenerate) and 2s-1s, two photon */
	ipISO = ipHYDROGEN;
	/* do remaining part of the he iso sequence */
	for( nelem=ipISO; nelem < LIMELM; nelem++ )
	{
		if( dense.lgElmtOn[nelem])
		{
			/* for H-like sequence want 2p-2s (energies degenerate) and 2s-1s, two photon */
			EmisLines[ipISO][nelem][ipH2p][ipH2s].ipCont = -1;
			EmisLines[ipISO][nelem][ipH2p][ipH2s].ipFine = -1;
			EmisLines[ipISO][nelem][ipH2s][ipH1s].ipCont = -1;
			EmisLines[ipISO][nelem][ipH2s][ipH1s].ipFine = -1;
		}
	}

	/* for He-like sequence the majority of the transitions are bogus - A set to special value in this case */
	ipISO = ipHELIUM;
	/* do remaining part of the he iso sequence */
	for( nelem=ipISO; nelem < LIMELM; nelem++ )
	{
		if( dense.lgElmtOn[nelem])
		{
			/* this is the two photon transition in the singlets */
			EmisLines[ipISO][nelem][ipHe2s1S][ipHe1s1S].ipCont = -1;
			EmisLines[ipISO][nelem][ipHe2s1S][ipHe1s1S].ipFine = -1;

			for( ipHi=1; ipHi < iso.numLevels[ipISO][nelem] ; ipHi++ )
			{
				for( ipLo=0; ipLo < ipHi; ipLo++ )
				{
					if( fabs(EmisLines[ipISO][nelem][ipHi][ipLo].Aul - iso.SmallA) < SMALLFLOAT )
					{
						/* iso.SmallA is value assigned to bogus transitions */
						EmisLines[ipISO][nelem][ipHi][ipLo].ipCont = -1;
						EmisLines[ipISO][nelem][ipHi][ipLo].ipFine = -1;
					}
				}
			}
		}
	}

	/* co carbon monoxide line array */
	for( i=0; i < nCORotate; i++ )
	{
		/* excitation energy of transition in degrees kelvin */
		C12O16Rotate[i].EnergyK = 
			(float)(T1CM)*C12O16Rotate[i].EnergyWN;
		C13O16Rotate[i].EnergyK = 
			(float)(T1CM)*C13O16Rotate[i].EnergyWN;

		/* energy of photon in ergs */
		C12O16Rotate[i].EnergyErg = 
			(float)(ERG1CM)*C12O16Rotate[i].EnergyWN;
		C13O16Rotate[i].EnergyErg = 
			(float)(ERG1CM)*C13O16Rotate[i].EnergyWN;

		/* put null terminated line label into chLab using line array*/
		chIonLbl(chLab, &C12O16Rotate[i]);
		chIonLbl(chLab, &C13O16Rotate[i]);

		C12O16Rotate[i].ipCont = 
			(int)ipLineEnergy(C12O16Rotate[i].EnergyWN * WAVNRYD, "12CO" ,0);
		C12O16Rotate[i].ipFine = 
			ipFineCont(C12O16Rotate[i].EnergyWN * WAVNRYD );
		C13O16Rotate[i].ipCont = 
			(int)ipLineEnergy(C13O16Rotate[i].EnergyWN * WAVNRYD, "13CO" ,0);
		C13O16Rotate[i].ipFine = 
			ipFineCont(C13O16Rotate[i].EnergyWN * WAVNRYD );

		if( C12O16Rotate[i].gf > 0. )
		{
			/* the gf value was entered
			 * derive the A, call to function is gf, wl (A), g_up */
			C12O16Rotate[i].Aul = 
				(float)(eina(C12O16Rotate[i].gf,
			  C12O16Rotate[i].EnergyWN,C12O16Rotate[i].gHi));
		}
		else if( C12O16Rotate[i].Aul > 0. )
		{
			/* the Einstion A value was entered
			 * derive the gf, call to function is A, wl (A), g_up */
			C12O16Rotate[i].gf = 
				(float)(GetGF(C12O16Rotate[i].Aul,
			  C12O16Rotate[i].EnergyWN,C12O16Rotate[i].gHi));
		}
		else
		{
			fprintf( ioQQQ, " 12CO line does not have valid gf or A\n" );
			fprintf( ioQQQ, " This is ContCreatePointers\n" );
			puts( "[Stop in ContCreatePointers]" );
			cdEXIT(EXIT_FAILURE);
		}
		if( C13O16Rotate[i].gf > 0. )
		{
			/* the gf value was entered
			 * derive the A, call to function is gf, wl (A), g_up */
			C13O16Rotate[i].Aul = 
				(float)(eina(C13O16Rotate[i].gf,
			  C13O16Rotate[i].EnergyWN,C13O16Rotate[i].gHi));
		}
		else if( C13O16Rotate[i].Aul > 0. )
		{
			/* the Einstion A value was entered
			 * derive the gf, call to function is A, wl (A), g_up */
			C13O16Rotate[i].gf = 
				(float)(GetGF(C13O16Rotate[i].Aul,
			  C13O16Rotate[i].EnergyWN,C13O16Rotate[i].gHi));
		}
		else
		{
			fprintf( ioQQQ, " 13CO line does not have valid gf or A\n" );
			fprintf( ioQQQ, " This is ContCreatePointers\n" );
			puts( "[Stop in ContCreatePointers]" );
			cdEXIT(EXIT_FAILURE);
		}

		/*line wavelength must be gt 0 */
		ASSERT( C12O16Rotate[i].WLAng > 0 );
		ASSERT( C13O16Rotate[i].WLAng > 0 );

		C12O16Rotate[i].dampXvel = 
			(float)(C12O16Rotate[i].Aul/
		  C12O16Rotate[i].EnergyWN/PI4);

		C13O16Rotate[i].dampXvel = 
			(float)(C13O16Rotate[i].Aul/
		  C13O16Rotate[i].EnergyWN/PI4);

		/* derive the abs coef, call to function is gf, wl (A), g_low */
		C12O16Rotate[i].opacity = 
			(float)(abscf(C12O16Rotate[i].gf,
		  C12O16Rotate[i].EnergyWN,C12O16Rotate[i].gLo));
		C13O16Rotate[i].opacity = 
			(float)(abscf(C13O16Rotate[i].gf,
		  C13O16Rotate[i].EnergyWN,C13O16Rotate[i].gLo));
	}

	/* inner shell transitions */
	for( i=0; i<nUTA; ++i )
	{
		if( UTALines[i].Aul > 0. )
		{
			UTALines[i].dampXvel = 
				(float)(UTALines[i].Aul/ UTALines[i].EnergyWN/PI4);

			/* derive the abs coef, call to function is gf, wl (A), g_low */
			UTALines[i].opacity = 
				(float)(abscf( UTALines[i].gf, UTALines[i].EnergyWN, UTALines[i].gLo));

			/* excitation energy of transition in degrees kelvin */
			UTALines[i].EnergyK = 
				(float)(T1CM*UTALines[i].EnergyWN);

			/* energy of photon in ergs */
			UTALines[i].EnergyErg = 
				(float)(ERG1CM*UTALines[i].EnergyWN);

			/* put null terminated line label into chLab using line array*/
			chIonLbl(chLab, &UTALines[i]);

			/* get pointer to energy in continuum mesh */
			UTALines[i].ipCont = (int)ipLineEnergy(UTALines[i].EnergyWN * WAVNRYD , chLab,0 );
			UTALines[i].ipFine = ipFineCont(UTALines[i].EnergyWN * WAVNRYD  );
		}
	}

	/* level 2 heavy element lines */
	for( i=0; i < nWindLine; i++ )
	{

		/* derive the A, call to function is gf, wl (A), g_up */
		TauLine2[i].Aul = 
			(float)(eina(TauLine2[i].gf,
		  TauLine2[i].EnergyWN,TauLine2[i].gHi));

		/* coef needed for damping constant */
		TauLine2[i].dampXvel = 
			(float)(TauLine2[i].Aul/
		  TauLine2[i].EnergyWN/PI4);

		/* derive the abs coef, call to function is gf, wl (A), g_low */
		TauLine2[i].opacity = 
			(float)(abscf(TauLine2[i].gf,
		  TauLine2[i].EnergyWN,TauLine2[i].gLo));

		/* excitation energy of transition in degrees kelvin */
		TauLine2[i].EnergyK = 
			(float)(T1CM*TauLine2[i].EnergyWN);

		/* energy of photon in ergs */
		TauLine2[i].EnergyErg = 
			(float)(ERG1CM*TauLine2[i].EnergyWN);

		/* put null terminated line label into chLab using line array*/
		chIonLbl(chLab, &TauLine2[i]);

		/* get pointer to energy in continuum mesh */
		TauLine2[i].ipCont = (int)ipLineEnergy(TauLine2[i].EnergyWN * WAVNRYD , chLab,0 );
		TauLine2[i].ipFine = ipFineCont(TauLine2[i].EnergyWN * WAVNRYD );
		/*if( TauLine2[i].ipCont==751 )
			fprintf(ioQQQ,"( atom_level2 %s\n", chLab);*/
	}

	/* hyperfine structure lines */
	for( i=0; i < nHFLines; i++ )
	{
		/* derive the A, call to function is gf, wl (A), g_up 
		HFLines[i].Aul = 
			(float)(eina(HFLines[i].gf,
		  HFLines[i].EnergyWN,HFLines[i].gHi));*/

		/* coef needed for damping constant */
		HFLines[i].dampXvel = 
			(float)(HFLines[i].Aul/
			HFLines[i].EnergyWN/PI4);

		/* derive the abs coef, call to function is gf, wl (A), g_low */
		HFLines[i].opacity = 
			(float)(abscf(HFLines[i].gf,
			HFLines[i].EnergyWN,
			HFLines[i].gLo));
		/* gf from this and 21 cm do not agree, A for HFS is 10x larger than level1 dat */
		/*fprintf(ioQQQ,"HFLinesss\t%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
			i,HFLines[i].opacity , HFLines[i].gf , HFLines[i].Aul , HFLines[i].EnergyWN,HFLines[i].gLo);*/

		/* excitation energy of transition in degrees kelvin */
		HFLines[i].EnergyK = 
			(float)(T1CM*HFLines[i].EnergyWN);

		/* energy of photon in ergs */
		HFLines[i].EnergyErg = 
			(float)(ERG1CM*HFLines[i].EnergyWN);

		/* put null terminated line label into chLab using line array*/
		chIonLbl(chLab, &HFLines[i]);

		/* get pointer to energy in continuum mesh */
		HFLines[i].ipCont = (int)ipLineEnergy(HFLines[i].EnergyWN * WAVNRYD , chLab,0 );
		HFLines[i].ipFine = ipFineCont(HFLines[i].EnergyWN * WAVNRYD );
		/*if( HFLines[i].ipCont==751 )
			fprintf(ioQQQ,"( atom_level2 %s\n", chLab);*/
	}


	/* Verner's FeII lines - do first so following labels will over write this
	 * only call if large feii atom is turned on */
	if( FeII.lgFeIION )
	{
		FeIIPoint();
	}

	/* the group of inner shell fluorescent lines */
	for( i=0; i<yield.nfl_lines; ++i )
	{
		strcpy( chLab , elementnames.chElementSym[yield.nfl_nelem[i]] );
		strcat( chLab , elementnames.chIonStage[yield.nfl_ion_emit[i]] );
		
		yield.nfl_ipoint[i] = ipLineEnergy( yield.fl_energy[i] , chLab , 0 );
	}

	/* ================================================================================== */
	/*        two photon two-photon 2-nu 2 nu 2 photon 2-photon                           */

	/* >>chng 01 jan 27, add iso-electronic dimension to following vectors */
	/* for induced two photon emission we need mirror image set
	 * of continuum indices for continuum between Lya and helf Lya */
	/* first MALLOC LIMELM dimension of space */
	/* >>chng 02 jun 28, malloc this NISO instead of 2.	*/
	if( (iso.ipSym2nu = (long ***)MALLOC(sizeof(long *)*(unsigned)NISO ) )==NULL )
		BadMalloc();

	if( (iso.As2nu = (float ***)MALLOC(sizeof(float *)*(unsigned)NISO ) )==NULL )
		BadMalloc();

	/* now loop over the two iso-sequences with two photon continua */
	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	{
		if( (iso.ipSym2nu[ipISO] = (long **)MALLOC(sizeof(long *)*(unsigned)LIMELM ) )==NULL )
			BadMalloc();
		if( (iso.As2nu[ipISO] = (float **)MALLOC(sizeof(float *)*(unsigned)LIMELM ) )==NULL )
			BadMalloc();
	
		/* set up two photon emission */
		for( nelem=ipISO; nelem<LIMELM; ++nelem )
		{
			if( dense.lgElmtOn[nelem] )
			{
				double E2nu = EmisLines[ipISO][nelem][TwoS][0].EnergyWN * WAVNRYD;
				double Aul = EmisLines[ipISO][nelem][TwoS][0].Aul;
				double SumShapeFunction = 0., Renorm= 0.;
				
				/* pointer to the Lya energy */
				iso.ipTwoPhoE[ipISO][nelem] = ipoint(E2nu);
				/* this half-energy is only used to get induced rates in twophoton */
				iso.ipHalfTwoPhoE[ipISO][nelem] = ipoint(E2nu / 2.);
				while( rfield.anu[iso.ipTwoPhoE[ipISO][nelem]] > E2nu )
				{
					--iso.ipTwoPhoE[ipISO][nelem];
				}
				while( rfield.anu[iso.ipHalfTwoPhoE[ipISO][nelem]] > E2nu / 2. )
				{
					--iso.ipHalfTwoPhoE[ipISO][nelem];
				}
				
				if( (iso.ipSym2nu[ipISO][nelem] = (long *)MALLOC(sizeof(long)*(unsigned)iso.ipTwoPhoE[ipISO][nelem] ) )==NULL )
					BadMalloc();
				if( (iso.As2nu[ipISO][nelem] = (float *)MALLOC(sizeof(float)*(unsigned)iso.ipTwoPhoE[ipISO][nelem] ) )==NULL )
					BadMalloc();

				/* >>chng 02 aug 14, change upper limit to full Lya energy */
				for( i=0; i < iso.ipTwoPhoE[ipISO][nelem]; i++ )
				{
					/* energy is symmetric energy, the other side of half E2nu */
					energy = E2nu - rfield.anu[i];
					/* this is needed since mirror image of cell next to two-nu energy
					 * may be slightly negative */
					energy = MAX2( energy, rfield.anu[0] + rfield.widflx[0]/2. );
					/* find index for this symmetric energy */
					iso.ipSym2nu[ipISO][nelem][i] = ipoint(energy);
					while( rfield.anu[iso.ipSym2nu[ipISO][nelem][i]] > E2nu ||
						iso.ipSym2nu[ipISO][nelem][i] >= iso.ipTwoPhoE[ipISO][nelem])
					{
						--iso.ipSym2nu[ipISO][nelem][i];
					}
					ASSERT( iso.ipSym2nu[ipISO][nelem][i] >= 0 );
				}

				/* ipTwoPhoE is the cell holding the 2nu energy itself, and we do not want
				 * to include that in the following sum */
				for( i=0; i<iso.ipTwoPhoE[ipISO][nelem] ; i++ )
				{
					double ShapeFunction;

					ASSERT( rfield.anu[i]<=E2nu );
					
					ShapeFunction = atmdat_2phot_shapefunction( rfield.AnuOrg[i]/E2nu, ipISO, nelem )*rfield.widflx[i]/E2nu;
					SumShapeFunction += ShapeFunction;

					/* >>refer	HI	2nu	Spitzer, L., & Greenstein, J., 1951, ApJ, 114, 407 */
					/* As2nu will add up to the A, so its units are s-1	*/ 
					iso.As2nu[ipISO][nelem][i] = (float)( Aul * ShapeFunction );
					/* the factor of widflx is not a mistake - when used for induced two-photon
					 * emission the sum will only go up to half of Lya, so only one factor
					 * of widflx will ever appear in the sum */
#if	0
					ASSERT( iso.As2nu[ipISO][nelem][i] >= 0.f );

					if( ipISO == ipH_LIKE && nelem == ipHYDROGEN )
					{
						fprintf(ioQQQ, "%.3e\t%.3e\t%.3e\n",rfield.anu[i],
							rfield.widflx[i], iso.As2nu[ipISO][nelem][i]);
					}
#endif
				}

				/* The spline function in twophoton.c causes a bit of an error in the integral of the
				 * shape function.  So we renormalize the integral to 1.	*/
				Renorm = 1./SumShapeFunction;
				
				for( i=0; i<iso.ipTwoPhoE[ipISO][nelem] ; i++ )
				{
					iso.As2nu[ipISO][nelem][i] *= (float)Renorm;
				}
				
				/* The result should be VERY close to 1.	*/
				ASSERT( fabs( SumShapeFunction*Renorm - 1. ) < 0.00001 );
			}
		}
	}

	{
		/* this is an option to print out one of the two photon continua */
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC )
		{	
#			define	NCRS	21
			double limit,ener[NCRS]={
			  0.,     0.03738,  0.07506,  0.1124,  0.1498,  0.1875,
			  0.225,  0.263,    0.300,    0.3373,  0.375,   0.4127,
			  0.4500, 0.487,    0.525,    0.5625,  0.6002,  0.6376,
			  0.6749, 0.7126,   0.75};

			nelem = ipHYDROGEN;
			ipISO = ipHYDROGEN;
			
			limit = iso.ipTwoPhoE[ipISO][nelem];
			
			for( i=0; i < NCRS; i++ )
			{
				fprintf(ioQQQ,"%.3e\t%.3e\n", ener[i] , 
					atmdat_2phot_shapefunction( ener[i]/0.75, ipISO, nelem ) );
			}

			xnew = 0.;
			/* TODO	2	what are we trying to print here?	*/
			for( i=0; i < limit; i++ )
			{
				double fach = iso.As2nu[ipISO][nelem][i]/2.*rfield.anu2[i]/rfield.widflx[i]*EN1RYD;
				fprintf(ioQQQ,"%.3e\t%.3e\t%.3e\n", 
					rfield.anu[i] , 
					iso.As2nu[ipISO][nelem][i] / rfield.widflx[i] , 
					fach );
				xnew += iso.As2nu[ipISO][nelem][i];
			}
			fprintf(ioQQQ," sum is %.3e\n", xnew );
			puts( "[Stop in ContCreatePointers]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* ================================================================================== */


	{
		/* this is an option to print out one of the two photon continua */
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC )
		{	
			for( i=0; i<11; ++i )
			{
				char chLsav[10];
				TauDummy.WLAng = (float)(PI * pow(10.,i));
				strcpy( chLsav, chLineLbl(&TauDummy) );
				fprintf(ioQQQ,"%.2f\t%s\n", TauDummy.WLAng , chLsav );
			}
			exit(0);
		}
	}

	/* option to print out whole thing with "trace lines" command */
	if( trace.lgTrLine )
	{
		fprintf( ioQQQ, "       WL(Ang)   E(RYD)   IP   gl  gu      gf       A        damp     abs K\n" );
		for( i=1; i <= nLevel1; i++ )
		{
			strcpy( chLab, chLineLbl(&TauLines[i]) );
			iWL_Ang = (long)TauLines[i].WLAng;
			if( iWL_Ang > 1000000 )
			{
				iWL_Ang /= 10000;
			}
			else if( iWL_Ang > 10000 )
			{
				iWL_Ang /= 1000;
			}

			fprintf( ioQQQ, " %10.10s%5ld%10.3e %4i%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n", 
			  chLab, iWL_Ang, RYDLAM/TauLines[i].WLAng, 
			  TauLines[i].ipCont, (long)(TauLines[i].gLo), 
			  (long)(TauLines[i].gHi), TauLines[i].gf, 
			  TauLines[i].Aul, TauLines[i].dampXvel, 
			  TauLines[i].opacity );
		}

		/* C12O16 lines */
		for( i=0; i < nCORotate; i++ )
		{
			strcpy( chLab, chLineLbl(&C12O16Rotate[i]) );

			iWL_Ang = (long)C12O16Rotate[i].WLAng;

			if( iWL_Ang > 1000000 )
			{
				iWL_Ang /= 10000;
			}
			else if( iWL_Ang > 10000 )
			{
				iWL_Ang /= 1000;
			}
			fprintf( ioQQQ, " %10.10s%5ld%10.3e %4i%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n", 
			  chLab, iWL_Ang, RYDLAM/C12O16Rotate[i].WLAng, 
			  C12O16Rotate[i].ipCont, (long)(C12O16Rotate[i].gLo), 
			  (long)(C12O16Rotate[i].gHi), C12O16Rotate[i].gf, 
			  C12O16Rotate[i].Aul, C12O16Rotate[i].dampXvel, 
			  C12O16Rotate[i].opacity );
		}

		/* C13O16 lines */
		for( i=0; i < nCORotate; i++ )
		{
			strcpy( chLab, chLineLbl(&C13O16Rotate[i]) );

			iWL_Ang = (long)C13O16Rotate[i].WLAng;

			if( iWL_Ang > 1000000 )
			{
				iWL_Ang /= 10000;
			}
			else if( iWL_Ang > 10000 )
			{
				iWL_Ang /= 1000;
			}
			fprintf( ioQQQ, " %10.10s%5ld%10.3e %4i%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n", 
			  chLab, iWL_Ang, RYDLAM/C13O16Rotate[i].WLAng, 
			  C13O16Rotate[i].ipCont, (long)(C13O16Rotate[i].gLo), 
			  (long)(C13O16Rotate[i].gHi), C13O16Rotate[i].gf, 
			  C13O16Rotate[i].Aul, C13O16Rotate[i].dampXvel, 
			  C13O16Rotate[i].opacity );
		}

		for( i=0; i < nWindLine; i++ )
		{
			strcpy( chLab, chLineLbl(&TauLine2[i]) );

			iWL_Ang = (long)TauLine2[i].WLAng;

			if( iWL_Ang > 1000000 )
			{
				iWL_Ang /= 10000;
			}
			else if( iWL_Ang > 10000 )
			{
				iWL_Ang /= 1000;
			}
			fprintf( ioQQQ, " %10.10s%5ld%10.3e %4i%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n", 
			  chLab, iWL_Ang, RYDLAM/TauLine2[i].WLAng, 
			  TauLine2[i].ipCont, (long)(TauLine2[i].gLo), 
			  (long)(TauLine2[i].gHi), TauLine2[i].gf, 
			  TauLine2[i].Aul, TauLine2[i].dampXvel, 
			  TauLine2[i].opacity );
		}
		for( i=0; i < nHFLines; i++ )
		{
			strcpy( chLab, chLineLbl(&HFLines[i]) );

			iWL_Ang = (long)HFLines[i].WLAng;

			if( iWL_Ang > 1000000 )
			{
				iWL_Ang /= 10000;
			}
			else if( iWL_Ang > 10000 )
			{
				iWL_Ang /= 1000;
			}
			fprintf( ioQQQ, " %10.10s%5ld%10.3e %4i%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n", 
			  chLab, iWL_Ang, RYDLAM/HFLines[i].WLAng, 
			  HFLines[i].ipCont, (long)(HFLines[i].gLo), 
			  (long)(HFLines[i].gHi), HFLines[i].gf, 
			  HFLines[i].Aul, HFLines[i].dampXvel, 
			  HFLines[i].opacity );
		}
	}

	/* this is an option to kill fine structure line optical depths */
	if( !rt.lgFstOn )
	{
		for( i=1; i <= nLevel1; i++ )
		{
			if( TauLines[i].EnergyWN < 10000. )
			{
				TauLines[i].opacity = 0.;
			}
		}
	}

#	ifdef DEBUG_FUN
	fputs( " <->ContCreatePointers()\n", debug_fp );
#	endif
	return;
}

/*fiddle adjust energy bounds of cell with index ipnt so that lower energy
 * matchs ionization edges exactly, called by createpoint */
/* ipnt is index on c scale */
static void fiddle(long int ipnt, 
  double exact)
{
	float Ehi, 
	  Elo,
	  OldEner;

#	ifdef DEBUG_FUN
	fputs( "<+>fiddle()\n", debug_fp );
#	endif

	/* make low edge of cell exact for photo integrals */
	ASSERT( ipnt >= 0 );
	ASSERT( ipnt < rfield.nupper-1 );
	
	/* upper edge of higher cell*/
	Ehi = rfield.anu[ipnt] + rfield.widflx[ipnt]*0.5f;
	/* lower edge of lower cell */
	Elo = rfield.anu[ipnt-1] - rfield.widflx[ipnt-1]*0.5f;

	/* >>chng 02 nov 11, do nothing if already very close,
	 * because VERY high res models had negative widths for some very close edges */
	if( fabs( Elo/exact - 1. ) < 0.001 ) 
		return;

	ASSERT( Elo <= exact );

	OldEner = rfield.anu[ipnt];

	/* centroid of desired lower energy and upper edge */
	rfield.anu[ipnt] = (float)((Ehi + exact)/2.);
	/* same for lower */
	rfield.anu[ipnt-1] = (float)((exact + Elo)/2.);

	rfield.widflx[ipnt] = (float)(Ehi - exact);
	rfield.widflx[ipnt-1] = (float)(exact - Elo);

	/* bring upper cell down a bit too, since we dont want large change in widflx */
	rfield.anu[ipnt+1] -= (OldEner - rfield.anu[ipnt])/2.f;

	/* sanity check */
	ASSERT( rfield.widflx[ipnt-1] > 0. );
	ASSERT( rfield.widflx[ipnt] > 0. );

#if	0
	/* This code spreads the tweak over some number of cells to each side of ipnt...
	 * May not be useful as it could cause errors in ipoint, and may be unimportant anyway.	*/

	/* gap is amount by which to tweak cell ipnt, SHIFT says 
	 * how many cells to each side of ipnt should be nudged.  */
	long n,SHIFT = 1;

	float gap;

	/* must make sure there is room to each side	*/
	ASSERT( ipnt >= SHIFT );
	ASSERT( ipnt < rfield.nupper-1-SHIFT );

	/* Amount by which to tweak cell ipnt	*/
	gap = rfield.anu[ipnt]-rfield.widflx[ipnt]*0.5f-(float)exact;
	
	/* We will adjust 2*SHIFT+1 cells, SHIFT each way plus ipnt itself...
	 * save upper and lower edge of this region to make sure they stay the same.	*/
	Elo = rfield.anu[ipnt-SHIFT] - rfield.widflx[ipnt-SHIFT]*0.5f;
	Ehi = rfield.anu[ipnt+SHIFT] + rfield.widflx[ipnt+SHIFT]*0.5f;

	if( gap > 0.f )
	{
		/* First slightly shorten each of the previous SHIFT cells by one nth of gap	*/
		for( n=0; n<=SHIFT-1; n++ )
		{
			rfield.anu[ipnt-SHIFT+n] -= gap*(2.f*n+1.f)/(2.f*SHIFT);
			rfield.widflx[ipnt-SHIFT+n] -= gap/(float)SHIFT;
		}

		/* We'll keep the ipnt cell itself the same width...so it must be shifted down by gap.	*/
		rfield.anu[ipnt] -= gap;

		/* Now slightly enlarge the next SHIFT cells by one nth of gap.	*/
		for( n=0; n<=SHIFT-1; n++ )
		{
			rfield.anu[ipnt+1+n] -= gap*(2.f*(SHIFT-1.f-n)+1.f)/(2.f*SHIFT);
			rfield.widflx[ipnt+1+n] += gap/(float)SHIFT;
		}
	}
	else if( gap < 0.f )
	{
		/* First slightly widen each of the previous SHIFT cells by one nth of gap	*/
		for( n=0; n<=SHIFT-1; n++ )
		{
			rfield.anu[ipnt-SHIFT+n] += gap*(2.f*n+1.f)/(2.f*SHIFT);
			rfield.widflx[ipnt-SHIFT+n] += gap/(float)SHIFT;
		}

		/* We'll keep the ipnt cell itself the same width...so it must be shifted up by gap.	*/
		rfield.anu[ipnt] += gap;

		/* Now slightly shorten the next SHIFT cells by one nth of gap.	*/
		for( n=0; n<=SHIFT-1; n++ )
		{
			rfield.anu[ipnt+1+n] += gap*(2.f*(SHIFT-1.f-n)+1.f)/(2.f*SHIFT);
			rfield.widflx[ipnt+1+n] -= gap/(float)SHIFT;
		}
	}

	/* Did it work correctly?	*/
	ASSERT( fabs( 1.f - (rfield.anu[ipnt-SHIFT] - rfield.widflx[ipnt-SHIFT]*0.5f)/Elo ) < 0.0001f );
	ASSERT( fabs( 1.f - (rfield.anu[ipnt+SHIFT] + rfield.widflx[ipnt+SHIFT]*0.5f)/Ehi ) < 0.0001f );
#endif
	
#	ifdef DEBUG_FUN
	fputs( " <->fiddle()\n", debug_fp );
#	endif
	return;
}

/*ipShells assign continuum energy pointers to shells for all atoms,
 * called by ContCreatePointers */
static void ipShells(
	/* nelem is the atomic number on the C scale, Li is 2 */
	long int nelem)
{
	char chLab[5];
	long int 
	  imax, 
	  ion, 
	  nelec, 
	  ns, 
	  nshell;
	/* following value cannot be used - will be set to proper threshold */
	double thresh=-DBL_MAX;

#	ifdef DEBUG_FUN
	fputs( "<+>ipShells()\n", debug_fp );
#	endif

	ASSERT( nelem >= 2);
	ASSERT( nelem < LIMELM );

	/* fills in pointers to valence shell ionization threshold
	 * PH1(a,b,c,d)
	 * a=1 => thresh, others fitting parameters
	 * b atomic number
	 * c number of electrons
	 * d shell number 7-1 */

	/* threshold in Ryd
	 * ion=0 for atom, up to nelem-1 for helium like, hydrogenic is elsewhere */
	for( ion=0; ion < nelem; ion++ )
	{
		/* generate label for ionization stage */
		/* this is short form of element name */
		strcpy( chLab, elementnames.chElementSym[nelem] );

		/* this is a number between 1 and 31 */
		strcat( chLab, elementnames.chIonStage[ion] );

		/* number of bound electrons */
		nelec = nelem+1 - ion ;

		/* nsShells(nelem,ion) is the number of shells for ion with nelec electrons,
		 * physical not c scale */
		imax = Heavy.nsShells[nelem][ion];

		/* loop on all inner shells, valence shell */
		for( nshell=0; nshell < imax; nshell++ )
		{
			/* ionization potential of this shell in rydbergs */
			thresh = (double)(PH1COM.PH1[nshell][nelec-1][nelem][0]/EVRYD* 0.9998787);
			if( thresh <= 0.1 )
			{
				/* negative ip shell does not exist, set upper limit
				 * to less than lower limit so this never looped upon
				 * these are used as flags by LimitSh to check whether
				 * this is a real shell - if 1 or 2 is changed - change LimitSh!! */
				opac.ipElement[nelem][ion][nshell][0] = 2;
				opac.ipElement[nelem][ion][nshell][1] = 1;
			}
			else
			{
				/* this is lower dbound to energy range for this shell */
				/* >>chng 02 may 27, change to version of ip with label, so that
				 * inner shell edges will appear */
				/*opac.ipElement[nelem][ion][nshell][0] = ipoint(thresh);*/
				opac.ipElement[nelem][ion][nshell][0] = 
					ipContEnergy( thresh , chLab );

				/* this is upper bound to energy range for this shell 
				 * LimitSh is an integer function, returns pointer
				 * to threshold of next major shell.  For k-shell it
				 * returns the values KshellLimit, default=7.35e4
				 * >>chng 96 sep 26, had been below, result zero cross sec at 
				 * many energies where opacity project did not produce state specific 
				 * cross section */
				opac.ipElement[nelem][ion][nshell][1] = 
					LimitSh(ion+1,  nshell+1,nelem+1);
				ASSERT( opac.ipElement[nelem][ion][nshell][1] > 0);
			}
		}

		/* this will be index pointing to valence edge */
		/* [0] is pointer to threshold in energy array */
		opac.ipElement[nelem][ion][imax-1][0] = 
			ipContEnergy(thresh, chLab);

		/* pointer to valence electron ionization potential */
		Heavy.ipHeavy[nelem][ion] = opac.ipElement[nelem][ion][imax-1][0];

		/* ionizaton potential of valence shell in ryd */
		/* thresh was eval above, now has last value, the valence shell */
		Heavy.Valence_IP_Ryd[nelem][ion] = thresh;

		/* this is set of 3/4 of valence shell IP, this is important
		 * source of ots deep in cloud */
		Heavy.ipLyHeavy[nelem][ion] = 
			ipLineEnergy(thresh*0.75,chLab , 0);
		Heavy.ipBalHeavy[nelem][ion] = 
			ipLineEnergy(thresh*0.25,chLab , 0);
	}

	/* above loop did up to hydrogenic, now do hydrogenic - 
	 * hydrogenic is special since arrays already set up */
	Heavy.nsShells[nelem][nelem] = 1;

	/* this is lower limit to range */
	/* hydrogenic photoionization set to special hydro array 
	 * this is pointer to threshold energy */
	/* this statement is in ContCreatePointers but has not been done when this routine called */
	/*iso.ipIsoLevNIonCon[ipH_LIKE][ipZ][ipLo] = ipContEnergy(iso.xIsoLevNIonRyd[ipH_LIKE][ipZ][ipLo],chLab);*/
	/*opac.ipElement[nelem][nelem][0][0] = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s];*/
	opac.ipElement[nelem][nelem][0][0] = ipoint( iso.xIsoLevNIonRyd[ipH_LIKE][nelem][ipH1s] );
	ASSERT( opac.ipElement[nelem][nelem][0][0] > 0 );

	/* this is the high energy limit */
	opac.ipElement[nelem][nelem][0][1] = KshllEnr.KshellLimit;

	Heavy.ipHeavy[nelem][nelem] = opac.ipElement[nelem][nelem][0][0];

	/* this is for backwards compatability with cambridge code */
	if( trace.lgTrace && trace.lgPointBug )
	{
		for( ion=0; ion < (nelem+1); ion++ )
		{
			fprintf( ioQQQ, "Ion:%3ld%3ld %2.2s%2.2s total shells:%3ld\n", 
			  nelem, ion+1, elementnames.chElementSym[nelem], elementnames.chIonStage[ion]
			  , Heavy.nsShells[nelem][ion] );
			for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
			{
				fprintf( ioQQQ, " shell%3ld %2.2s range eV%10.2e-%8.2e\n", 
				  ns+1, Heavy.chShell[ns], rfield.anu[opac.ipElement[nelem][ion][ns][0]-1]*
				  EVRYD, rfield.anu[opac.ipElement[nelem][ion][ns][1]-1]*EVRYD );
			}
		}
	}

#	ifdef DEBUG_FUN
	fputs( " <->ipShells()\n", debug_fp );
#	endif
	return;
}

/*LimitSh sets upper energy limit to subshell integrations */
static long LimitSh(long int ion, 
  long int nshell, 
  long int nelem)
{
	long int LimitSh_v;

#	ifdef DEBUG_FUN
	fputs( "<+>LimitSh()\n", debug_fp );
#	endif

	/* this routine returns the high energy limit to the energy range
	 * for photoionizatin of a given shell
	 * */
	if( nshell == 1 )
	{
		/* this limit is high energy limit to code unless changed with set kshell */
		LimitSh_v = KshllEnr.KshellLimit;

	}
	else if( nshell == 2 )
	{
		/* this is 2s shell, upper limit is 1s
		 * >>chng 96 oct 08, up to high energy limit
		 * LimitSh = ipElement(nelem,ion , 1,1)-1 */
		LimitSh_v = KshllEnr.KshellLimit;

	}
	else if( nshell == 3 )
	{
		/* this is 2p shell, upper limit is 1s
		 * >>chng 96 oct 08, up to high energy limit
		 * LimitSh = ipElement(nelem,ion , 1,1)-1 */
		LimitSh_v = KshllEnr.KshellLimit;

	}
	else if( nshell == 4 )
	{
		/* this is 3s shell, upper limit is 2p
		 * >>chng 96 oct 08, up to K-shell edge
		 * LimitSh = ipElement(nelem,ion , 3,1)-1 */
		LimitSh_v = opac.ipElement[nelem-1][ion-1][0][0] - 1;

	}
	else if( nshell == 5 )
	{
		/* this is 3p shell, upper limit is 2p
		 * >>chng 96 oct 08, up to K-shell edge
		 * LimitSh = ipElement(nelem,ion , 3,1)-1 */
		LimitSh_v = opac.ipElement[nelem-1][ion-1][0][0] - 1;

	}
	else if( nshell == 6 )
	{
		/* this is 3d shell, upper limit is 2p
		 * >>chng 96 oct 08, up to K-shell edge
		 * LimitSh = ipElement(nelem,ion , 3,1)-1 */
		LimitSh_v = opac.ipElement[nelem-1][ion-1][0][0] - 1;

	}
	else if( nshell == 7 )
	{
		/* this is 4s shell, upper limit is 3d */
		if( opac.ipElement[nelem-1][ion-1][5][0] < 3 )
		{
			/* this is check for empty shell 6, 3d
			 * if so then set to 3p instead */
			LimitSh_v = opac.ipElement[nelem-1][ion-1][4][0] - 
			  1;
		}
		else
		{
			LimitSh_v = opac.ipElement[nelem-1][ion-1][5][0] - 
			  1;
		}
		/* >>chng 96 sep 26, set upper limit down to 2s */
		LimitSh_v = opac.ipElement[nelem-1][ion-1][2][0] - 1;

	}
	else
	{
		fprintf( ioQQQ, " LimitSh cannot handle nshell as large as%4ld\n", 
		  nshell );
		puts( "[Stop in limitsh]" );
		cdEXIT(EXIT_FAILURE);
	}

#	ifdef DEBUG_FUN
	fputs( " <->LimitSh()\n", debug_fp );
#	endif
	return( LimitSh_v );
}

