/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * For conditions of distribution and use see copyright notice in license.txt */
/*FeII_Colden maintain H2 column densities within X */
/*FeIILevelPops  main feii routine, called by CoolIron to evaluate iron cooling */
/*FeIICreate read in needed data from file 
 * convert form of feii data, as read in from file within routine FeIICreate
 * into physical form.  called by atmdat_readin  */
/*FeIIPunPop - punch level populations */
/*AssertFeIIDep called by assert feii depart coef command */
/*FeIIPrint print feii information */
/*FeIICollRatesBoltzmann evaluate collision strenths, 
 * both interpolating on r-mat and creating g-bar 
 * find Boltzmann factors, evaluate collisional rate coefficients */
/*FeIIPrint print output from large feii atom, called by prtzone */
/*FeIISumBand sum up large FeII emission over certain bands, called in lineset4 */
/*FeIITauInc called once per zone in RT_tau_inc to increment large FeII atom line optical depths */
/*FeII_RT_tau_reset reset optical depths for large FeII atom, called by update after each iteration */
/*FeIIPoint called by ContCreatePointers to create pointers for lines in large FeII atom */
/*FeIIAccel called by rt_line_driving to compute radiative acceleration due to FeII lines */
/*FeIIAddLines save accumulated FeII intensities called by lineset4 */
/*FeIIPunchLines punch feii lines at end of calculation, if punch verner set, called by dopunch*/
/*FeIIPunchLevels punch feii levels and energies */
/*FeII_LineZero zero out storage for large FeII atom, called by tauout */
/*FeIIIntenZero zero out intensity of FeII atom */
/*FeIIZero initialize some variables, called by zero one time before commands parsed */
/*FeIIReset reset some variables, called by zero */
/*FeIIPunData punch line data */ 
/*FeIIPunDepart punch some departure coef for large atom, 
 * set with punch feii departure command*/
/*FeIIPun1Depart send the departure coef for physical level nPUN to unit ioPUN */
/*FeIIContCreate create FeII continuum bins to add lines into ncell cells between wavelengths lambda low and high,
 * returns number of cells used */
/*FeIIBandsCreate returns number of feii bands */
/*FeII_RTOut do outward rates for FeII, called by RT_diffuse */
/*FeII_OTS do ots rates for FeII, called by RT_OTS */
/*FeIIRTMake called by RT_line_all, does large FeII atom radiative transfer */
/*FeIILyaPump find rate of Lya excitation of the FeII atom */
/*FeIIOvrLap handle overlapping FeII lines */
/*ParseAtomFeII parse the atom FeII command */
/*FeIIPunchLineStuff include FeII lines in punched optical depths, etc, called from PunchLineStuff */
#include "cddefines.h"
#include "cddrive.h"
#include "thermal.h"
#include "physconst.h"
#include "doppvel.h"
#include "taulines.h"
#include "dense.h"
#include "rfield.h"
#include "radius.h"
#include "lines_service.h"
#include "ipoint.h"
#include "lapack.h"
#include "hydrogenic.h"
#include "lines.h"
#include "rt.h"
#include "path.h"
#include "trace.h"
#include "punch.h"
#include "phycon.h"
#include "atomfeii.h"

/* FeIIOvrLap handle overlapping FeII lines */
static void FeIIOvrLap(void);

/* add FeII lines into ncell cells between wavelengths lambda low and high */
static void FeIIContCreate(double xLamLow , double xLamHigh , long int ncell );

/* read in the FeII Bands file, and sets nFeIIBands, the number of bands,
 * if argument is "" then use default file with bands, 
 * if filename within "" then use it instead,
 * return value is 0 if success, 1 if failed */
static int FeIIBandsCreate( const char chFile[] );

/* this will be the smallest collision strength we will permit with the gbar
 * collision strengths, or for the data that have zeroes */
/* >>chng 99 jul 24, this was 1e-9 in the old fortran version and in baldwin et al. 96,
 * and has been changed several times, and affects results.  this is the smallest
 * non-zero collision strength in the r-matrix calculations */
#define CS2SMALL 1e-5f
/* routines used only within this file */

/*FeIICollRatesBoltzmann evaluate collision strenths, 
 * both interpolating on r-mat and creating g-bar 
 * find Boltzmann factors, evaluate collisional rate coefficients */
static void FeIICollRatesBoltzmann(void);
/* find rate of Lya excitation of the FeII atom */
static void FeIILyaPump(void);

/*extern float Fe2LevN[NFE2LEVN][NFE2LEVN][NTA];*/
/*extern float Fe2LevN[ipHi][ipLo].t[NTA];*/
/*float ***Fe2LevN;*/
EmLine **Fe2LevN;
/* flag for the collision strength */
int **ncs1;

/* all following variables have file scope
#define	NFEIILINES	68635 */

/* this is size of nnPradDat array */
#define NPRADDAT 159

/* band wavelength, lower and upper bounds, in vacuum Angstroms */
/* FeII_Bands[n][3], where n is the number of bands in fe2Bands.dat
 * these bands are defined in fe2Bands.dat and read in at startup
 * of calculation */
float **FeII_Bands; 

/* continuum wavelengths, lower and upper bounds, in vacuum Angstroms,
 * third is integrated intensity */
/* FeII_Cont[n][3], where n is the number of cells needed
 * these bands are defined in FeIIContCreate */
float **FeII_Cont; 

/* this is the number of bands read in from bands_Fe2.dat */
long int nFeIIBands;

/* number of bands in continuum array */
long int nFeIIConBins;

/* the dim of this vector this needs one extra since there are 20 numbers per line,
 * highest not ever used for anything */
/*long int nnPradDat[NPRADDAT+1];*/
static long int *nnPradDat;

/* malloced in feiidata */
/* float sPradDat[NPRADDAT][NPRADDAT][8];*/
/* float sPradDat[ipHi][ipLo].t[8];*/
static float ***sPradDat;

/* array used to integrate feii line intensities over model,
 * Fe2SavN[upper][lower],
 *static float Fe2SavN[NFE2LEVN][NFE2LEVN];*/
static double **Fe2SavN;

/* save effective transition rates */
static double **Fe2A;

/* induced transition rates */
static double **Fe2LPump, **Fe2CPump;

/* collision rates */
static float **Fe2Coll;

/* Fe2DepCoef[NFE2LEVN];
float cli[NFEIILINES], 
  cfe[NFEIILINES];*/
static double 
	/* departure coefficients */
	*Fe2DepCoef , 
	/* level populations */
	*Fe2LevelPop ,
	/* column densities */
	*Fe2ColDen ,
	/* this will become array of Boltzmann factors */
	*FeIIBoltzmann ;
	/*FeIIBoltzmann[NFE2LEVN] ,*/

static double EnerLyaProf1, 
  EnerLyaProf4, 
  PhotOccNumLyaCenter;
static double 
		/* the yVector - will become level populations after matrix inversion */
		*yVector,
	  /* this is used to call matrix routines */
	  /*xMatrix[NFE2LEVN][NFE2LEVN] ,*/
	  **xMatrix , 
	  /* this will become the very large 1-D array that 
	   * is passed to the matrix inversion routine*/
	  *amat;


/*FeII_Colden maintain H2 column densities within X */
void FeII_Colden( const char *chLabel )
{
	long int n;

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

	/* >>chng 05 nov 29, FeII always on, always want to evaluate this */
	/* nothing to do if no FeII 
	if( !FeII. lgFeIION )
		return;*/

	if( strcmp(chLabel,"ZERO") == 0 )
	{
		/* zero out column density */
		for( n=0; n < FeII.nFeIILevel; ++n )
		{
			/* space for the rotation quantum number */
			Fe2ColDen[n] = 0.;
		}
	}

	else if( strcmp(chLabel,"ADD ") == 0 )
	{
		/*  add together column densities */
		for( n=0; n < FeII.nFeIILevel; ++n )
		{
			/* state-specific FeII column density */
			Fe2ColDen[n] += Fe2LevelPop[n]*radius.drad_x_fillfac;
		}
	}

	/* check for the print option, which we can't do, else  we have a problem */
	else if( strcmp(chLabel,"PRIN") != 0 )
	{
		fprintf( ioQQQ, " FeII_Colden does not understand the label %s\n", 
		  chLabel );
		puts( "[Stop in FeII_Colden]" );
		cdEXIT(EXIT_FAILURE);
	}

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

	return;
}

/*
 *=====================================================================
 */
/* FeIICreate read in feii data from files on disk.  called by atmdat_readin 
 * but only if FeII. lgFeIION is true, set with atom feii verner command */
void FeIICreate(void)
{
	FILE *ioDATA;

	char chLine[FILENAME_PATH_LENGTH_2] , 
		chFilename[FILENAME_PATH_LENGTH_2] ;

	long int i, 
	  ipHi ,
	  ipLo,
	  lo,
	  ihi,
	  k, 
	  m1, 
	  m2, 
	  m3;

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

	if( lgFeIIMalloc )
	{
		/* we have already been called one time, just bail out */
#		ifdef DEBUG_FUN
		fputs( " <->FeIICreate()\n", debug_fp );
#		endif

		return;
	}

	/* now set flag so never done again - this is set FALSE in cdDefines
	 * when this is true it is no longer possible to change the number of levels
	 * in the model atom with the atom feii levels command */
	lgFeIIMalloc = TRUE;

	/* remember how many levels this was, so that in future calculations
	 * we can reset the atom to full value */
	FeII.nFeIILevelAlloc = FeII.nFeIILevel;

	/* set up array to save FeII transition probabilities */
	if( (Fe2A = (double **)MALLOC(sizeof(double *)*(unsigned long)FeII.nFeIILevel )) ==NULL)
	{
		fprintf( ioQQQ, 
			" FeIICreate could not MALLOC Fe2A 1\n" );
		puts( "[Stop in FeIICreate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* second dimension, lower level, for line save array */
	for( ipHi=0; ipHi < FeII.nFeIILevel; ++ipHi )
	{
		if( (Fe2A[ipHi]=(double*)MALLOC(sizeof(double )*(unsigned long)FeII.nFeIILevel)) ==NULL)
		{
			fprintf( ioQQQ, 
				" FeIICreate could not MALLOC Fe2A 2\n" );
			puts( "[Stop in FeIICreate]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* set up array to save FeII pumping rates */
	if( (Fe2CPump = (double **)MALLOC(sizeof(double *)*(unsigned long)FeII.nFeIILevel )) ==NULL)
	{
		fprintf( ioQQQ, 
			" FeIICreate could not MALLOC Fe2CPump 1\n" );
		puts( "[Stop in FeIICreate]" );
		cdEXIT(EXIT_FAILURE);
	}
	/* set up array to save FeII pumping rates */
	if( (Fe2LPump = (double **)MALLOC(sizeof(double *)*(unsigned long)FeII.nFeIILevel )) ==NULL)
	{
		fprintf( ioQQQ, 
			" FeIICreate could not MALLOC Fe2LPump 1\n" );
		puts( "[Stop in FeIICreate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* second dimension, lower level, for line save array */
	for( ipHi=0; ipHi < FeII.nFeIILevel; ++ipHi )
	{
		if( (Fe2CPump[ipHi]=(double*)MALLOC(sizeof(double )*(unsigned long)FeII.nFeIILevel)) ==NULL)
		{
			fprintf( ioQQQ, 
				" FeIICreate could not MALLOC Fe2CPump 2\n" );
			puts( "[Stop in FeIICreate]" );
			cdEXIT(EXIT_FAILURE);
		}
		if( (Fe2LPump[ipHi]=(double*)MALLOC(sizeof(double )*(unsigned long)FeII.nFeIILevel)) ==NULL)
		{
			fprintf( ioQQQ, 
				" FeIICreate could not MALLOC Fe2LPump 2\n" );
			puts( "[Stop in FeIICreate]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* set up array to save FeII collision rates */
	if( (Fe2Coll = (float **)MALLOC(sizeof(float *)*(unsigned long)FeII.nFeIILevel )) ==NULL)
	{
		fprintf( ioQQQ, 
			" FeIICreate could not MALLOC Fe2Coll 1\n" );
		puts( "[Stop in FeIICreate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* second dimension, lower level, for line save array */
	for( ipHi=0; ipHi < FeII.nFeIILevel; ++ipHi )
	{
		if( (Fe2Coll[ipHi]=(float*)MALLOC(sizeof(float )*(unsigned long)FeII.nFeIILevel)) ==NULL)
		{
			fprintf( ioQQQ, 
				" FeIICreate could not MALLOC Fe2Coll 2\n" );
			puts( "[Stop in FeIICreate]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* MALLOC space for the 2D matrix array */

	if(  (xMatrix = (double **)MALLOC(sizeof(double *)*(unsigned long)FeII.nFeIILevel ))==NULL  )
	{
		fprintf( ioQQQ, " could not MALLOC xMatrix arrays in 1D\n" );
		puts( "[Stop in FeIICreate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* now do the second dimension */
	for( i=0; i<FeII.nFeIILevel; ++i )
	{
		if( (xMatrix[i] = (double *)MALLOC(sizeof(double)*(unsigned long)FeII.nFeIILevel ) )==NULL )
		{
			fprintf( ioQQQ, " could not MALLOC 2D xMatrix\n" );
			puts( "[Stop in FeIICreate]" );
			cdEXIT(EXIT_FAILURE);
		}
	}
	/* MALLOC space for the  1-yVector array */
	if( (amat=(double*)MALLOC( (sizeof(double)*(unsigned long)(FeII.nFeIILevel*FeII.nFeIILevel) ))) == NULL )
	{
		fprintf( ioQQQ, " FeIICreate MALLOC amat error\n" );
		puts( "[Stop in FeIICreate]" );
		cdEXIT(EXIT_FAILURE);
	}
	/* MALLOC space for the  1-yVector array */
	if( (yVector=(double*)MALLOC( (sizeof(double)*(unsigned long)(FeII.nFeIILevel) ))) == NULL )
	{
		fprintf( ioQQQ, " FeIICreate MALLOC yVector error\n" );
		puts( "[Stop in FeIICreate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* set up array to save FeII line intensities */
	if( (Fe2SavN = (double **)MALLOC(sizeof(double *)*(unsigned long)FeII.nFeIILevel )) ==NULL)
	{
		fprintf( ioQQQ, 
			" FeIICreate could not MALLOC Fe2SavN 1\n" );
		puts( "[Stop in FeIICreate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* second dimension, lower level, for line save array */
	for( ipHi=1; ipHi < FeII.nFeIILevel; ++ipHi )
	{
		if( (Fe2SavN[ipHi]=(double*)MALLOC(sizeof(double )*(unsigned long)ipHi)) ==NULL)
		{
			fprintf( ioQQQ, 
				" FeIICreate could not MALLOC Fe2SavN 2\n" );
			puts( "[Stop in FeIICreate]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* now MALLOC space for energy level table*/
	if( (nnPradDat = (long*)MALLOC( (NPRADDAT+1)*sizeof(long) )) == NULL )
	{
		fprintf(ioQQQ," FeIICreate could not MALLOC nnPradDat\n");
		puts( "[Stop in FeIICreate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/*Fe2DepCoef[NFE2LEVN];*/
	if( ( Fe2DepCoef = (double*)MALLOC( (unsigned long)FeII.nFeIILevel*sizeof(double) )) == NULL )
	{
		fprintf(ioQQQ," FeIICreate could not MALLOC Fe2DepCoef\n");
		puts( "[Stop in FeIICreate]" );
		cdEXIT(EXIT_FAILURE);
	}
	/*Fe2LevelPop[NFE2LEVN];*/
	if( ( Fe2LevelPop = (double*)MALLOC( (unsigned long)FeII.nFeIILevel*sizeof(double) )) == NULL )
	{
		fprintf(ioQQQ," FeIICreate could not MALLOC Fe2LevelPop\n");
		puts( "[Stop in FeIICreate]" );
		cdEXIT(EXIT_FAILURE);
	}
	/*Fe2ColDen[NFE2LEVN];*/
	if( ( Fe2ColDen = (double*)MALLOC( (unsigned long)FeII.nFeIILevel*sizeof(double) )) == NULL )
	{
		fprintf(ioQQQ," FeIICreate could not MALLOC Fe2ColDen\n");
		puts( "[Stop in FeIICreate]" );
		cdEXIT(EXIT_FAILURE);
	}
	/*FeIIBoltzmann[NFE2LEVN];*/
	if( ( FeIIBoltzmann = (double*)MALLOC( (unsigned long)FeII.nFeIILevel*sizeof(double) )) == NULL )
	{
		fprintf(ioQQQ," FeIICreate could not MALLOC FeIIBoltzmann\n");
		puts( "[Stop in FeIICreate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* MALLOC the float sPradDat[NPRADDAT][NPRADDAT][8] array */
	/* MALLOC the float sPradDat[ipHi][ipLo].t[8] array */
	if( (sPradDat = ((float ***)MALLOC(NPRADDAT*sizeof(float **)))) == NULL )
	{
		fprintf(ioQQQ," FeIICreate could not MALLOC 1 sPradDat\n");
		puts( "[Stop in FeIICreate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* >>chng 00 dec 06, changed lower limit of loop to 1, Tru64 alpha's will not allocate 0 bytes!, PvH */
	sPradDat[0] = NULL;
	for (ipHi=1; ipHi < NPRADDAT; ipHi++) 
	{
		/* >>chng 00 dec 06, changed sizeof(float) to sizeof(float*), PvH */
		if(  (sPradDat[ipHi] = (float **)MALLOC((unsigned long)ipHi*sizeof(float *))) == NULL )
		{
			fprintf(ioQQQ," FeIICreate could not MALLOC 2 sPradDat\n");
			puts( "[Stop in FeIICreate]" );
			cdEXIT(EXIT_FAILURE);
		}
		/* now make space for the second dimension */
		for( ipLo=0; ipLo< ipHi; ipLo++ )
		{ 
			if(  (sPradDat[ipHi][ipLo] = (float *)MALLOC(8*sizeof(float ))) == NULL )
			{
				fprintf(ioQQQ," FeIICreate could not MALLOC 3 sPradDat\n");
				puts( "[Stop in FeIICreate]" );
				cdEXIT(EXIT_FAILURE);
			}
		}
	}

	/* now set junk to make sure reset before used */
	for (ipHi=0; ipHi < NPRADDAT; ipHi++) 
	{
		for( ipLo=0; ipLo< ipHi; ipLo++ )
		{ 
			for( k=0; k<8; ++k )
			{
				sPradDat[ipHi][ipLo][k] = -FLT_MAX;
			}
		}
	}

	/* now create the main emission line array and a helper array for the cs flag  */
	if( (Fe2LevN=(EmLine**)MALLOC(sizeof(EmLine *)*(unsigned long)FeII.nFeIILevel)) ==NULL)
		BadMalloc();
	if( (ncs1=(int**)MALLOC(sizeof(int *)*(unsigned long)FeII.nFeIILevel)) ==NULL)
		BadMalloc();

	for( ipHi=1; ipHi < FeII.nFeIILevel; ++ipHi )
	{
		if( (Fe2LevN[ipHi]=(EmLine*)MALLOC(sizeof(EmLine)*(unsigned long)ipHi))==NULL)
			BadMalloc();
		if( (ncs1[ipHi]=(int*)MALLOC(sizeof(int)*(unsigned long)ipHi))==NULL)
			BadMalloc();
	}

	/* now that its created, set to junk */
	for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel; ipHi++ )
		{
			EmLineJunk( &Fe2LevN[ipHi][ipLo] );
		}
	}


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

	if( trace.lgTrace )
	{
		fprintf( ioQQQ," FeIICreate opening fe2nn.dat:");
	}

	if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
	{
		fprintf( ioQQQ, " FeIICreate could not open fe2nn.dat\n" );
		if( lgDataPathSet == TRUE )
		{
			fprintf( ioQQQ, " FeIICreate could not open fe2nn.dat, even tried path.\n");
			fprintf( ioQQQ, " path is *%s*\n",chDataPath );
			fprintf( ioQQQ, " final path is *%s*\n",chFilename );
		}
		puts( "[Stop in FeIICreate]" );
		cdEXIT(EXIT_FAILURE);
	}

	ASSERT( ioDATA !=NULL );
	for( i=0; i < 8; i++ )
	{
		if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
		{
			fprintf( ioQQQ, " fe2nn.dat error reading data\n" );
			puts( "[Stop in FeIICreate]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* get these integers from fe2nn.dat */
		k = 20*i;
		/* NPRADDAT is size of nnPradDat array, 159+1, make sure we do not exceed it */
		ASSERT( k+19 < NPRADDAT+1 );
		sscanf( chLine ,
			"%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld",
			&nnPradDat[k+0], &nnPradDat[k+1],  &nnPradDat[k+2], &nnPradDat[k+3], &nnPradDat[k+4],
			&nnPradDat[k+5], &nnPradDat[k+6],  &nnPradDat[k+7], &nnPradDat[k+8], &nnPradDat[k+9],
			&nnPradDat[k+10],&nnPradDat[k+11], &nnPradDat[k+12],&nnPradDat[k+13],&nnPradDat[k+14],
			&nnPradDat[k+15],&nnPradDat[k+16], &nnPradDat[k+17],&nnPradDat[k+18],&nnPradDat[k+19]
			);
#		if !defined(NDEBUG)
		for( m1=0; m1<20; ++m1 )
		{
			ASSERT( nnPradDat[k+m1] >= 0 && nnPradDat[k+m1] <= NFE2LEVN );
		}
#		endif
	}
	fclose(ioDATA);

	/* now get radiation data.
	 * These are from the following data sources:
	 >>refer	fe2	as	Nahar, S., 1995, A&A 293, 967
	 >>refer	fe2	as	Quinet, P., LeDourneuf, M., & Zeipppen C.J., 1996, A&AS, 120, 361 
	 >>refer	fe2	as	Furh, J.R., Martin, G.A., & Wiese, W.L., 1988; J Phys Chem Ref Data 17, Suppl 4
	 >>refer	fe2	as	Giridhar, S., & Arellano Ferro, A., 1995; Ref Mexicana Astron Astrofis 31, 23
	 >>refer	fe2	as	Kurucz, R.L., 1995, SAO CD ROM 23
	 */
	if( lgDataPathSet == TRUE )
	{
		/*path set, so look only there */
		strcpy( chFilename , chDataPath );
		strcat( chFilename , "fe2rad.dat" );
	}
	else
	{
		/* path not set, check local space only */
		strcpy( chFilename , "fe2rad.dat" );
	}

	if( trace.lgTrace )
	{
		fprintf( ioQQQ," FeIICreate opening fe2rad.dat:");
	}

	if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
	{
		fprintf( ioQQQ, " FeIICreate could not open fe2rad.dat\n" );
		if( lgDataPathSet == TRUE )
		{
			fprintf( ioQQQ, " FeIICreate could not open fe2rad.dat, even tried path.\n");
			fprintf( ioQQQ, " path is *%s*\n",chDataPath );
			fprintf( ioQQQ, " final path is *%s*\n",chFilename );
		}
		puts( "[Stop in FeIICreate]" );
		cdEXIT(EXIT_FAILURE);
	}
	ASSERT( ioDATA !=NULL );

	/* get the first line, this is a version number */
	if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
	{
		fprintf( ioQQQ, " fe2rad.dat error reading data\n" );
		puts( "[Stop in FeIICreate]" );
		cdEXIT(EXIT_FAILURE);
	}
	/* scan off three integers */
	sscanf( chLine ,"%ld%ld%ld",&lo, &ihi, &m1);
	if( lo!=3 || ihi!=1 || m1!=28 )
	{
		fprintf( ioQQQ, " fe2rad.dat has the wrong magic numbers, expected 3 1 28\n" );
		puts( "[Stop in FeIICreate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* file now open, read the data */
	for( ipHi=1; ipHi < FeII.nFeIILevel; ipHi++ )
	{
		for( ipLo=0; ipLo < ipHi; ipLo++ )
		{
			float save[2];
			/* keep reading until have non-comment line, one that does not
			 * start with # */
			chLine[0] = '#';
			while( chLine[0] == '#' )
			{
				if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
				{
					fprintf( ioQQQ, " fe2nn.dat error reading data\n" );
					puts( "[Stop in FeIICreate]" );
					cdEXIT(EXIT_FAILURE);
				}
			}

			/* first and last number on this line */
			sscanf( chLine ,
				"%ld%ld%ld%ld%f%f%ld",
				&lo, &ihi, &m1, &m2 ,
				&save[0], &save[1] , &m3);
			/* the pointeres ihi and lo within this array were 
			 * read in from the line above */
			Fe2LevN[ihi-1][lo-1].gLo = (float)m1;
			Fe2LevN[ihi-1][lo-1].gHi = (float)m2;
			Fe2LevN[ihi-1][lo-1].Aul = save[0];
			Fe2LevN[ihi-1][lo-1].EnergyWN = save[1];
			if( Fe2LevN[ihi-1][lo-1].Aul == 1e3f )
			{
				/*float xmicron;
				xmicron = 1e4f/ Fe2LevN[ihi-1][lo-1].EnergyWN;*/
				/* these are made-up intercombination lines, set gf to 1e-5
				 * then get better A */
				Fe2LevN[ihi-1][lo-1].gf = 1e-5f;

				Fe2LevN[ihi-1][lo-1].Aul = Fe2LevN[ihi-1][lo-1].gf*(float)TRANS_PROB_CONST*
					POW2(Fe2LevN[ihi-1][lo-1].EnergyWN)*Fe2LevN[ihi-1][lo-1].gLo/Fe2LevN[ihi-1][lo-1].gHi;

				/*fprintf(ioQQQ," %e from 1e3 to %e\n",xmicron , Fe2LevN[ihi-1][lo-1].Aul );*/
			}
			/* this is the last column of fe2rad.dat, and is 1, 2, or 3.  
			 * 1 signifies that transition is permitted,
			 * 2 is semi-forbidden
			 * 3 forbidden, within lowest 63 levels are forbidden, first permitted
			 * transition is from 64 */
			ncs1[ihi-1][lo-1] = (int)m3;
		}
	}
	fclose(ioDATA);

	/* now read collision data in fe2col.dat 
	 * These are from the following sources
	 >>refer	fe2	cs	Zhang, H.L., & Pradhan, A.K., 1995, A&A 293, 953 
	 >>refer	fe2	cs	Bautista, M., (private communication), 
	 >>refer	fe2	cs	Mewe, R., 1972, A&AS 20, 215 (the g-bar approximation)
	 */
	if( lgDataPathSet == TRUE )
	{
		/*path set, so look only there */
		strcpy( chFilename , chDataPath );
		strcat( chFilename , "fe2col.dat" );
	}
	else
	{
		/* path not set, check local space only */
		strcpy( chFilename , "fe2col.dat" );
	}

	if( trace.lgTrace )
	{
		fprintf( ioQQQ," FeIICreate opening fe2col.dat:");
	}

	if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
	{
		fprintf( ioQQQ, " FeIICreate could not open fe2col.dat\n" );
		if( lgDataPathSet == TRUE )
		{
			fprintf( ioQQQ, " FeIICreate could not open fe2col.dat, even tried path.\n");
			fprintf( ioQQQ, " path is *%s*\n",chDataPath );
			fprintf( ioQQQ, " final path is *%s*\n",chFilename );
		}
		puts( "[Stop in FeIICreate]" );
		cdEXIT(EXIT_FAILURE);
	}

	ASSERT( ioDATA !=NULL);
	for( ipHi=1; ipHi<NPRADDAT; ++ipHi )
	{
		for( ipLo=0; ipLo<ipHi; ++ipLo )
		{
			float save[8];
 			if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
			{
				fprintf( ioQQQ, " fe2col.dat error reading data\n" );
				puts( "[Stop in FeIICreate]" );
				cdEXIT(EXIT_FAILURE);
			}
			sscanf( chLine ,
				"%ld%ld%f%f%f%f%f%f%f%f",
				&m1, &m2,
				&save[0], &save[1] , &save[2],&save[3], &save[4] , &save[5],
				&save[6], &save[7] 
				);
			for( k=0; k<8; ++k )
			{
				/* the max is here because there are some zeroes in the data file.
				 * this is unphysical but is part of their distribution.  As a result
				 * don't let the cs fall below 0.01 */
				sPradDat[m2-1][m1-1][k] = MAX2(CS2SMALL , save[k] );
			}
		}
	}

	/*generate needed opacity data for the large FeII atom */

	/* this routine is called only one time when cloudy is init
	 * for the very first time.  It converts the FeII data stored
	 * in the large feii arrays into the array storage needed by cloudy
	 * for its line optical depth arrays
	 */

	/* convert large feii line arrays into standard heavy el ar */
	for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel; ipHi++ )
		{
			/* pull information out of existing data arrays */

			ASSERT( Fe2LevN[ipHi][ipLo].EnergyWN > 0. );
			ASSERT( Fe2LevN[ipHi][ipLo].Aul> 0. );

			/* now put into standard internal line format
			Fe2LevN[ipHi][ipLo].WLAng = (float)(1.e8/Fe2LevN[ipHi][ipLo].EnergyWN); */
			/* >>chng 03 oct 28, above neglected index of refraction of air -
			 * convert to below */
			Fe2LevN[ipHi][ipLo].WLAng = 
				(float)(1.0e8/
			Fe2LevN[ipHi][ipLo].EnergyWN/
			RefIndex( Fe2LevN[ipHi][ipLo].EnergyWN ));

			/* generate gf from A */
			Fe2LevN[ipHi][ipLo].gf = 
				(float)(Fe2LevN[ipHi][ipLo].Aul*Fe2LevN[ipHi][ipLo].gHi/
				TRANS_PROB_CONST/POW2(Fe2LevN[ipHi][ipLo].EnergyWN));

			Fe2LevN[ipHi][ipLo].IonStg = 2;
			Fe2LevN[ipHi][ipLo].nelem = 26;
			/* which redistribution function??  
			 * For resonance line use K2 (-1),
			 * for subordinate lines use CRD with wings */
			/* >>chng 01 mar 09, all had been 1, inc redis with wings */
			/* to reset this, so that the code works as it did pre 01 mar 29,
			 * use command 
			 * atom feii redistribution resonance pdr 
			 * atom feii redistribution subordinate pdr */
			if( ipLo == 0 )
			{
				Fe2LevN[ipHi][ipLo].iRedisFun = FeII.ipRedisFcnResonance ;
			}
			else
			{
				/* >>chng 01 feb 27, had been -1, crd with core only,
				 * change to crd with wings as per discussion with Ivan Hubeny */
				Fe2LevN[ipHi][ipLo].iRedisFun = FeII.ipRedisFcnSubordinate ;
			}
			Fe2LevN[ipHi][ipLo].phots = 0.;
			Fe2LevN[ipHi][ipLo].FracInwd = 1.;
		}
	}

	/* finally get FeII bands, this sets  */
	if( FeIIBandsCreate("") )
	{
		fprintf( ioQQQ," failed to open FeII bands file \n");
		puts( "[Stop in FeIICreate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* now establish the FeII continuum, these are set to
	 * 1000, 7000, and 1000 in FeIIZero in this file, and
	 * are reset with the atom feii continuum command */
	FeIIContCreate( FeII.fe2con_wl1 , FeII.fe2con_wl2 , FeII.nfe2con );

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

	return;
}

/*
 *=====================================================================
 */
/***********************************************************************
 *** version of FeIILevelPops.f with overlap in procces 05.28.97 ooooooo
 ******change in common block *te* sqrte 05.28.97
 *******texc is fixed 03.28.97
 *********this version has work on overlap 
 *********this version has # of zones (ZoneCnt) 07.03.97
 *********taux - optical depth depends on iter correction 03.03.97
 *********calculate cooling (Fe2_large_cool) and output fecool from Cloudy 01.29.97god
 *********and cooling derivative (ddT_Fe2_large_cool)
 ************ this program for 371 level iron model 12/14/1995
 ************ this program for 371 level iron model 1/11/1996
 ************ this program for 371 level iron model 3/21/1996
 ************ this program without La 3/27/1996
 ************ this program for 371 level iron model 4/9/1996
 ************ includes:FeIICollRatesBoltzmann;cooling;overlapping in lines */
void FeIILevelPops( void )
{
	long int  i, 
		ipHi ,
		ipLo ,
		n;
	/* used in test for non-positive level populations */
	int lgPopNeg ;

	double 
	  EnergyWN,
	  pop ,
	  sum;

	int32 info;
	int32 ipiv[NFE2LEVN];

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

	if( trace.lgTrace )
	{
		fprintf( ioQQQ,"   FeIILevelPops fe2 pops called\n");
	}

	/* FeII.lgSimulate was set true with simulate flag on atom feii command,
	 * for bebugging without actually calling the routine */
	if( FeII.lgSimulate )
	{

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

		return;
	}

	/* zero out some arrays */
	for( n=0; n<FeII.nFeIILevel; ++n)
	{
		for( ipLo=0; ipLo<FeII.nFeIILevel; ++ipLo )
		{
			Fe2CPump[ipLo][n] = 0.;
			Fe2LPump[ipLo][n] = 0.;
			Fe2A[ipLo][n] = 0.;
			xMatrix[ipLo][n] = 0.;
		}
	}

	/* make the g-bar collision strengths and do linear interpolation on r-matrix data.
	 * this also sets Boltzmann factors for all levels,
	 * sets values of FeColl used below, but only if temp has changed */
	FeIICollRatesBoltzmann();

	/* pumping and spontantous decays */
	for( n=0; n<FeII.nFeIILevel; ++n)
	{
		for( ipHi=n+1; ipHi<FeII.nFeIILevel; ++ipHi )
		{
			/* continuum pumping rate from n to upper ipHi */
			Fe2CPump[n][ipHi] = Fe2LevN[ipHi][n].pump;

			/* continuum pumping rate from ipHi to lower n */
			Fe2CPump[ipHi][n] = Fe2LevN[ipHi][n].pump*
				Fe2LevN[ipHi][n].gHi/Fe2LevN[ipHi][n].gLo;

			/* spontaneous decays */
			Fe2A[ipHi][n] = Fe2LevN[ipHi][n].Aul*(Fe2LevN[ipHi][n].Pesc + Fe2LevN[ipHi][n].Pelec_esc +
			  Fe2LevN[ipHi][n].Pdest);
		}
	}

	/* now do pumping of atom by Lya */
	FeIILyaPump();

	/* ************************************************************************** 
	 *
	 * final rates into matrix 
	 *
	 ***************************************************************************/

	for( n=0; n<FeII.nFeIILevel; ++n)
	{
		/* all processes leaving level n going down*/
		for( ipLo=0; ipLo<n; ++ipLo )
		{
			xMatrix[n][n] = xMatrix[n][n] + Fe2CPump[n][ipLo] + Fe2LPump[n][ipLo]+ Fe2A[n][ipLo] + 
				Fe2Coll[n][ipLo]*dense.eden;
		}
		/* all processes leaving level n going up*/
		for( ipHi=n+1; ipHi<FeII.nFeIILevel; ++ipHi )
		{
			xMatrix[n][n] = xMatrix[n][n] + Fe2CPump[n][ipHi] + Fe2LPump[n][ipHi] + Fe2Coll[n][ipHi]*dense.eden;
		}
		/* all processes entering level n from below*/
		for( ipLo=0; ipLo<n; ++ipLo )
		{
			xMatrix[ipLo][n] = xMatrix[ipLo][n] - Fe2CPump[ipLo][n] - Fe2LPump[ipLo][n] - Fe2Coll[ipLo][n]*dense.eden;
		}
		/* all processes entering level n from above*/
		for( ipHi=n+1; ipHi<FeII.nFeIILevel; ++ipHi )
		{
			xMatrix[ipHi][n] = xMatrix[ipHi][n] - Fe2CPump[ipHi][n] - Fe2LPump[ipHi][n] - Fe2Coll[ipHi][n]*dense.eden - 
				Fe2A[ipHi][n];
		}

		/* the top row of the matrix is the sum of level populations */
		xMatrix[n][0] = 1.0;
	}

	{
		/* option to print out entire matrix */
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC )
		{
			/* print the matrices */
			for( n=0; n<FeII.nFeIILevel; ++n)
			{
				/*fprintf(ioQQQ,"\n");*/
				/* now print the matrix*/
				for( ipLo=0; ipLo<FeII.nFeIILevel; ++ipLo )
				{
					fprintf(ioQQQ," %.1e", xMatrix[n][ipLo]);
				}
				fprintf(ioQQQ,"\n");
			}
		}
	}

	/* define the Y Vector.  The oth element is the sum of all level populations
	 * adding up to the total population.  The remaining elements are the level
	 * balance equations adding up to zero */
	yVector[0] = 1.0;
	for( n=1; n < FeII.nFeIILevel; n++ )
	{
		yVector[n] = 0.0;
	}

	/* create the 1-yVector array that will save vector,
	 * this is the macro trick */
#	ifdef AMAT
#		undef AMAT
#	endif
#	define AMAT(I_,J_)	(*(amat+(I_)*FeII.nFeIILevel+(J_)))

	/* copy current contents of xMatrix array over to special amat array*/
	for( ipHi=0; ipHi < FeII.nFeIILevel; ipHi++ )
	{
		for( i=0; i < FeII.nFeIILevel; i++ )
		{
			/*amat[i][ipHi] = z[i][ipHi];*/
			AMAT(i,ipHi) = xMatrix[i][ipHi];
		}
	}

	info = 0;

	/* do the linear algebra to find the level populations */
	getrf_wrapper(FeII.nFeIILevel, FeII.nFeIILevel, amat, FeII.nFeIILevel, ipiv, &info);
	getrs_wrapper('N', FeII.nFeIILevel, 1, amat, FeII.nFeIILevel, ipiv, yVector, FeII.nFeIILevel, &info);

	if( info != 0 )
	{
		fprintf( ioQQQ, "DISASTER FeIILevelPops: dgetrs finds singular or ill-conditioned matrix\n" );
		puts( "[Stop in FeIILevelPops" );
		cdEXIT(EXIT_FAILURE);
	}

	/* yVector now contains the level populations */

	/* this better be false after this loop - if not then non-positive level pops */
	lgPopNeg = FALSE;
	/* copy all level pops over to Fe2LevelPop */
	for( ipLo=0; ipLo < FeII.nFeIILevel; ipLo++ )
	{
		if(yVector[ipLo] <= 0. )
		{
			lgPopNeg = TRUE;
			fprintf(ioQQQ,"PROBLEM FeIILevelPops finds non-positive level population, level is %ld pop is %g\n",
				ipLo , yVector[ipLo] );
		}
		/* this is now correct level population, cm^-3 */
		Fe2LevelPop[ipLo] = yVector[ipLo] * dense.xIonDense[ipIRON][1];
	}
	if( lgPopNeg )
	{
		fprintf(ioQQQ," PROBLEM FeIILevelPops finds non-positive level population\n" );
	}

	/* now set line opacities, intensities, and level populations*/
	for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ )
	{
		for( ipHi=ipLo+1; ipHi < FeII.nFeIILevel; ipHi++ )
		{
			Fe2LevN[ipHi][ipLo].PopOpc = ((yVector[ipLo] - 
				yVector[ipHi]*Fe2LevN[ipHi][ipLo].gLo/Fe2LevN[ipHi][ipLo].gHi)*
				dense.xIonDense[ipIRON][1]);

			Fe2LevN[ipHi][ipLo].PopLo = yVector[ipLo]*dense.xIonDense[ipIRON][1];

			Fe2LevN[ipHi][ipLo].PopHi = yVector[ipHi]*dense.xIonDense[ipIRON][1];

			Fe2LevN[ipHi][ipLo].phots = yVector[ipHi]*
			  Fe2LevN[ipHi][ipLo].Aul*(Fe2LevN[ipHi][ipLo].Pesc + Fe2LevN[ipHi][ipLo].Pelec_esc )*
			  dense.xIonDense[ipIRON][1];

			Fe2LevN[ipHi][ipLo].xIntensity = Fe2LevN[ipHi][ipLo].phots*
			  Fe2LevN[ipHi][ipLo].EnergyErg;

			/* ratio of collisional (new) to pumped excitations */
			/* >>chng 02 mar 04, add test MAX2 to prevent div by zero */
			Fe2LevN[ipHi][ipLo].ColOvTot = (float)(Fe2Coll[ipLo][ipHi]*dense.eden /
				SDIV( Fe2Coll[ipLo][ipHi]*dense.eden + Fe2CPump[ipLo][ipHi] + Fe2LPump[ipLo][ipHi] ) );
		}
	}

	/* only do this if large atom is on and more than ground terms computed - 
	 * do not if only lowest levels are computed */
	if( FeII.lgFeIILargeOn )
	{
		/* the hydrogen Lya destruction rate, then probability */
		hydro.dstfe2lya = 0.;
		EnergyWN = 0.;
		/* count how many photons were removed-added */
		for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ )
		{
			for( ipHi=ipLo+1; ipHi < FeII.nFeIILevel; ipHi++ )
			{
				EnergyWN += Fe2LPump[ipLo][ipHi] + Fe2LPump[ipHi][ipLo];
				hydro.dstfe2lya += (float)(
					Fe2LevN[ipHi][ipLo].PopLo*Fe2LPump[ipLo][ipHi] -
					Fe2LevN[ipHi][ipLo].PopHi*Fe2LPump[ipHi][ipLo]); 
			}
		}
		/* the destruction prob comes from
		 * dest rate = n(2p) * A(21) * PDest */
		pop = EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].PopHi*dense.xIonDense[ipHYDROGEN][1];
		if( pop > SMALLFLOAT )
		{
			hydro.dstfe2lya /= (float)(pop * EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Aul);
		}
		else
		{
			hydro.dstfe2lya = 0.;
		}
		/* NB - do not update hydro.dstfe2lya here if large feii not on since
		 * done in feii overlap */
	}

	{
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC)
		{
			/*lint -e644 EnergyWN not init */
			fprintf(ioQQQ," sum all %.1e dest rate%.1e escR= %.1e\n", 
				EnergyWN,hydro.dstfe2lya, 
				EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pesc);
			/*lint +e644 EnergyWN not init */
		}
	}

	/* next two blocks determine departure coefficients for the atom */
	
	/* first sum up partition function for the model atom  */
	Fe2DepCoef[0] = 1.0;
	sum = 1.0;
	for( i=1; i < FeII.nFeIILevel; i++ )
	{
		/* Boltzmann factor relative to ground times ratio of statistical weights */
		Fe2DepCoef[i] = Fe2DepCoef[0]*FeIIBoltzmann[i]*
			Fe2LevN[i][0].gHi/ Fe2LevN[1][0].gLo;
		/* this sum is the partition function for the atom */
		sum += Fe2DepCoef[i];
	}

	/* now renormalize departure coefficients with partition function */
	for( i=0; i < FeII.nFeIILevel; i++ )
	{
		/* divide by total partition function, Fe2DepCoef is now the fraction of the
		 * population that is in this level in TE */
		Fe2DepCoef[i] /= sum;

		/* convert to true departure coefficient */
		if( Fe2DepCoef[i]>SMALLFLOAT )
		{
			Fe2DepCoef[i] = yVector[i]/Fe2DepCoef[i];
		}
		else
		{
			Fe2DepCoef[i] = 0.;
		}
	}

	/*calculate cooling, heating, and cooling derivative */
	/* this is net cooling, cooling minus heating */
	FeII.Fe2_large_cool = 0.0f;
	/* this is be heating, not heating minus cooling */
	FeII.Fe2_large_heat = 0.f;
	/* derivative of cooling */
	FeII.ddT_Fe2_large_cool = 0.0f;

	/* compute heating and cooling due to model atom */
	for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel; ipHi++ )
		{
			double OneLine;

			/* net cooling due to single line */
			OneLine = (Fe2Coll[ipLo][ipHi]*Fe2LevelPop[ipLo] - Fe2Coll[ipHi][ipLo]*Fe2LevelPop[ipHi])*
				dense.eden*Fe2LevN[ipHi][ipLo].EnergyErg;

			/* net cooling due to this line */
			FeII.Fe2_large_cool += (float)MAX2(0., OneLine);

			/* net heating due to line */
			FeII.Fe2_large_heat += (float)MAX2(0., -OneLine);

			/* derivative of FeII cooling */
			if( OneLine >= 0. )
			{
				/* net coolant, exponential dominates deriv */
				FeII.ddT_Fe2_large_cool += (float)OneLine*
					(Fe2LevN[ipHi][0].EnergyK*thermal.tsq1 - thermal.halfte);
			}
			else
			{
				/* net heating, te^-1/2 dominates */
				FeII.ddT_Fe2_large_cool -= (float)OneLine*thermal.halfte;
			}
		}
	}

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

	return;
}

/*FeIICollRatesBoltzmann evaluate collision strenths, 
 * both interpolating on r-mat and creating g-bar 
 * find Boltzmann factors, evaluate collisional rate coefficients */
/* >>05 dec 03, verified that this rountine only changes temperature dependent
 * quantities - nothing with temperature dependence is done */
static void FeIICollRatesBoltzmann(void)
{
	/* this will be used to only reevaluate cs when temp changes */
	/* >>chng 00 jun 02, demoted next to float, PvH */
	static double OldTemp = -1.;
	long int i,
		ipLo ,
		ipHi;
	float ag, 
	  cg, 
	  dg, 
	  gb, 
	  y;
	float FracLowTe , FracHighTe;
	static float tt[8]={1e3f,3e3f,5e3f,7e3f,1e4f,12e3f,15e3f,2e4f};
	long int ipTemp,
		ipTempp1 , 
		ipLoFe2, 
		ipHiFe2;

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

	/* return if temperature has not changed */
	/*lint -e777 test floats for equality*/
	if( phycon.te == OldTemp )
	{
#		ifdef DEBUG_FUN
		fputs( " <->FeIICollRatesBoltzmann()\n", debug_fp );
#		endif

		return;
	}
	/*lint +e777 test floats for equality*/
	OldTemp = phycon.te;

	for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel; ipHi++ )
		{
			if( ncs1[ipHi][ipLo] == 3 ) 
			{
				/************************forbidden tr**************************/
				ag = 0.15f;
				cg = 0.f;
				dg = 0.f;
			}
			/************************allowed tr*******************************/
			else if( ncs1[ipHi][ipLo] == 1 )
			{
				ag = 0.6f;
				cg = 0.f;
				dg = 0.28f;
			}
			/************************semiforbidden****************************/
			else if( ncs1[ipHi][ipLo] == 2 )
			{
				ag = 0.0f;
				cg = 0.1f;
				dg = 0.0f;
			}
			else
			{
				/* this branch is impossible, since cs must be 1, 2, or 3 */
				ag = 0.0f;
				cg = 0.1f;
				dg = 0.0f;
				fprintf(ioQQQ,">>>PROBLEM impossible ncs1 in FeIILevelPops\n");
			}

			/*y = 1.438826f*Fe2LevN[ipHi][ipLo].EnergyWN/ phycon.te;*/
			y = Fe2LevN[ipHi][ipLo].EnergyWN/ (float)phycon.te_wn;

			gb = (float)(ag + (-cg*POW2(y) + dg)*(log(1.0+1.0/y) - 0.4/
			  POW2(y + 1.0)) + cg*y);

			Fe2LevN[ipHi][ipLo].cs = 23.861f*1e5f*gb*
			  Fe2LevN[ipHi][ipLo].Aul*Fe2LevN[ipHi][ipLo].gHi/
			  POW3(Fe2LevN[ipHi][ipLo].EnergyWN);

			/* confirm that collision strength is positive */
			ASSERT( Fe2LevN[ipHi][ipLo].cs > 0.);

			/* g-bar cs becomes unphysically small for forbidden transitions -
			 * this sets a lower limit to the final cs - CS2SMALL is defined above */
			Fe2LevN[ipHi][ipLo].cs = MAX2( CS2SMALL, Fe2LevN[ipHi][ipLo].cs);
			/* this was the logic used in the old fortran version,
			 * and reproduces the results in Baldwin et al '96
			 if( Fe2LevN[ipHi][ipLo].cs < 1e-10 )
			{
				Fe2LevN[ipHi][ipLo].cs = 0.01f;
			}
			*/
		}
	}
	
	/* we will interpolate on the set of listed collision strengths -
	 * where in this set are we? */
	if( phycon.te <= tt[0] )
	{
		/* temperature is lower than lowest tabulated, use the
		 * lowest tabulated point */
		/* ipTemp usually points to the cell cooler than te, ipTemp+1 to one higher,
		 * here both are lowest */
		ipTemp = 0;
		ipTempp1 = 0;
		FracHighTe = 0.;
	}
	else if( phycon.te > tt[7] )
	{
		/* temperature is higher than highest tabulated, use the
		 * highest tabulated point */
		/* ipTemp usually points to the cell cooler than te, ipTemp+1 to one higher,
		 * here both are highest */
		ipTemp = 7;
		ipTempp1 = 7;
		FracHighTe = 0.;
	}
	else
	{
		/* where in the array is the temperature we need? */
		ipTemp = -1;
		for( i=0; i < 8-1; i++ )
		{
			if( phycon.te <= tt[i+1] )
			{
				ipTemp = i;
				break;
			}

		}
		/* this cannot possibly happen */
		if( ipTemp < 0 )
		{
			fprintf( ioQQQ, " Insanity while looking for temperature in coll str array, te=%g.\n", 
			  phycon.te );
			puts( "[Stop in FeIILevelPops]" );
			cdEXIT(EXIT_FAILURE);
		}
		/* ipTemp points to the cell cooler than te, ipTemp+1 to one higher,
		 * do linear interpolation between */
		ipTemp = i;
		ipTempp1 = i+1;
		FracHighTe = ((float)phycon.te - tt[ipTemp])/(tt[ipTempp1] - tt[ipTemp]);
	}

	/* this is fraction of final linear interpolated collision strength that
	 * is weighted by the lower bound cs */
	FracLowTe = 1.f - FracHighTe;

	for( ipHi=1; ipHi < NPRADDAT; ipHi++ )
	{
		for( ipLo=0; ipLo < ipHi; ipLo++ )
		{
			/* ipHiFe2 should point to upper level of this pair, and
			 * ipLoFe2 should point to lower level */
			ipHiFe2 = MAX2( nnPradDat[ipHi] , nnPradDat[ipLo] );
			ipLoFe2 = MIN2( nnPradDat[ipHi] , nnPradDat[ipLo] );
			ASSERT( ipHiFe2-1 < NFE2LEVN );
			ASSERT( ipHiFe2-1 >= 0 );
			ASSERT( ipLoFe2-1 < NFE2LEVN );
			ASSERT( ipLoFe2-1 >= 0 );

			/* do linear interpolation for CS, these are dimensioned NPRADDAT = 159 */
			if( ipHiFe2-1 < FeII.nFeIILevel )
			{
				/* do linear interpolation */
				Fe2LevN[ipHiFe2-1][ipLoFe2-1].cs = 
					sPradDat[ipHi][ipLo][ipTemp] * FracLowTe + 
					sPradDat[ipHi][ipLo][ipTempp1] * FracHighTe ;

				/* confirm that this is positive */
				ASSERT( Fe2LevN[ipHiFe2-1][ipLoFe2-1].cs > 0. );
			}
		}
	}

	/* create Boltzmann factors for all levels */
	FeIIBoltzmann[0] = 1.0;
	for( ipHi=1; ipHi < FeII.nFeIILevel; ipHi++ )
	{
		/*FeIIBoltzmann[ipHi] = sexp( 1.438826f*Fe2LevN[ipHi][0].EnergyWN/phycon.te );*/
		/* >>chng 99 may 21, from above to following slight different number from phyconst.h 
		 * >>chng 05 dec 01, from sexp to dsexp to avoid all 0 Boltzmann factors in low temps
		 * now that FeII is on all the time */
		FeIIBoltzmann[ipHi] = dsexp( Fe2LevN[ipHi][0].EnergyWN/phycon.te_wn );
	}

	/* now possibly trim down atom if Boltzmann factors for upper levels are zero */
	ipHi = FeII.nFeIILevel - 1;
	while( FeIIBoltzmann[ipHi] == 0. && ipHi > 0 )
	{
		--ipHi;
	}
	/* ipHi now is the highest level with finite Boltzmann factor - 
	 * use this as the number of levels 
	 *>>chng 05 dec 01, from <=1 to <=0 - func_map does 10K, and large FeII had never
	 * been tested in that limit.  1 is ok. */
	if( ipHi <= 0 )
	{
		/* this is basically a sanity check and can't possibly happen */
		fprintf( ioQQQ, " FeIILevelPops trimmed Boltzmann factors too low, Te=%.3e\n",
			phycon.te );
		puts( "[Stop in FeIILevelPops]" );
		cdEXIT(EXIT_FAILURE);
	}
	else
	{
		/* is nearly all cases this does nothing since ipHi and FeII.nFeIILevel
		 * are equal . . . */
		FeII.nFeIILevel = ipHi+1;
		/*fprintf(ioQQQ," levels reset to %li\n",FeII.nFeIILevel);*/
	}

	/* debug code to print out the collision strengths for some levels */
	{
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC)
		{
			for( ipLo=0; ipLo < 52; ipLo++ )
			{
				fprintf(ioQQQ,"%e %e\n", 
					Fe2LevN[51][ipLo].cs,Fe2LevN[52][ipLo].cs);
			}
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* collisional excitation and deexcitation */
	for( ipLo=0; ipLo<FeII.nFeIILevel; ++ipLo)
	{
		for( ipHi=ipLo+1; ipHi<FeII.nFeIILevel; ++ipHi )
		{
			/* collisional deexcitation rate coefficient from ipHi to lower ipLo
			 * note that it needs eden to become rate
			 * units cm3 s-1 */
			Fe2Coll[ipHi][ipLo] = (float)(COLL_CONST/phycon.sqrte*Fe2LevN[ipHi][ipLo].cs/
			  Fe2LevN[ipHi][ipLo].gHi);

			/* collisional excitation rate coefficient from ipLo to upper ipHi 
			 * units cm3 s-1 */
			Fe2Coll[ipLo][ipHi] = (float)(Fe2Coll[ipHi][ipLo]*FeIIBoltzmann[ipHi]/FeIIBoltzmann[ipLo]*
				Fe2LevN[ipHi][ipLo].gHi/Fe2LevN[ipHi][ipLo].gLo);
		}
	}

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

	return;
}

/*
 *=====================================================================
 */
/*subroutine FeIIPrint PhotOccNum_at_nu raspechatki naselennostej v cloudy.out ili v
 * otdel'nom file unit=33
 *!!nado takzhe vklyuchit raspechatku iz perekrytiya linii */
/*FeIIPrint - print output from large feii atom, called by prtzone */
void FeIIPrint(void)
{

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

	return;
}

/*
 *=====================================================================
 */
/*FeIISumBand, sum up large FeII emission over certain bands, called in lineset4
 * units are erg s-1 cm-3, same units as xIntensity in line structure */
double FeIISumBand(
	/* lower and upper range to wavelength */
	float wl1, 
	float wl2)
{
	long int ipHi, 
	  ipLo;
	double SumBandFe2_v;
	/* >>chng 00 jun 02, demoted next two to float, PvH */
	float ahi, 
	  alo;

#	ifdef DEBUG_FUN
	fputs( "<+>FeIISumBand()\n", debug_fp );
#	endif
	/*sum up large FeII emission over certain bands */

	if( dense.xIonDense[ipIRON][1] == 0. )
	{
#		ifdef DEBUG_FUN
		fputs( " <->FeIISumBand()\n", debug_fp );
#		endif

		return( 0. );
	}
	else
	{

		/* convert to line energy in ergs */
		ahi = 1.99e-8f/wl1;
		alo = 1.99e-8f/wl2;

		SumBandFe2_v = 0.;
		for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ )
		{
			for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel; ipHi++ )
			{
				if( Fe2LevN[ipHi][ipLo].EnergyErg >= alo && 
				  Fe2LevN[ipHi][ipLo].EnergyErg <= ahi )
					SumBandFe2_v += Fe2LevN[ipHi][ipLo].xIntensity;
			}
		}

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

		return( SumBandFe2_v );
	}
}

/*
 *=====================================================================
 */
/*FeIITauInc called once per zone in RT_tau_inc to increment large FeII atom line optical depths */
void FeIITauInc(void)
{
	long int ipHi, 
	  ipLo;

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

	for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel; ipHi++ )
		{
			/* >>chng 04 aug 28, add test on bogus line */
			if( Fe2LevN[ipHi][ipLo].ipCont > 0 )
				RT_line_one_tauinc( &Fe2LevN[ipHi][ipLo] ,
				-8 , -8 , ipHi , ipLo );
		}
	}

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

	return;
}

/*
 *=====================================================================
 */
/*FeII_RT_tau_reset reset optical depths for large FeII atom, called by update after each iteration */
void FeII_RT_tau_reset(void)
{
	long int ipHi, 
	  ipLo;

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

	/*fprintf(ioQQQ,"DEBUG FeIITauAver1 %li %.2e %.2e nFeIILevel %li \n",
		iteration,
		Fe2LevN[9][0].TauIn ,
		Fe2LevN[9][0].TauTot, 
		FeII.nFeIILevel);*/

	/* called at end of iteration */
	/* >>chng 05 dec 07, had been FeII.nFeIILevel but this may have been trimmed down
	 * if previous iteration went very deep - not reset until FeIIReset called */
	for( ipHi=1; ipHi < FeII.nFeIILevelAlloc; ipHi++ )
	{
		for( ipLo=0; ipLo < ipHi; ipLo++ )
		{
			RT_line_one_tau_reset( &Fe2LevN[ipHi][ipLo],0.5 );
		}
	}
	/*fprintf(ioQQQ,"DEBUG FeIITauAver2 %li %.2e %.2e\n",
		iteration,
		Fe2LevN[9][0].TauIn ,
		Fe2LevN[9][0].TauTot);*/

	/* call the line overlap routine */
	FeIIOvrLap();
	/*fprintf(ioQQQ,"DEBUG FeIITauAver3 %li %.2e %.2e\n",
		iteration,
		Fe2LevN[9][0].TauIn ,
		Fe2LevN[9][0].TauTot);*/

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

	return;
}

/*
 *=====================================================================
 */
/*FeIIPoint called by ContCreatePointers to create pointers for lines in large FeII atom */
void FeIIPoint(void)
{
	long int ipHi, 
	  ip, 
	  ipLo;

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

	/* routine called when cloudy sets continuum array indices for Fe2 lines, 
	 * once per coreload */
	for( ipLo=0; ipLo < FeII.nFeIILevel-1; ipLo++ )
	{
		for( ipHi=ipLo+1; ipHi < FeII.nFeIILevel; ipHi++ )
		{

			/* >>chng 02 feb 11, set continuum index to negative value for fake transition */
			if( fabs(Fe2LevN[ipHi][ipLo].Aul- 1e-5 ) > 1e-8 ) 
			{
				ip = ipoint(Fe2LevN[ipHi][ipLo].EnergyWN * WAVNRYD);
				Fe2LevN[ipHi][ipLo].ipCont = ip;

				/* do not over write other pointers with feii since feii is everywhere */
				if( strcmp( rfield.chLineLabel[ip-1], "    " ) == 0 )
					strcpy( rfield.chLineLabel[ip-1], "FeII" );

				Fe2LevN[ipHi][ipLo].ipFine = ipFineCont( Fe2LevN[ipHi][ipLo].EnergyWN * WAVNRYD);
			}
			else
			{
				Fe2LevN[ipHi][ipLo].ipCont = -1;
				Fe2LevN[ipHi][ipLo].ipFine = -1;
			}

			Fe2LevN[ipHi][ipLo].dampXvel = 
				(float)(Fe2LevN[ipHi][ipLo].Aul/
			  Fe2LevN[ipHi][ipLo].EnergyWN/PI4);

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

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

			/* energy of photon in ergs */
			Fe2LevN[ipHi][ipLo].EnergyErg = 
				(float)(ERG1CM*Fe2LevN[ipHi][ipLo].EnergyWN);
		}
	}

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

	return;
}

/*
 *=====================================================================
 */
/*FeIIAccel called by rt_line_driving to compute radiative acceleration due to FeII lines */
void FeIIAccel(double *fe2drive)
{
	long int ipHi, 
	  ipLo;

#	ifdef DEBUG_FUN
	fputs( "<+>FeIIAccel()\n", debug_fp );
#	endif
	/*compute acceleration due to large FeII atom */

	/* this routine computes the line driven radiative acceleration
	 * due to large FeII atom*/

	*fe2drive = 0.;
	for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ )
	{
		for( ipHi=ipLo+1; ipHi < FeII.nFeIILevel; ipHi++ )
		{
			*fe2drive += Fe2LevN[ipHi][ipLo].pump*
			  Fe2LevN[ipHi][ipLo].EnergyErg*Fe2LevN[ipHi][ipLo].PopOpc;
		}
	}

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

	return;
}

/*
 *=====================================================================
 */
/*FeIIRTMake called by RT_line_all, does large FeII atom radiative transfer */
void FeIIRTMake( int lgDoEsc , int lgUpdateFineOpac )
{
	long int ipHi, 
	  ipLo;

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

	if( trace.lgTrace )
	{
		fprintf( ioQQQ,"   FeIIRTMake called\n");
	}

	/* this routine drives calls to make RT relations with large FeII atom */
	for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel; ipHi++ )
		{
			/* only evaluate real transitions */
			if( Fe2LevN[ipHi][ipLo].ipCont > 0 ) 
			{
				RT_line_one( &Fe2LevN[ipHi][ipLo] , lgDoEsc , lgUpdateFineOpac,TRUE);
			}
		}
	}

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

	return;
}

/*
 *=====================================================================
 */
/* called in LineSet4 to add FeII lines to save array */
void FeIIAddLines( void )
{
	long int ipHi, 
	  ipLo;

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

	/* this routine puts the emission from the large FeII atom
	 * into an array to save and integrate them*/

	/* add lines option called from lines, add intensities into storage array */

	/* routine is called three different ways, ipass < 0 before space allocated,
	 * =0 when time to generate labels (and we zero out array here), and ipass>0
	 * when time to save intensities */
	if( LineSave.ipass == 0 )
	{
		/* we were called by lines, and we want to zero out Fe2SavN */
		for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ )
		{
			for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel; ipHi++ )
			{
				Fe2SavN[ipHi][ipLo] = 0.;
			}
		}
	}

	/* this call during calculations, save intensities */
	else if( LineSave.ipass == 1 )
	{
		/* total emission from vol */
		for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ )
		{
			for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel; ipHi++ )
			{
				Fe2SavN[ipHi][ipLo] += 
					radius.dVeff*Fe2LevN[ipHi][ipLo].xIntensity;
			}
		}
	}

	{
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC /*&& (iteration==2)*/ )
		{
			fprintf(ioQQQ," 69-1\t%li\t%e\n", nzone , Fe2SavN[68][0] );
		}
	}

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

	return;
}

/*FeIIPunchLevels punch level energies and stat weights */
void FeIIPunchLevels(
  /* file we will punch to */
  FILE *ioPUN )
{

	long int ipHi;

#	ifdef DEBUG_FUN
	fputs( "<+>FeIIPunchLevels()\n", debug_fp );
#	endif
	
	fprintf(ioPUN , "%.2f\t%li\n", 0., (long)Fe2LevN[1][0].gLo );
	for( ipHi=1; ipHi < FeII.nFeIILevel; ++ipHi )
	{
		fprintf(ioPUN , "%.2f\t%li\n", 
			Fe2LevN[ipHi][0].EnergyWN ,
			(long)Fe2LevN[ipHi][0].gHi);
	}

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

	return;
}

/*FeIIPunchColden punch level energies, stat weights, column density */
void FeIIPunchColden(
  /* file we will punch to */
  FILE *ioPUN )
{

	long int n;

#	ifdef DEBUG_FUN
	fputs( "<+>FeIIPunchColden()\n", debug_fp );
#	endif
	
	fprintf(ioPUN , "%.2f\t%li\t%.3e\n", 0., (long)Fe2LevN[1][0].gLo , Fe2ColDen[0]);
	for( n=1; n < FeII.nFeIILevel; ++n )
	{
		fprintf(ioPUN , "%.2f\t%li\t%.3e\n", 
			Fe2LevN[n][0].EnergyWN ,
			(long)Fe2LevN[n][0].gHi,
			Fe2ColDen[n]);
	}

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

	return;
}

/*
 *=====================================================================
 */
/*FeIIPunchLines punch accumulated FeII intensities, at end of calculation,
	called by dopunch to punch them out,
	punch turned on with punch lines command */
void FeIIPunchLines(
  /* file we will punch to */
  FILE *ioPUN )
{
	long int MaseHi, 
	  MaseLow, 
	  ipHi, 
	  ipLo;
	double hbeta, absint , renorm;
	/* >>chng 00 jun 02, demoted next two to float, PvH */
	float TauMase, 
	  thresh;

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

	/* this routine puts the emission from the large FeII atom
	 * into a line array, and eventually will punch it out */
	 
	/* get the normalization line */
	if( LineSv[LineSave.ipNormWavL].sumlin > 0. )
	{
		renorm = LineSave.ScaleNormLine/LineSv[LineSave.ipNormWavL].sumlin;
	}
	else
	{
		renorm = 1.;
	}

	fprintf( ioPUN, " up low log I, I/I(LineSave), Tau\n" );

	/* first look for any masing lines */
	MaseLow = -1;
	MaseHi = -1;
	TauMase = 0.f;
	for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel; ipHi++ )
		{
			if( Fe2LevN[ipHi][ipLo].TauIn < TauMase )
			{
				TauMase = Fe2LevN[ipHi][ipLo].TauIn;
				MaseLow = ipLo;
				MaseHi = ipHi;
			}
		}
	}

	if( TauMase < 0.f )
	{
		fprintf( ioPUN, " Most negative optical depth was %4ld%4ld%10.2e\n", 
		  MaseHi, MaseLow, TauMase );
	}

	/* now print actual line intensities, Hbeta first */
	if( cdLine("TOTL", 4861 , &hbeta , &absint)<=0 )
	{
		fprintf( ioQQQ, " FeIILevelPops could not find Hbeta with cdLine.\n" );
		puts( "[Stop in FeIILevelPops]" );
		cdEXIT(EXIT_FAILURE);
	}

	fprintf( ioPUN, "Hbeta=%7.3f %10.2e\n", 
	  absint , 
	  hbeta );

	if( renorm > SMALLFLOAT )
	{
		/* this is threshold for faintest line, normally 0, set with 
		 * number on punch verner command */
		thresh = FeII.fe2thresh/(float)renorm;
	}
	else
	{
		thresh = 0.f;
	}

	for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel; ipHi++ )
		{
			/* fe2ener(1) and (2) are lower and upper limits to range of
			 * wavelengths of interest.  entered in ryd with
			 * punch feii verner command, where they are converted
			 * to wavenumbers, */
			if( (Fe2SavN[ipHi][ipLo] > thresh && 
				Fe2LevN[ipHi][ipLo].EnergyWN > FeII.fe2ener[0]) &&
				Fe2LevN[ipHi][ipLo].EnergyWN < FeII.fe2ener[1] )
			{
				if( FeII.lgShortFe2 )
				{
					/* short printout does not include rel intensity or optical dep */
					fprintf( ioPUN, "%3ld%3ld%7.3f\n", 
					  ipHi+1, ipLo+1, 
					  log10(MAX2(1e-37,Fe2SavN[ipHi][ipLo])) + radius.Conv2PrtInten );
				}
				else
				{
					/* long printout does */
					fprintf( ioPUN, "%3ld%3ld %8.3f %10.2e %10.2e\n", 
					  ipHi+1, ipLo+1, 
					  log10(MAX2(1e-37,Fe2SavN[ipHi][ipLo])) + radius.Conv2PrtInten, 
					  Fe2SavN[ipHi][ipLo]*renorm ,
					  Fe2LevN[ipHi][ipLo].TauIn );
				}
			}
		}
	}

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

	return;
}


/*
 *=====================================================================
 */
/*FeII_LineZero zero out storage for large FeII atom, called in tauout */
void FeII_LineZero(void)
{
	long int ipHi, 
	  ipLo;

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

	/* this routine is called in routine zero and it
	* zero's out various elements of the FeII arrays
	* it is called on every iteration
	* */
	for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel; ipHi++ )
		{
			/* >>chng 03 feb 14, use EmLineZero rather than explicit sets */
			EmLineZero(&Fe2LevN[ipHi][ipLo] );
		}
	}

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

	return;
}
/*
 *=====================================================================
 */
/*FeIIIntenZero zero out storage for large FeII atom, called in tauout */
void FeIIIntenZero(void)
{
	long int ipHi, 
	  ipLo;

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

	/* this routine is called by routine cool_iron and it
	 * zero's out various elements of the FeII arrays
	 * */
	Fe2LevelPop[0] = 0.;
	for( ipHi=1; ipHi < FeII.nFeIILevel; ipHi++ )
	{
		/* >>chng 05 dec 14, zero Fe2LevelPop */
		Fe2LevelPop[ipHi] = 0.;
		for( ipLo=0; ipLo < ipHi; ipLo++ )
		{

			/* population of lower level with correction for stim emission */
			Fe2LevN[ipHi][ipLo].PopOpc = 0.;

			/* population of lower level */
			Fe2LevN[ipHi][ipLo].PopLo = 0.;

			/* population of upper level */
			Fe2LevN[ipHi][ipLo].PopHi = 0.;

			/* following two heat exchange excitation, deexcitation */
			Fe2LevN[ipHi][ipLo].cool = 0.;
			Fe2LevN[ipHi][ipLo].heat = 0.;

			/* intensity of line */
			Fe2LevN[ipHi][ipLo].xIntensity = 0.;

			Fe2LevN[ipHi][ipLo].phots = 0.;
			Fe2LevN[ipHi][ipLo].ots = 0.;
			Fe2LevN[ipHi][ipLo].ColOvTot = 0.;
		}
	}

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

	return;
}


/*
 *=====================================================================
 */
/* fill in IR lines from lowest 16 levels - these are predicted in Fe2Lev16 */
void FeIIFillLow16(void)
{

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

	/* this is total cooling due to 16 level atom that would have been computed there */
	/*FeII.Fe2_16levl_cool = 0.;*/

	/* all transitions within 16 levels of ground term */
	FeII.fe21308 = Fe2LevN[12][7].xIntensity;
	FeII.fe21207 = Fe2LevN[11][6].xIntensity;
	FeII.fe21106 = Fe2LevN[10][5].xIntensity;
	FeII.fe21006 = Fe2LevN[9][5].xIntensity;
	FeII.fe21204 = Fe2LevN[11][3].xIntensity;
	FeII.fe21103 = Fe2LevN[10][2].xIntensity;
	FeII.fe21104 = Fe2LevN[10][3].xIntensity;
	FeII.fe21001 = Fe2LevN[9][0].xIntensity;
	FeII.fe21002 = Fe2LevN[9][1].xIntensity;
	FeII.fe20201 = Fe2LevN[1][0].xIntensity;
	FeII.fe20302 = Fe2LevN[2][1].xIntensity;
	FeII.fe20706 = Fe2LevN[6][5].xIntensity;
	FeII.fe20807 = Fe2LevN[7][6].xIntensity;
	FeII.fe20908 = Fe2LevN[8][7].xIntensity;
	FeII.fe21007 = Fe2LevN[9][6].xIntensity;
	FeII.fe21107 = Fe2LevN[10][6].xIntensity;
	FeII.fe21108 = Fe2LevN[10][7].xIntensity;
	FeII.fe21110 = Fe2LevN[10][9].xIntensity;
	FeII.fe21208 = Fe2LevN[11][7].xIntensity;
	FeII.fe21209 = Fe2LevN[11][8].xIntensity;
	FeII.fe21211 = Fe2LevN[11][10].xIntensity;
	FeII.fe21406 = Fe2LevN[13][5].xIntensity;
	FeII.fe21507 = Fe2LevN[14][6].xIntensity;
	FeII.fe21508 = Fe2LevN[14][7].xIntensity;
	FeII.fe21609 = Fe2LevN[15][8].xIntensity;

	/* these are unique to this large model atom */
	if( FeII.nFeIILevel > 43 )
	{
		/* NB do not forget to decrement by one when going from physical (in left variables)
		 * to c scale (in array) */
		FeII.fe25to6 = Fe2LevN[25-1][5].xIntensity;
		FeII.fe27to7 = Fe2LevN[27-1][6].xIntensity;
		FeII.fe28to8 = Fe2LevN[28-1][7].xIntensity;
		FeII.fe29to9 = Fe2LevN[29-1][8].xIntensity;
 		FeII.fe32to6 = Fe2LevN[32-1][5].xIntensity;
 		FeII.fe33to7 = Fe2LevN[33-1][6].xIntensity;
 		FeII.fe37to7 = Fe2LevN[37-1][6].xIntensity;
 		FeII.fe39to8 = Fe2LevN[39-1][7].xIntensity;
 		FeII.fe40to9 = Fe2LevN[40-1][8].xIntensity;
 		FeII.fe37to6 = Fe2LevN[37-1][5].xIntensity;
 		FeII.fe39to7 = Fe2LevN[39-1][6].xIntensity;
 		FeII.fe40to8 = Fe2LevN[40-1][7].xIntensity;
 		FeII.fe41to9 = Fe2LevN[41-1][8].xIntensity;
 		FeII.fe39to6 = Fe2LevN[39-1][5].xIntensity;
 		FeII.fe40to7 = Fe2LevN[40-1][6].xIntensity;
 		FeII.fe41to8 = Fe2LevN[41-1][7].xIntensity;

		FeII.fe42to6 = Fe2LevN[42-1][5].xIntensity;
		FeII.fe43to7 = Fe2LevN[43-1][6].xIntensity;
		FeII.fe42to7 = Fe2LevN[42-1][6].xIntensity;
		FeII.fe36to2 = Fe2LevN[36-1][1].xIntensity;
		FeII.fe36to3 = Fe2LevN[36-1][2].xIntensity;
		/*lint -e778 const expression eval to 0 */
		FeII.fe32to1 = Fe2LevN[32-1][0].xIntensity;
		/*lint +e778 const expression eval to 0 */
		FeII.fe33to2 = Fe2LevN[33-1][1].xIntensity;
		FeII.fe36to5 = Fe2LevN[36-1][4].xIntensity;
		FeII.fe32to2 = Fe2LevN[32-1][1].xIntensity;
		FeII.fe33to3 = Fe2LevN[33-1][2].xIntensity;
		FeII.fe30to3 = Fe2LevN[30-1][2].xIntensity;
		FeII.fe33to6 = Fe2LevN[33-1][5].xIntensity;
		FeII.fe24to2 = Fe2LevN[24-1][1].xIntensity;
		FeII.fe32to7 = Fe2LevN[32-1][6].xIntensity;
		FeII.fe35to8 = Fe2LevN[35-1][7].xIntensity;
		FeII.fe34to8 = Fe2LevN[34-1][7].xIntensity;
		FeII.fe27to6 = Fe2LevN[27-1][5].xIntensity;
		FeII.fe28to7 = Fe2LevN[28-1][6].xIntensity;
		FeII.fe30to8 = Fe2LevN[30-1][7].xIntensity;
		FeII.fe24to6 = Fe2LevN[24-1][5].xIntensity;
		FeII.fe29to8 = Fe2LevN[29-1][7].xIntensity;
		FeII.fe24to7 = Fe2LevN[24-1][6].xIntensity;
		FeII.fe22to7 = Fe2LevN[22-1][6].xIntensity;
		FeII.fe38to11 =Fe2LevN[38-1][11-1].xIntensity;
		FeII.fe19to8 = Fe2LevN[19-1][7].xIntensity;
		FeII.fe17to6 = Fe2LevN[17-1][5].xIntensity;
		FeII.fe18to7 = Fe2LevN[18-1][6].xIntensity;
		FeII.fe18to8 = Fe2LevN[18-1][7].xIntensity;
		FeII.fe17to7 = Fe2LevN[17-1][6].xIntensity;
	}
	else
	{
		FeII.fe25to6 = 0.;
		FeII.fe27to7 = 0.;
		FeII.fe28to8 = 0.;
		FeII.fe29to9 = 0.;
 		FeII.fe32to6 = 0.;
 		FeII.fe33to7 = 0.;
 		FeII.fe37to7 = 0.;
 		FeII.fe39to8 = 0.;
 		FeII.fe40to9 = 0.;
 		FeII.fe37to6 = 0.;
 		FeII.fe39to7 = 0.;
 		FeII.fe40to8 = 0.;
 		FeII.fe41to9 = 0.;
 		FeII.fe39to6 = 0.;
 		FeII.fe40to7 = 0.;
 		FeII.fe41to8 = 0.;

		FeII.fe42to6 = 0.;
		FeII.fe43to7 = 0.;
		FeII.fe42to7 = 0.;
		FeII.fe36to2 = 0.;
		FeII.fe36to3 = 0.;
		FeII.fe32to1 = 0.;
		FeII.fe33to2 = 0.;
		FeII.fe36to5 = 0.;
		FeII.fe32to2 = 0.;
		FeII.fe33to3 = 0.;
		FeII.fe30to3 = 0.;
		FeII.fe33to6 = 0.;
		FeII.fe24to2 = 0.;
		FeII.fe32to7 = 0.;
		FeII.fe35to8 = 0.;
		FeII.fe34to8 = 0.;
		FeII.fe27to6 = 0.;
		FeII.fe28to7 = 0.;
		FeII.fe30to8 = 0.;
		FeII.fe24to6 = 0.;
		FeII.fe29to8 = 0.;
		FeII.fe24to7 = 0.;
		FeII.fe22to7 = 0.;
		FeII.fe38to11 =0.;
		FeII.fe19to8 = 0.;
		FeII.fe17to6 = 0.;
		FeII.fe18to7 = 0.;
		FeII.fe18to8 = 0.;
		FeII.fe17to7 = 0.;
	}

	if( FeII.nFeIILevel > 80 )
	{
		FeII.fe80to28 = Fe2LevN[80-1][28-1].xIntensity;
	}
	else
	{
		FeII.fe80to28 =0.;
	}

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

	return;
}
/*
 *=====================================================================
 * punch line data for FeII atom */
void FeIIPunData(
	/* io unit for punch */
	FILE* ioPUN ,
	/* punch all levels if true, only subset if false */
	int lgDoAll )
{
	long int ipLo , ipHi;

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

	if( lgDoAll )
	{
		fprintf( ioQQQ, 
			" FeIIPunData ALL option not implemented yet 1\n" );
		puts( "[Stop in FeIIPunData]" );
		cdEXIT(EXIT_FAILURE);
	}
	else
	{
		long int nSkip=0, limit=MIN2(64, FeII.nFeIILevel);
		/* false, only punch subset in above init */
		/* first 64 do all lines */
		for( ipHi=1; ipHi<limit; ++ipHi )
		{
			for( ipLo=0; ipLo<ipHi; ++ipLo )
			{
				fprintf(ioPUN,"%4li%4li ",ipLo,ipHi );
				Punch1LineData( &Fe2LevN[ipHi][ipLo] , ioPUN);
			}
		}
		fprintf( ioPUN , "\n");

		if( limit==64 )
		{
			/* higher than 64 only do real transitions */
			for( ipHi=limit; ipHi<FeII.nFeIILevel; ++ipHi )
			{
				for( ipLo=limit; ipLo<ipHi; ++ipLo )
				{
					if( ncs1[ipHi][ipLo] == 3 && 
						(fabs(Fe2LevN[ipHi][ipLo].Aul-1e-5) < 1e-8 ) )
					{
						++nSkip;
					}
					else
					{
						/* add one so that one atomic, not C, scale */
						fprintf(ioPUN,"%4li%4li ",ipLo+1,ipHi+1 );
						Punch1LineData( &Fe2LevN[ipHi][ipLo] , ioPUN);
					}
				}
			}
			fprintf( ioPUN , " %li lines skiped\n",nSkip);
		}
	}

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

	return;
}

/*
 *=====================================================================
 */
void FeIIPunDepart(
	/* io unit for punch */
	FILE* ioPUN ,
	/* punch all levels if true, only subset if false */
	int lgDoAll )
{
	/* numer of levels with dep coef that we will punch out */
#	define NLEVDEP 11
	/* these are the levels on the physical, not c, scale (count from 1) */
	const int LevDep[NLEVDEP]={1,10,25,45,64,124,206,249,295,347,371};
	long int i;
	static int lgFIRST=TRUE;

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

	/* on first call only, print levels that we will punch later */
	if( lgFIRST && !lgDoAll )
	{
		/* but all the rest do */
		for( i=0; i<NLEVDEP; ++i )
		{
			fprintf( ioPUN , "%i\t", LevDep[i] );
		}
		fprintf( ioPUN , "\n");
		lgFIRST = FALSE;
	}

	if( lgDoAll )
	{
		/* true, punch all levels, one per line */
		for( i=1; i<=FeII.nFeIILevel; ++i )
		{
			FeIIPun1Depart( ioPUN , i );
			fprintf( ioPUN , "\n");
		}
	}
	else
	{
		/* false, only punch subset in above init */
		for( i=0; i<NLEVDEP; ++i )
		{
			FeIIPun1Depart( ioPUN , LevDep[i] );
			fprintf( ioPUN , "\t");
		}
		fprintf( ioPUN , "\n");
	}

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

	return;
}

/*
 *=====================================================================
 */
void FeIIPun1Depart( 
	/* the io unit where the print should be directed */
	FILE * ioPUN , 
	/* the physical (not c) number of the level */
	long int nPUN )
{

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

	ASSERT( nPUN > 0 );
	/* need real assert to keep lint happy */
	assert( ioPUN != NULL );

	/* print the level departure coef on ioPUN if the level was computed,
	 * print a zero if it was not */
	if( nPUN <= FeII.nFeIILevel )
	{
		/*lint -e668 ioPUN possibly null */
		fprintf( ioPUN , "%e ",Fe2DepCoef[nPUN-1] );
	}
	else
	{
		fprintf( ioPUN , "%e ",0. );
		/*lint +e668 ioPUN possibly null */
	}

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

	return;
}

/*
 *=====================================================================
 */
void FeIIReset(void)
{

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

	/* space has been allocated, so reset number of levels to whatever it was */
	FeII.nFeIILevel = FeII.nFeIILevelAlloc;

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

/*
 *=====================================================================
 */
/*FeIIZero initialize some variables, called by zero one time before commands parsed*/
void FeIIZero(void)
{

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

	/* heating, cooling, and deriv wrt temperature */
	FeII.Fe2_large_cool = 0.;
	FeII.ddT_Fe2_large_cool = 0.;
	FeII.Fe2_large_heat = 0.;

	/* flag saying that lya pumping of feii in large atom is turned on */
	FeII.lgLyaPumpOn = TRUE;

	/*FeII. lgFeIION = FALSE;*/
	/* >>chng 05 nov 29, test effects of always including FeII ground term with large atom
	 * but if ground term only is done, still also do simple UV approximation */
	/*FeII. lgFeIION = TRUE;*/
	/* will not compute large feii atom */
	FeII.lgFeIILargeOn = FALSE;

	/* energy range of large FeII atom is zero to infinity */
	FeII.fe2ener[0] = 0.;
	FeII.fe2ener[1] = 1e8;

	/* pre mar 01, these had both been 1, ipPRD */
	/* resonance lines, ipCRD is -1 */
	FeII.ipRedisFcnResonance = ipCRD;
	/* subordinate lines, ipCRDW is 2 */
	FeII.ipRedisFcnSubordinate = ipCRDW;

	/* set zero for the threshold of weakest large FeII atom line to punch */
	FeII.fe2thresh = 0.;

	/* normally do not constantly reevaluate the atom, set TRUE with
	 * SLOW key on atom feii command */
	FeII.lgSlow = FALSE;

	/* option to print each call to FeIILevelPops, set with print option on atom feii */
	FeII.lgPrint = FALSE;

	/* option to only simulate calls to FeIILevelPops */
	FeII.lgSimulate = FALSE;

	/* set number of levels for the atom */
		/* number of levels for the large FeII atom, changed with the atom feii levels command */
	if( lgFeIIMalloc )
	{
		/* space has been allocated, so reset number of levels to whatever it was */
		FeII.nFeIILevel = FeII.nFeIILevelAlloc;
	}
	else
	{
		/* space not allocated yet, set to largest possible number of levels */
		FeII.nFeIILevel = NFE2LEVN;
		/* >>chng 05 nov 29, test effects of always including FeII ground term with large atom
		 * but if ground term only is done, still also do simple UV approximation 
		 * set this to only ground term, - will reset to NFE2LEVN when atom feii parsed if levels not set */
		FeII.nFeIILevel = 16;
	}

	/* lower and upper wavelength bounds, Angstroms, for the FeII continuum */
	FeII.fe2con_wl1 = 1000.;
	FeII.fe2con_wl2 = 7000.;
	FeII.nfe2con = 1000;

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

/*FeIIPunPop - punch level populations */
void FeIIPunPop(
	/* io unit for punch */
	FILE* ioPUN ,
	/* punch all levels if true, only subset if false */
	int lgDoAll )
{
	/* numer of levels with dep coef that we will punch out */
#	define NLEVPOP 11
	/* these are the levels on the physical, not c, scale (count from 1) */
	const int LevPop[NLEVPOP]={1,10,25,45,64,124,206,249,295,347,371};
	long int i;
	static int lgFIRST=TRUE;

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

	/* on first call only, print levels that we will punch later,
	 * but only if we will only punch selected levels*/
	if( lgFIRST && !lgDoAll )
	{
		/* but all the rest do */
		for( i=0; i<NLEVPOP; ++i )
		{
			fprintf( ioPUN , "%i\t", LevPop[i] );
		}
		fprintf( ioPUN , "\n");
		lgFIRST = FALSE;
	}

	if( lgDoAll )
	{
		/* true, punch all levels, one per line,
		 * both call with physical level so that list is physical */
		for( i=1; i<=FeII.nFeIILevel; ++i )
		{
			FeIIPun1Pop( ioPUN , i );
			fprintf( ioPUN , "\n");
		}
	}
	else
	{
		/* false, only punch subset in above init,
		 * both call with physical level so that list is physical  */
		for( i=0; i<NLEVPOP; ++i )
		{
			FeIIPun1Pop( ioPUN , LevPop[i] );
			fprintf( ioPUN , "\t");
		}
		fprintf( ioPUN , "\n");
	}

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

	return;
}

/*
 *=====================================================================
 */
void FeIIPun1Pop( 
	/* the io unit where the print should be directed */
	FILE * ioPUN , 
	/* the physical (not c) number of the level */
	long int nPUN )
{
#	ifdef DEBUG_FUN
	fputs( "<+>FeIIPun1Pop()\n", debug_fp );
#	endif

	ASSERT( nPUN > 0 );
	/* need real assert to keep lint happy */
	assert( ioPUN != NULL );

	/* print the level population on ioPUN if the level was computed,
	 * print a zero if it was not, 
	 * note that nPUN is on physical scale, so test is <=*/
	if( nPUN <= FeII.nFeIILevel )
	{
		fprintf( ioPUN , "%e ",Fe2LevelPop[nPUN-1] );
	}
	else
	{
		fprintf( ioPUN , "%e ",0. );
	}

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

	return;
}

/*
 *=====================================================================
 */
static int FeIIBandsCreate(
	/* chFile is optional filename, if void then use default bands,
	 * if not void then use file specified,
	 * return value is 0 for success, 1 for failure */
	 const char chFile[] )
{

	char chLine[FILENAME_PATH_LENGTH_2] , 
		chFilename[FILENAME_PATH_LENGTH_2] ,
		chFile1[FILENAME_PATH_LENGTH_2];
	FILE *ioDATA;

	int lgEOL;
	long int i,k;

	/* keep track of whether we have been called - want to be
	 * called a total of one time */
	static int lgCalled=FALSE;

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

	/* return previous number of bands if this is second or later call*/
	if( lgCalled )
	{
#		ifdef DEBUG_FUN
		fputs( " <->FeIIBandsCreate()\n", debug_fp );
#		endif
		/* success */
		return 0;
	}
	lgCalled = TRUE;

	/* use default filename if void string, else use file specified */
	if( strlen( chFile )==0 )
	{
		/* void string, use default name */
		strcpy( chFile1 , "bands_Fe2.dat" );
	}
	else
	{
		/* not void, use filename given us */
		strcpy( chFile1 , chFile );
	}

	/* get FeII band data */
	/* check on path if path set */
	/* path was parsed in getset */
	if( lgDataPathSet == TRUE )
	{
		/*path set, so look only there */
		strcpy( chFilename , chDataPath );
		strcat( chFilename , chFile1 );
	}
	else
	{
		/* path not set, check local space only */
		strcpy( chFilename , chFile1 );
	}

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, " FeIICreate opening %s:", chFile1 );
	}

	if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
	{
		if( lgDataPathSet == TRUE )
		{
			fprintf( ioQQQ, " FeIICreate could not open %s, even tried path.\n" , chFile1 );
			fprintf( ioQQQ, " path is *%s*\n",chDataPath );
			fprintf( ioQQQ, " final path is *%s*\n",chFilename );
		}
		else
		{
			fprintf( ioQQQ, " FeIICreate could not open %s\n" , chFile1 );
		}
		return 1;
	}

	ASSERT( ioDATA !=NULL);

	/* now count how many bands are in the file */
	nFeIIBands = 0;

	/* first line is a version number and does not count */
	if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
	{
		fprintf( ioQQQ, " FeIICreate could not read first line of %s.\n", chFile1 );
		return 1;
	}
	while( fgets( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
	{
		/* we want to count the lines that do not start with #
		 * since these contain data */
		if( chLine[0] != '#')
			++nFeIIBands;
	}

	/* now rewind the file so we can read it a second time*/
	if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
	{
		fprintf( ioQQQ, " FeIICreate could not rewind %s.\n", chFile1 );
		return 1;
	}

	if( (FeII_Bands = (float **)MALLOC(sizeof(float *)*(unsigned)(nFeIIBands) ) )==NULL )
	{
		fprintf( ioQQQ, 
			" FeIICreate could not MALLOC FeII_Bands 1\n" );
		return 1;
	}
	/* now make second dim, id wavelength, and lower and upper bounds */
	for( i=0; i<nFeIIBands; ++i )
	{
		if( (FeII_Bands[i] = (float *)MALLOC(sizeof(float)*(unsigned)(3) ) )==NULL )
		{
			fprintf( ioQQQ, 
				" FeIICreate could not MALLOC FeII_Bands 2\n" );
			return 1;
		}
	}

	/* first line is a version number - now confirm that it is valid */
	if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
	{
		fprintf( ioQQQ, " FeIICreate could not read first line of %s.\n", chFile1 );
		return 1;
	}
	i = 1;
	if( ( (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL) != 99 ) ||
	  ( (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL) != 12 ) ||
	  ( (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL) != 1 ) )
	{
		fprintf( ioQQQ, 
			" FeIICreate: the version of %s is not the current version.\n", chFile1 );
		return 1;
	}

	/* now read in data again, but save it this time */
	k = 0;
	while( fgets( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
	{
		/* we want to count the lines that do not start with #
		 * since these contain data */
		if( chLine[0] != '#')
		{
			i = 1;
			FeII_Bands[k][0] = (float)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				fprintf( ioQQQ, " There should have been a number on this band line 1.   Sorry.\n" );
				fprintf( ioQQQ, "string==%s==\n" ,chLine );
				return 1;
			}
			FeII_Bands[k][1] = (float)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				fprintf( ioQQQ, " There should have been a number on this band line 2.   Sorry.\n" );
				fprintf( ioQQQ, "string==%s==\n" ,chLine );
				return 1;
			}
			FeII_Bands[k][2] = (float)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				fprintf( ioQQQ, " There should have been a number on this band line 3.   Sorry.\n" );
				fprintf( ioQQQ, "string==%s==\n" ,chLine );
				return 1;
			}
			/*fprintf(ioQQQ,
			" band data %f %f %f \n", FeII_Bands[k][0],FeII_Bands[k][1],FeII_Bands[k][2]);*/
			++k;
		}
	}
	/* now validate this incoming data */
	for( i=0; i<nFeIIBands; ++i )
	{
		/* make sure all are positive */
		if( FeII_Bands[i][0] <=0. || FeII_Bands[i][1] <=0. || FeII_Bands[i][2] <=0. )
		{
			fprintf( ioQQQ, " FeII band %li has none positive entry.\n",i );
			return 1;
		}
		/* make sure bands bounds are in correct order, shorter - longer wavelength*/
		if( FeII_Bands[i][1] >= FeII_Bands[i][2] )
		{
			fprintf( ioQQQ, " FeII band %li has improper bounds.\n" ,i);
			return 1;
		}

	}

	fclose(ioDATA);
	
#	ifdef DEBUG_FUN
	fputs( " <->FeIIBandsCreate()\n", debug_fp );
#	endif

	/* return success */
	return 0;
}
/*
 *=====================================================================
 */
void AssertFeIIDep( double *pred , double *BigError , double *StdDev )
{
	long int n;
	double arg ,
		error ,
		sum2 ;

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

	if( FeII.lgSimulate )
	{

		*pred = 0.;
		*BigError = 0.;
		*StdDev = 0.;

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

		return;
	}

	/* sanity check */
	ASSERT( FeII.nFeIILevel > 0 );

	/* find sum of deviation of departure coef from unity */
	*BigError = 0.;
	*pred = 0.;
	sum2 = 0;
	for( n=0; n<FeII.nFeIILevel; ++n )
	{
		/* get mean */
		*pred += Fe2DepCoef[n];

		/* error is the largest deviation from unity for any single level*/
		error = fabs( Fe2DepCoef[n] -1. );
		/* remember biggest deviation */
		*BigError = MAX2( *BigError , error );

		/* get standard deviation */
		sum2 += POW2( Fe2DepCoef[n] );
	}

	/* now get standard deviation */
	arg = sum2 - POW2( *pred ) / (double)FeII.nFeIILevel ;
	ASSERT( (arg >= 0.) );
	*StdDev = sqrt( arg / (double)(FeII.nFeIILevel - 1.) );

	/* this is average value, should be unity */
	*pred /= (double)(FeII.nFeIILevel);

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

/*
 *=====================================================================
 */
/* do ots rates for FeII, called by RT_OTS */
void FeII_OTS( void )
{
	long int ipLo , 
		ipHi;

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

	for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel; ipHi++ )
		{
			/* >>chng 02 feb 11, skip bogus transitions */
			if( Fe2LevN[ipHi][ipLo].ipCont < 1)
				continue;

			/* ots rates, the destp prob was set in hydropesc */
			Fe2LevN[ipHi][ipLo].ots = 
				Fe2LevN[ipHi][ipLo].Aul*
				Fe2LevN[ipHi][ipLo].PopHi*
				Fe2LevN[ipHi][ipLo].Pdest;

			ASSERT( Fe2LevN[ipHi][ipLo].ots >= 0. );

			/* finally dump the ots rate into the stack */
			RT_OTS_AddLine(Fe2LevN[ipHi][ipLo].ots,
				Fe2LevN[ipHi][ipLo].ipCont );
		}
	}

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

	return;
}

/*
 *=====================================================================
 */
/* do outward rates for FeII, 
 * called by RT_diffuse, which is called by cloudy */
void FeII_RTOut(void)
{
	long int ipLo , ipHi;

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

	/* only do this if Fe+ exists */
	if( dense.xIonDense[ipIRON][1] > 0. )
	{
		/* outward line photons */
		for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ )
		{
			for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel; ipHi++ )
			{
				/* >>chng 02 feb 11, skip bogus transitions */
				if( Fe2LevN[ipHi][ipLo].ipCont < 1)
					continue;
				/* >>chng 03 sep 09, change to outline, ouutlin is
				 * not exactly the same in the two - tmn missing in outline */
				outline( &Fe2LevN[ipHi][ipLo] );

			}
		}
	}

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

/*
 *=====================================================================
 */
static void FeIIContCreate(
	/* wavelength of low-lambda end */
	double xLamLow , 
	/* wavelength of high end */
	double xLamHigh , 
	/* number of cells to break this into */
	long int ncell )
{

	double step , wl1;
	long int i;

	/* keep track of whether we have been called - want to be
	 * called a total of one time */
	static int lgCalled=FALSE;

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

	/* return previous number of bands if this is second or later call*/
	if( lgCalled )
	{
#		ifdef DEBUG_FUN
		fputs( " <->FeIIContCreate()\n", debug_fp );
#		endif
		/* return value of number of bands, may be used by calling program*/
		return ;
	}
	lgCalled = TRUE;

	/* how many cells will be needed to go from xLamLow to xLamHigh in ncell steps */
	nFeIIConBins = ncell;

	if( (FeII_Cont = (float **)MALLOC(sizeof(float *)*(unsigned)(nFeIIConBins) ) )==NULL )
	{
		fprintf( ioQQQ, 
			" FeIIContCreate could not MALLOC FeII_Cont 1\n" );
		puts( "[Stop in FeIICreate]" );
		cdEXIT(EXIT_FAILURE);
	}
	/* now make second dim, id wavelength, and lower and upper bounds */
	for( i=0; i<nFeIIConBins; ++i )
	{
		if( (FeII_Cont[i] = (float *)MALLOC(sizeof(float)*(unsigned)(3) ) )==NULL )
		{
			fprintf( ioQQQ, 
				" FeIIContCreate could not MALLOC FeII_Cont 2\n" );
			puts( "[Stop in FeIICreate]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	step = log10( xLamHigh/xLamLow)/ncell;
	wl1 = log10( xLamLow);
	FeII_Cont[0][1] = (float)pow(10. ,wl1);
	for( i=1; i<ncell; ++i)
	{
		/* lower bound of cell */
		FeII_Cont[i][1] = (float)pow(10. , ( wl1 + i*step ) );
	}

	for( i=0; i<(ncell-1); ++i)
	{
		/* upper bound of cell */
		FeII_Cont[i][2] = FeII_Cont[i+1][1];
	}
	FeII_Cont[ncell-1][2] = (float)(pow(10., log10(FeII_Cont[ncell-2][2]) + step) );

	for( i=0; i<ncell; ++i)
	{
		/* wavelength as it will appear in the printout */
		FeII_Cont[i][0] = ( FeII_Cont[i][1] + FeII_Cont[i][2] ) /2.f;
	}
	{
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC )
		{
			FILE *ioDEBUG;
			ioDEBUG = fopen( "fe2cont.txt", "w");
			if( ioDEBUG==NULL)
			{
				fprintf( ioQQQ," fe2con could not open fe2cont.txt for writing\n");
				cdEXIT(EXIT_FAILURE);
			}
			for( i=0; i<ncell; ++i)
			{
				/* wavelength as it will appear in the printout */
				fprintf(ioDEBUG,"%.1f\t%.1f\t%.1f\n",
					FeII_Cont[i][0] , FeII_Cont[i][1] , FeII_Cont[i][2] ) ;
			}
			fclose( ioDEBUG);
		}
	}
	
#	ifdef DEBUG_FUN
	fputs( " <->FeIIContCreate()\n", debug_fp );
#	endif

	return ;
}

/* FeIIOvrLap handle overlapping FeII lines */
static void FeIIOvrLap(void)
{

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

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

	return;
}

/*ParseAtomFeII parse the atom FeII command */
void ParseAtomFeII(char *chCard )
{
	long int i;
	int lgEOL=FALSE;

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

	/* turn on the large verner atom */
	FeII.lgFeIILargeOn = TRUE;
	/* >>chng 05 nov 29, reset number of levels when called, so increased above default of 16 */
	if( lgFeIIMalloc )
	{
		/* space has been allocated, so reset number of levels to whatever it was */
		FeII.nFeIILevel = FeII.nFeIILevelAlloc;
	}
	else
	{
		/* space not allocated yet, set to largest possible number of levels */
		FeII.nFeIILevel = NFE2LEVN;
	}

	/* levels keyword is to adjust number of levels.  But this only has effect
	 * BEFORE space is allocated for the FeII arrays */
	if( lgMatch("LEVE",chCard) )
	{
		/* do option only if space not yet allocated */
		if( !lgFeIIMalloc )
		{
			i = 5;
			/* number of levels for hydrogen at, 2s is this plus one */
			FeII.nFeIILevel = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				/* hit eol with no number - this is a problem */
				NoNumb(chCard);
			}
			if( FeII.nFeIILevel <16 )
			{
				fprintf( ioQQQ, " This would be too few levels, must have at least 16.\n" );
				puts( "[Stop in ParseAtomFeII]" );
				cdEXIT(EXIT_FAILURE);
			}
			else if( FeII.nFeIILevel > NFE2LEVN )
			{
				fprintf( ioQQQ, " This would be too many levels.\n" );
				puts( "[Stop in ParseAtomFeII]" );
				cdEXIT(EXIT_FAILURE);
			}
		}
	}

	/* slow keyword means do not try to avoid evaluating atom */
	else if( lgMatch("SLOW",chCard) )
	{
		FeII.lgSlow = TRUE;
	}

	/* redistribution keyword changes form of redistribution function */
	else if( lgMatch("REDI",chCard) )
	{
		int ipRedis=0;
		/* there are three functions, PRD_, CRD_, and CRDW,
		 * representing partial redistribution, 
		 * complete redistribution (doppler core only, no wings)
		 * and complete with wings */
		/* partial redistribution */
		if( lgMatch(" PRD",chCard) )
		{
			ipRedis = ipPRD;
		}
		/* complete redistribution */
		else if( lgMatch(" CRD",chCard) )
		{
			ipRedis = ipCRD;
		}
		/* complete redistribution with wings */
		else if( lgMatch("CRDW",chCard) )
		{
			ipRedis = ipCRDW;
		}

		/* if not SHOW option (handled below) then we have a problem */
		else if( !lgMatch("SHOW",chCard) )
		{
			fprintf(ioQQQ," There should have been a second keyword on this command.\n");
			fprintf(ioQQQ," Options are _PRD, _CRD, CRDW (_ is space).  Sorry.\n");
			puts( "[Stop in ParseAtomFeII]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* resonance lines */
		if( lgMatch("RESO",chCard) )
		{
			FeII.ipRedisFcnResonance = ipRedis;
		}
		/* subordinate lines */
		else if( lgMatch("SUBO",chCard) )
		{
			FeII.ipRedisFcnSubordinate = ipRedis;
		}
		/* the show option, say what we are assuming */
		else if( lgMatch("SHOW",chCard) )
		{
			fprintf(ioQQQ," FeII resonance lines are ");
			if( FeII.ipRedisFcnResonance ==ipCRDW )
			{
				fprintf(ioQQQ,"complete redistribution with wings\n");
			}
			else if( FeII.ipRedisFcnResonance ==ipCRD )
			{
				fprintf(ioQQQ,"complete redistribution with core only.\n");
			}
			else if( FeII.ipRedisFcnResonance ==ipPRD )
			{
				fprintf(ioQQQ,"partial redistribution.\n");
			}
			else
			{
				fprintf(ioQQQ," PROBLEM Impossible value for ipRedisFcnResonance.\n");
				TotalInsanity();
			}

			fprintf(ioQQQ," FeII subordinate lines are ");
			if( FeII.ipRedisFcnSubordinate ==ipCRDW )
			{
				fprintf(ioQQQ,"complete redistribution with wings\n");
			}
			else if( FeII.ipRedisFcnSubordinate ==ipCRD )
			{
				fprintf(ioQQQ,"complete redistribution with core only.\n");
			}
			else if( FeII.ipRedisFcnSubordinate ==ipPRD )
			{
				fprintf(ioQQQ,"partial redistribution.\n");
			}
			else
			{
				fprintf(ioQQQ," PROBLEM Impossible value for ipRedisFcnSubordinate.\n");
				TotalInsanity();
			}
		}
		else
		{
			fprintf(ioQQQ," here should have been a second keyword on this command.\n");
			fprintf(ioQQQ," Options are RESONANCE, SUBORDINATE.  Sorry.\n");
			puts( "[Stop in ParseAtomFeII]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* trace keyword - print info for each call to FeIILevelPops */
	else if( lgMatch("TRAC",chCard) )
	{
		FeII.lgPrint = TRUE;
	}

	/* only simulate the feii atom, do not actually call it */
	else if( lgMatch("SIMU",chCard) )
	{
		/* option to only simulate calls to FeIILevelPops */
		FeII.lgSimulate = TRUE;
	}

	/* only simulate the feii atom, do not actually call it */
	else if( lgMatch("CONT",chCard) )
	{
		/* the continuum output with the punch feii continuum command */
		i=5;

		/* number of levels for hydrogen at, 2s is this plus one */
		FeII.fe2con_wl1 = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		FeII.fe2con_wl2 = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		FeII.nfe2con = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( lgEOL )
		{
			fprintf(ioQQQ," there are three numbers on the FeII continuum command, start and end wavelengths, and number of intervals.\n");
			/* hit eol with no number - this is a problem */
			NoNumb(chCard);
		}
		/* check that all are positive */
		if( FeII.fe2con_wl1<=0. || FeII.fe2con_wl2<=0. || FeII.nfe2con<= 0 )
		{
			fprintf(ioQQQ," there are three numbers on the FeII continuum command, start and end wavelengths, and number of intervals.\n");
			fprintf(ioQQQ," all three must be greater than zero, sorry.\n");
			puts( "[Stop in ParseAtomFeII]" );
			cdEXIT(EXIT_FAILURE);
		}
		/* make sure that wl1 is less than wl2 */
		if( FeII.fe2con_wl1>= FeII.fe2con_wl2 )
		{
			fprintf(ioQQQ," there are three numbers on the FeII continuum command, start and end wavelengths, and number of intervals.\n");
			fprintf(ioQQQ," the second wavelength must be greater than the first, sorry.\n");
			puts( "[Stop in ParseAtomFeII]" );
			cdEXIT(EXIT_FAILURE);
		}
	}
	/* note no fall-through error since routine can be called with no options,
	 * to turn on the large atom */

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

	return;
}

void PunFeII( FILE * io )
{
	long int n, ipHi ;

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

	for( n=0; n<FeII.nFeIILevel-1; ++n)
	{
		for( ipHi=n+1; ipHi<FeII.nFeIILevel; ++ipHi )
		{
			if( Fe2LevN[ipHi][n].ipCont > 0 ) 
				fprintf(io,"%li\t%li\t%.2e\n", n , ipHi , 
					Fe2LevN[ipHi][n].TauIn );
		}
	}

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

	return;
}

/* include FeII lines in punched optical depths, etc, called from PunchLineStuff */
void FeIIPunchLineStuff( FILE * io , float xLimit  , long index)
{
	long int n, ipHi ;

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

	for( n=0; n<FeII.nFeIILevel-1; ++n)
	{
		for( ipHi=n+1; ipHi<FeII.nFeIILevel; ++ipHi )
		{
			pun1Line( &Fe2LevN[ipHi][n] , io , xLimit  , index , 1. );
		}
	}

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

	return;
}

/* rad pre due to FeII lines called in PresTotCurrent*/
double FeIIRadPress(void)
{

	/* will be used to check on size of opacity, was capped at this value */
	float smallfloat=SMALLFLOAT*10.f;
	double press,
		RadPres1 ;
#	undef	DEBUGFE
#	ifdef DEBUGFE
	double RadMax;
	long ipLoMax , ipHiMax;
#	endif
	long int ipLo, ipHi ;

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

	press = 0.;
#	ifdef DEBUGFE
	RadMax = 0;
	ipLoMax = -1;
	ipHiMax = -1;
#	endif
	/* loop over all levels to find radiation pressure */
	for( ipHi=1; ipHi<FeII.nFeIILevel; ++ipHi )
	{
		for( ipLo=0; ipLo<ipHi; ++ipLo)
		{
			/* >>chng 04 aug 27, add test on ipCont for bogus line */
			if( Fe2LevN[ipHi][ipLo].ipCont > 0  &&
				Fe2LevN[ipHi][ipLo].PopHi > 1e-30 )
			{
				if( Fe2LevN[ipHi][ipLo].PopHi > smallfloat &&
					Fe2LevN[ipHi][ipLo].PopOpc > smallfloat )
				{
					RadPres1 = 5.551e-2*(powi(Fe2LevN[ipHi][ipLo].EnergyWN/
						1.e6,4))*(Fe2LevN[ipHi][ipLo].PopHi/Fe2LevN[ipHi][ipLo].gHi)/
						(Fe2LevN[ipHi][ipLo].PopLo/Fe2LevN[ipHi][ipLo].gLo)*
						RT_LineWidth(&Fe2LevN[ipHi][ipLo]);
#					ifdef DEBUGFE
					if( RadPres1 > RadMax )
					{
						RadMax = RadPres1;
						ipLoMax = ipLo;
						ipHiMax = ipHi;
					}
#					endif
					press += RadPres1;
				}
			}
		}
	}

#	ifdef	DEBUGFE
	/* option to print radiation pressure */
	if( iteration > 1 || nzone > 1558 )
	{
		fprintf(ioQQQ,"DEBUG FeIIRadPress finds P(FeII) = %.2e %.2e %li %li width %.2e\n",
			press  ,
			RadMax ,
			ipLoMax ,
			ipHiMax ,
			RT_LineWidth(&Fe2LevN[9][0]) );
		DumpLine( &Fe2LevN[9][0] );
	}
#	endif
#	undef	DEBUGFE

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

	return press;
}

/* internal energy of FeII called in PresTotCurrent */
double FeII_InterEnergy(void)
{
	double energy;
	long int n, ipHi ;

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

	energy = 0.;
	for( n=0; n<FeII.nFeIILevel-1; ++n)
	{
		for( ipHi=n+1; ipHi<FeII.nFeIILevel; ++ipHi )
		{
			if( Fe2LevN[ipHi][n].PopHi > 1e-30 )
			{
				energy += 
					Fe2LevN[ipHi][n].PopHi* Fe2LevN[ipHi][n].EnergyErg;
			}
		}
	}

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

	return energy;
}

/* HP cc cannot compile this routine with any optimization */
#if defined(__HP_cc)
#	pragma OPT_LEVEL 1
#endif
/*FeIILyaPump find rate of Lya excitation of the FeII atom */
static void FeIILyaPump(void)
{

	long int ipLo ,
		ipHi;
	double EnerLyaProf2, 
	  EnerLyaProf3, 
	  EnergyWN,
	  Gup_ov_Glo, 
	  PhotOccNum_at_nu, 
	  PumpRate, 
	  de, 
	  FeIILineWidthHz, 
	  taux;

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

	/* lgLyaPumpOn is false if no Lya pumping, with no feii command */
	/* get rates FeII atom is pumped */

	/* elsewhere in this file the dest prob hydro.dstfe2lya is defined from
	 * quantites derived here, and the resulting populations */
	if( FeII.lgLyaPumpOn )
	{

		/*************trapeze form La profile:de,EnerLyaProf1,EnerLyaProf2,EnerLyaProf3,EnerLyaProf4*************************
		 * */
		/* width of Lya in cm^-1 */
		/* HLineWidth has units of cm/s, as was evaluated in PresTotCurrent */
		/* the factor is 1/2 of E(Lya, cm^-1_/c */
		de = 1.37194e-06*hydro.HLineWidth*2.0/3.0;
		/* 82259 is energy of Lya in wavenumbers, so these are the form of the trapezoid */
		EnerLyaProf1 = 82259.0 - de*2.0;
		EnerLyaProf2 = 82259.0 - de;
		EnerLyaProf3 = 82259.0 + de;
		EnerLyaProf4 = 82259.0 + de*2.0;

		/* find Lya photon occupation number */
		if( EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].PopHi > SMALLFLOAT )
		{
			/* This is the photon occupation number at the Lya line center */
			PhotOccNumLyaCenter = 
				MAX2(0.,1.0- EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pelec_esc - 
				EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pesc)/
			  (EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].PopLo/EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].PopHi*3. - 1.0);
		}
		else
		{
			/* lya excitation temperature not available */
			PhotOccNumLyaCenter = 0.;
		}
	}
	else
	{
		PhotOccNumLyaCenter = 0.;
		de = 0.;
		EnerLyaProf1 = 0.;
		EnerLyaProf2 = 0.;
		EnerLyaProf3 = 0.;
		EnerLyaProf4 = 0.;
	}

	/* is Lya pumping enabled?  */
	if( FeII.lgLyaPumpOn )
	{
		for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ )
		{
			for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel; ipHi++ )
			{
				/* on first iteration optical depth in line is inward only, on later
				 * iterations is total optical depth */
				if( iteration == 1 )
				{
					taux = Fe2LevN[ipHi][ipLo].TauIn;
				}
				else
				{
					taux = Fe2LevN[ipHi][ipLo].TauTot;
				}

				/* Gup_ov_Glo is ratio of g values */
				Gup_ov_Glo = Fe2LevN[ipHi][ipLo].gHi/Fe2LevN[ipHi][ipLo].gLo;

				/* the energy of the FeII line */
				EnergyWN = Fe2LevN[ipHi][ipLo].EnergyWN;

				if( EnergyWN >= EnerLyaProf1 && EnergyWN <= EnerLyaProf4  &&  taux > 1 )
				{
					/* this branch, line is within the Lya profile */

					/*
					 * Lya source function, at peak is PhotOccNumLyaCenter,
					 *
					 *     Prof2    Prof3
					 *       ----------
					 *      /          \
					 *     /            \
					 *    /              \
					 *  ======================
					 * Prof1              Prof4
					 *
					 */

					if( EnergyWN < EnerLyaProf2 )
					{
						/* linear interpolation on edge of trapazoid */
						PhotOccNum_at_nu = PhotOccNumLyaCenter*(EnergyWN - EnerLyaProf1)/ de;
					}
					else if( EnergyWN < EnerLyaProf3 )
					{
						/* this is the central plateau */
						PhotOccNum_at_nu = PhotOccNumLyaCenter;
					}
					else
					{
						/* linear interpolation on edge of trapazoid */
						PhotOccNum_at_nu = PhotOccNumLyaCenter*(EnerLyaProf4 - EnergyWN)/de;
					}

					/* at this point Lya source function at FeII line energy is defined, but
					 * we need to multiply by line width in Hz,
					 * >>refer	fe2	pump rate	Netzer, H., Elitzur, M., & Ferland, G.J., 1985, ApJ, 299, 752-762*/

					/* width of feii line in Hz  */
					FeIILineWidthHz = 1.e8/(EnergyWN*299792.5)*sqrt(log(taux))*DoppVel.doppler[ipIRON];

					/* final Lya pumping rate, s^-1*/
					PumpRate = FeIILineWidthHz * PhotOccNum_at_nu * Fe2LevN[ipHi][ipLo].Aul*
					  powi(82259.0f/EnergyWN,3);
					/* above must be bogus, use just occ num times A */
					PumpRate = Fe2LevN[ipHi][ipLo].Aul* PhotOccNum_at_nu ;

					/* Lya pumping rate from ipHi to lower n */
					Fe2LPump[ipHi][ipLo] += PumpRate;

					/* Lya pumping rate from n to upper ipHi */
					Fe2LPump[ipLo][ipHi] += PumpRate*Gup_ov_Glo;
				}
			}
		}
	}
#	ifdef DEBUG_FUN
	fputs( " <->FeIILyaPump()\n", debug_fp );
#	endif

	return;
}

/* end work around bugs in HP compiler */
#if defined(__HP_cc)
#pragma OPTIMIZE OFF
#pragma OPTIMIZE ON
#endif

