/*CreatePoint set up pointers for lines and continua called by cloudy after input read in */
#include <string.h>
#include <limits.h>
#include "cddefines.h"
#include "hydrogenic.h"
#include "secondaries.h"
#include "nfeii.h"
#include "nhe1lvl.h"
#include "taulines.h"
#include "opacpoint.h"
#include "physconst.h"
#include "elementnames.h"
#include "chionlbl.h"
#include "getgf.h"
#include "abscf.h"
#include "fston.h"
#include "eina.h"
#include "chlinelbl.h"
#include "ipvfil.h"
#include "ipfe2.h"
#include "fe2pmp.h"
#include "elmton.h"
#include "xry.h"
#include "pumpfe.h"
#include "hhe2phtots.h"
#include "pbowen.h"
#include "ip626.h"
#include "bounds.h"
#include "hehp.h"
#include "ircoil.h"
#include "rfield.h"
#include "nh.h"
#include "he1nionryd.h"
#include "iphe1l.h"
#include "nhe1.h"
#include "he.h"
#include "ip2gam.h"
#include "i3727.h"
#include "pop371.h"
#include "ipoi.h"
#include "fekems.h"
#include "trace.h"
#include "iphmin.h"
#include "nxray.h"
#include "egray.h"
#include "kshllenr.h"
#include "fillnu.h"
#include "linelabl.h"
#include "ip2phtcom.h"
#include "heavy.h"
#include "he3lines.h"
#include "nsshells.h"
#include "predcont.h"
#include "ph1com.h"
#include "ipconsafe.h"
#include "iplinsafe.h"
#include "ipoint.h"
#include "ipshells.h"
#include "chckfill.h"
#include "fill.h"
#include "fiddle.h"
#include "createpoint.h"

void CreatePoint(void)
{
	char chLab[5];
	long int 
	  i, 
	  ipHi, 
	  ipLo, 
	  ipZ, 
	  ipnt, 
	  iWavLen,
	  j, 
	  n0, 
	  nelem,
	  nshells;
	double energy;
	/* flag to say whether pointers have ever been evaluated */
	static int lgPntEval = FALSE;

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

	/* lgPntEval is local static variable defined FALSE when defined. 
	 * it is set TRUE below, so that pointers only created one time in the
	 * history of this coreload. */
	if( lgPntEval )
	{
		if( trace.lgTrace )
		{
			fprintf( ioQQQ, " CreatePoint called, not evaluating.\n" );
		}
		/* now save current form of energy array */
		for( i=0; i < NCELL; i++ )
		{
			rfield.anu[i] = rfield.AnuOrg[i];
		}
		
#		ifdef DEBUG_FUN
		fputs( " <->CreatePoint()\n", debug_fp );
#		endif
		return;
	}
	else
	{
		if( trace.lgTrace )
		{
			fprintf( ioQQQ, " CreatePoint called first time.\n" );
		}
		lgPntEval = TRUE;
	}

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

	/* create the hydrogen atom for this coreload */
	HydroCreate();

	/* fill in continuum with freq points
	 * arg are range pointer, 2 energy limits, resolution
	 * first argument is lower energy of the continuum range to be defined
	 * second argument is upper range; both are in Rydbergs
	 * third number is the relative energy resolution, dnu/nu,
	 * for that range of the continuum
	 * last two numbers are internal book keeping
	 * N0 is number of energy cells used so far
	 * IPNT is a counter for the number of fills used
	 *
	 * if this is changed, then also change warning in GetTable about using
	 * transmitted continuum - it says version number where continuum changed
	 * */
	n0 = 1;
	ipnt = 0;
	fillnu.nrange = 0;/* this is number of ranges that will be introduced below*/
	/* >>>chng 99 Oct 13, need finer resolution in FIR */
	/*fill(bounds.emm, 1e-3,  0.1,&n0,&ipnt);*/
	/*fill(1e-3      , 0.05,  0.02,&n0,&ipnt);*/
	fill(bounds.emm, 1e-4,  0.1,&n0,&ipnt);
	fill(1e-4      , 0.05,  0.02,&n0,&ipnt);
	fill(     0.05 , 600.,  0.01,&n0,&ipnt);
	fill(600.,bounds.egamry, 0.03,&n0,&ipnt);
	/* at this point debugger shows that anu and widflx are defined 
	 * up through n0-2 - due to c offset -1 from fortran! */
	rfield.nupper = n0 - 1;
	rfield.nflux = rfield.nupper;

	/* there must be a cell above nflux for us to pass unity through the vol integrator
	 * as a sanity check.  assert that this is true so we will crash if ever changed */
	assert( rfield.nupper +1 <= NCELL );

	/* finish off rest of array, even if not used, N0 is next element to fill
	 * NUPPER is number of NCELL cells actually used for freq mesh */
	if( n0 <= NCELL )
	{
		for( i=rfield.nupper; i < NCELL; i++ )
		{
			rfield.widflx[i] = rfield.widflx[i-1];
			rfield.anu[i] = rfield.anu[i-1] + rfield.widflx[i];
		}
	}

	/* this is a sanity check for results produced above by fill */
	ChckFill();

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

	/* pointer to excited state of O+2 */
	OpacPoint.ipo3exc[0] = ipConSafe(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 */
	for( ipZ=0; ipZ < 2; ipZ++ )
	{
		/* generate label for this ion */
		sprintf( chLab, "%2s%2ld",elementnames.chElementSym[ipZ], ipZ+1 );
		/* Lya itself is the only transition below n=3 - we explicitly do not
		 * want to generate pointers for 2s-1s or 2p-2s */
		HydroLines[ipZ][IP2P][IP1S].ipCont = 
			ipLinSafe(HydroLines[ipZ][IP2P][IP1S].EnergyRyd , chLab , 0);
		/* set large negative pointers for 2p-2s and 2s-1s so we blow
		 * if we try to address 2s-2p or 1s-2s */
		HydroLines[ipZ][IP2S][IP1S].ipCont = INT_MIN;
		/* 2p-2s set to zero so that finite but not important */
		HydroLines[ipZ][IP2P][IP2S].ipCont = 1;
		for( ipLo=IP1S; ipLo <= nhlevel; ipLo++ )
		{
			/* HLevNIonRyd is level energy ryd 
			 * this is pointer to continuum ionization threshold 
			 * also puts label into continuum array */
			nh.ipHn[ipZ][ipLo] = ipConSafe(hydro.HLevNIonRyd[ipZ][ipLo],chLab);

			/* define all hydrogenic line pointers, MAX2 is to avoid creating
			 * pointers for two photon emission or 2s-2p transition */
			for( ipHi=MAX2(3,ipLo+1); ipHi <= nhlevel; ipHi++ )
			{
				HydroLines[ipZ][ipHi][ipLo].ipCont = 
					ipLinSafe(HydroLines[ipZ][ipHi][ipLo].EnergyRyd , chLab ,
					nh.ipHn[ipZ][ipLo]);

				/* do not let pointer go into next higher continuum */
				/*HydroLines[ipZ][ipHi][ipLo].ipCont = 
					MIN2( j , (nh.ipHn[ipZ][ipLo]-1) );*/
			}
		}
		/* set pointer to two photon continuum as emission line to 
		 * completely bad value - if ever used should cause code to crash */
		HydroLines[ipZ][IP2S][IP1S].ipCont = INT_MIN;
	}

	/* need to increase the cell for the HeII balmer continuum by one, so that
	 * hydrogen lyman continuum will pick it up. */
	ipZ = 1;
	/* copy label over to next higher cell */
	strcpy( LineLabl.chContLabel[nh.ipHn[ipZ][IP2S]], 
		     LineLabl.chContLabel[nh.ipHn[ipZ][IP2S]-1] );
	/* set previous spot to blank */
	strcpy( LineLabl.chContLabel[nh.ipHn[ipZ][IP2S]-1] , "    ");
	/* finally increment the two pointers */
	++nh.ipHn[ipZ][IP2S];
	++nh.ipHn[ipZ][IP2P];

	/* 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 */
	OpacPoint.ih2pnt[0] = ipConSafe(0.502,"H2+ ");
	OpacPoint.ih2pnt[1] = ipoint(1.03);

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

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

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

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

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

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


	/* fix pointers for hydrogen and helium
	 *
	 * following are pointers to places to save locally destroyed
	 * two photon continua */
	hhe2phtots.ipH12pht = ipConSafe(0.6,"H1p2");
	hhe2phtots.ipHe12pht = ipConSafe(1.2,"He12");
	hhe2phtots.ipHe22pht = ipConSafe(1.5,"He22");

	/* pointers to He I triplet lines */
	he3lines.ipHe3l[0] = ipLinSafe(0.08416,"HeI3",0);
	he3lines.ipHe3l[1] = ipLinSafe(RYDLAM/5876.,"HeI3",0);
	he3lines.ipHe3l[2] = ipLinSafe(RYDLAM/7065.,"HeI3",0);
	he3lines.ipHe3l[3] = ipLinSafe(RYDLAM/3889.,"HeI3",0);

	/* pointer to Bowen 374A resonance line */
	pbowen.ip374 = ipLinSafe((1.92),"He 2",0);
	pbowen.ip660 = ipLinSafe((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, from block data scalar */
		ipPredCont[i] = ipoint(EnrPredCont[i]);
	}

	/* for induced two photon emission we need mirror image set
	 * of pointers for continuum between Lya and helf Lya */
	ip2phtCom.ipHalfLya = ipoint(0.375);
	if( ip2phtCom.ipHalfLya > NIP2PHT )
	{
		fprintf( ioQQQ, " CreatePoint: too many continuum points needed to fill ip2pht array - increase NIP2PHT to at least%7ld\n", 
		  ip2phtCom.ipHalfLya );
		puts( "[Stop in CreatePoint]" );
		exit(1);
	}

	for( i=0; i < ip2phtCom.ipHalfLya; i++ )
	{
		/* energy is symmetric energy, the other side of half lya */
		energy = 0.75 - rfield.anu[i];
		ip2phtCom.ip2pht[i] = ipoint(energy);
	}

	/* do HeI */
	for( i=0; i < (NHE1LVL - 1); i++ )
	{
		nhe1Com.nhe1[i] = ipConSafe(He1NIonRyd.He1IonRyd[i],"He 1" );
		for( j=i + 1; j < NHE1LVL; j++ )
		{
			/* get energies from table, in zero */
			iphe1lCom.iphe1l[i][j] = ipLinSafe((He1NIonRyd.He1IonRyd[i]-
			  He1NIonRyd.He1IonRyd[j]),"He 1" , nhe1Com.nhe1[i]);
			/* make sure pointer not in next higher continuum */
			/*if( iphe1lCom.iphe1l[i][j] >= nhe1Com.nhe1[i] )
			{
				fprintf( ioQQQ, " He1 line%3ld%3ld spills over into%3ld continuum, energy=%11.3e\n", 
				  j, i, i, rfield.anu[iphe1lCom.iphe1l[i][j]-1] );
			}*/
		}
	}

	nhe1Com.nhe1[NHE1LVL-1] = ipConSafe(He1NIonRyd.He1IonRyd[NHE1LVL-1],
	  "He 1");

	nhe1Com.nhe1[NHE1LVL] = nhe1Com.nhe1[1];

	/* pointer to helium 2 trip s to ground decay */
	ip626.ipHe23sGrnd = ipLinSafe(1.457,"HeI3",0);

	/* photoionization for HeI triplet 23s */
	he.nhei3 = ipConSafe(0.3506,"HeI3");

	/* 666 error! remove following two, replace with above */
	he.nhei1 = nhe1Com.nhe1[0];
	he.nheii = nh.ipHn[1][IP1S];

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

	/* double photoionization
	 * 666 error! most of these have been removed */
	ip2gam.ip2ePhoto = ipoint(7.3);

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

	/* pointer to CaII XYZ at 8700A
	 * ipxyz = ipoint(0.1)
	 * */
	hehpCom.iheh1 = ipoint(1.6);
	hehpCom.iheh2 = ipoint(2.3);

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

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

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

	/* pointers to bound compton electron scattering of H, He */
	ircoil.ipRecoilIon = ipoint(194.);
	ircoil.irche = ipoint(260.);

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

	/* pointers to hot and cold iron */
	fekems.ipkcld = ipLinSafe(6400./13.54,"Fekc",0);
	fekems.ipkhot = ipLinSafe(6600./13.54,"Fekh",0);
	he3lines.ipHe3l[0] = ipLinSafe(0.08416,"HeI3",0);

	/* make low energy edges of some cells exact */
	fiddle(nh.ipHn[0][IP1S],0.99946);
	fiddle(nh.ipHn[0][IP2P],0.24986);

	fiddle(nhe1Com.nhe1[0],1.8072);
	fiddle(nh.ipHn[1][IP1S],4.00000);

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

	/* now save current form of array */
	for( i=0; i < NCELL; 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
	 * ipConSafe, 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 = 0;
	ipZ = 0;

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

	/*pointer to ionization threshold in energy array*/
	Heavy.ipHeavy[0][0] = nh.ipHn[nelem][IP1S];
	OpacPoint.ipElement[0][0][ipZ][nelem] = nh.ipHn[nelem][IP1S];

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

	if( elmton.lgElmtOn[1] )
	{
		/* set up shell pointers for helium */
		nelem = 1;
		ipZ = 0;

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

		/*pointer to ionization threshold in energy array*/
		Heavy.ipHeavy[0][1] = nhe1Com.nhe1[0];
		OpacPoint.ipElement[0][0][ipZ][nelem] = nhe1Com.nhe1[0];

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

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

		/*pointer to ionization threshold in energy array*/
		Heavy.ipHeavy[1][1] = nh.ipHn[nelem][IP1S];
		OpacPoint.ipElement[0][0][ipZ][nelem] = nh.ipHn[nelem][IP1S];

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

	/* begin sanity check
	 * check that ionization potential of neutral carbon valence shell is
	 * positive - check that block data has been linked in */
	if( PH1COM.PH1[2][5][5][0] <= 0. )
	{
		fprintf( ioQQQ, " PH1 returned non-positive threshold for Co.\n" );
		fprintf( ioQQQ, " Was the photo common block linked in?\n" );
		puts( "[Stop in CreatePoint]" );
		exit(1);
	}
	/*end sanity check*/

	/* fill in all shells for elements heavier than He */
	for( i=2; i<LIMELM; ++i )
	{
		/* i is the atomic number on the c scale, 2 for Li */
		ipShells(i);
	}

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

	/* finish off excited state photoionization
	 * do not use ipConSafe since only upper limit */
	OpacPoint.ipo3exc[1] = ipoint(5.0);

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

	/* following are various pointers for OI - Lybeta pump problem
	 * these are delta energies for boltzmann factors
	 * 666 error! 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 */
	ipoi.ipoiex[0] = ipoint(0.7005);
	ipoi.ipoiex[1] = ipoint(0.10791);
	ipoi.ipoiex[2] = ipoint(0.06925);
	ipoi.ipoiex[3] = ipoint(0.01151);
	ipoi.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 */
	OpacPoint.in1[0] = ipConSafe(0.893,"N1ex");

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

	if( (trace.lgTrace && trace.lgConBug) || (trace.lgTrace && trace.lgPointBug) )
	{
		fprintf( ioQQQ, "   CreatePoint:%5ldcells possible,%5ld used. N(1R):%4ld N(1.8):%4ld  N(4Ryd):%4ld N(O3)%4ld  N(x-ray):%5ld N(rcoil)%5ld\n", 
		  NCELL, rfield.nupper, nh.ipHn[0][IP1S], nhe1Com.nhe1[0], 
		  nh.ipHn[1][IP1S], OpacPoint.ipo3exc[0], nxray.ipCKshell, 
		  ircoil.ipRecoilIon );

		fprintf( ioQQQ, "   CreatePoint: HeI trip: %5ld\n", he.nhei3 );

		fprintf( ioQQQ, "   CreatePoint: ipEnerGammaRay: %5ld IPPRpari produc%5ld\n", 
		  egray.ipEnerGammaRay, OpacPoint.ippr );

		fprintf( ioQQQ, "   CreatePoint: H pointers;" );
		for( i=0; i <= 6; i++ )
		{
			fprintf( ioQQQ, "%4ld%4ld", i, nh.ipHn[0][i] );
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "   CreatePoint: HeI pnters;" );
		for( i=1; i <= 6; i++ )
		{
			fprintf( ioQQQ, "%4ld%4ld", i, nhe1Com.nhe1[i-1] );
		}
		fprintf( ioQQQ, "\n" );
		fprintf( ioQQQ, "   CreatePoint: Oxy pnters;" );

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

	/* Magnesium
	 * following is energy for phot of MG+ from exc state producing 2798 */
	OpacPoint.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) */
	OpacPoint.ica2ex[0] = ipConSafe(0.72,"Ca2x");
	OpacPoint.ica2ex[1] = ipoint(1.);

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

	/* center of V filter */
	ipvfil.ipVFilter = ipoint(0.164);

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

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

	/* this is 341.1 f=0.012 */
	pumpfe.pfe11b = (float)(2.75e-19/rfield.widflx[pumpfe.ipfe10-1]);
	pumpfe.pfe14 = (float)(1.15e-18/rfield.widflx[pumpfe.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++ )
	{
		TauLines[i].EnergyRyd = 
			(float)(TauLines[i].EnergyWN* WAVNRYD);

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

		TauLines[i].ipCont = 
			ipLinSafe(TauLines[i].EnergyRyd, chLab ,0);
		/* for debugging pointer - recall this number is cell+1
		   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 SetLine\n" );
			puts( "[Stop in setline]" );
			exit(1);
		}

		TauLines[i].damprel = 
			(float)(TauLines[i].Aul*
		  TauLines[i].WLAng*1e-8/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));

		/* 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 );
	}

	/* this is rest of hydrogenic species, after level 1 but before level 2,
	 * so that line labels remain meaningful, to go back to first time this
	 * block occurs, searh on ipZ=0 */
	for( ipZ=2; ipZ < LIMELM; ipZ++ )
	{
		if( elmton.lgElmtOn[ipZ])
		{
			/* generate label for this ion */
			sprintf( chLab, "%2s%2ld",elementnames.chElementSym[ipZ], ipZ+1 );
			/* Lya itself is the only transition below n=3 - we explicitly do not
			 * want to generate pointers for 2s-1s or 2p-2s */
			HydroLines[ipZ][IP2P][IP1S].ipCont = 
				ipLinSafe(HydroLines[ipZ][IP2P][IP1S].EnergyRyd , chLab ,0);
			/* set large negative pointers for 2p-2s and 2s-1s so we blow
			 * if we try to address 2s-2p or 1s-2s */
			HydroLines[ipZ][IP2S][IP1S].ipCont = INT_MIN;
			HydroLines[ipZ][IP2P][IP2S].ipCont = 1;
			for( ipLo=IP1S; ipLo <= nhlevel; ipLo++ )
			{
				/* HLevNIonRyd is level energy ryd 
				 * this is pointer to continuum ionization threshold 
				 * also puts label into continuum array */
				nh.ipHn[ipZ][ipLo] = ipConSafe(hydro.HLevNIonRyd[ipZ][ipLo],chLab);

				/* define all hydrogenic line pointers */
				for( ipHi=MAX2(3,ipLo+1); ipHi <= nhlevel; ipHi++ )
				{
					HydroLines[ipZ][ipHi][ipLo].ipCont = 
						ipLinSafe(HydroLines[ipZ][ipHi][ipLo].EnergyRyd , chLab , 
						nh.ipHn[ipZ][ipLo]);

					/* do not let pointer go into next higher continuum */
					/*HydroLines[ipZ][ipHi][ipLo].ipCont = 
						MIN2( j , (nh.ipHn[ipZ][ipLo]-1) );*/
				}
			}
			/* set pointer to two photon continuum as emission line to 
			 * completely bad value - if ever used should cause code to crash */
			HydroLines[ipZ][IP2S][IP1S].ipCont = INT_MIN;
		}
	}

	/* level 2 heavy element lines */
	for( i=0; i < nWindLine; i++ )
	{
		/* energy in ryd scale from wavenumber */
		TauLine2[i].EnergyRyd = (float)(TauLine2[i].EnergyWN* WAVNRYD);

		/* 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].damprel = 
			(float)(TauLine2[i].Aul*
		  TauLine2[i].WLAng*1e-8/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 = ipLinSafe(TauLine2[i].EnergyRyd , chLab,0 );
		/*if( TauLine2[i].ipCont==751 )
			fprintf(ioQQQ,"( 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();
	}

	/* 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]) );
			iWavLen = (long)TauLines[i].WLAng;
			if( iWavLen > 1000000 )
			{
				iWavLen /= 10000;
			}
			else if( iWavLen > 10000 )
			{
				iWavLen /= 1000;
			}

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

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

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

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

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

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