/* This file is part of Cloudy and is copyright (C) 1978-2004 by Gary J. Ferland.
 * For conditions of distribution and use, see copyright notice in license.txt */
/*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  */
/*AssertFeIIDep called by assert feii depart coef command */
/*FeIIPrint print feii information */
/*FeIICollStrength make up collision data for feii */
/*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 */
/*FeIITauAver 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*/
/*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(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 */
/* evaluate collision strenths, both interpolating on r-mat and creating g-bar */
void FeIICollStrength(void);
/* find rate of Lya excitation of the FeII atom */
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 Fe2Bands.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 float **Fe2A;

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

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

/* Fe2DepCoef[NFE2LEVN];
float cli[NFEIILINES], 
  cfe[NFEIILINES];*/
static double 
	/* departure coefficients */
	*Fe2DepCoef , 
	/* level populations */
	*Fe2LevelPop ,
	/* 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;

/*
 *=====================================================================
 */
/* 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 */
	nFeIILevelAlloc = nFeIILevel;

	/* set up array to save FeII transition probabilities */
	if( (Fe2A = (float **)MALLOC(sizeof(float *)*(unsigned long)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 < nFeIILevel; ++ipHi )
	{
		if( (Fe2A[ipHi]=(float*)MALLOC(sizeof(float )*(unsigned long)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 = (float **)MALLOC(sizeof(float *)*(unsigned long)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 = (float **)MALLOC(sizeof(float *)*(unsigned long)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 < nFeIILevel; ++ipHi )
	{
		if( (Fe2CPump[ipHi]=(float*)MALLOC(sizeof(float )*(unsigned long)nFeIILevel)) ==NULL)
		{
			fprintf( ioQQQ, 
				" FeIICreate could not MALLOC Fe2CPump 2\n" );
			puts( "[Stop in FeIICreate]" );
			cdEXIT(EXIT_FAILURE);
		}
		if( (Fe2LPump[ipHi]=(float*)MALLOC(sizeof(float )*(unsigned long)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)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 < nFeIILevel; ++ipHi )
	{
		if( (Fe2Coll[ipHi]=(float*)MALLOC(sizeof(float )*(unsigned long)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)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<nFeIILevel; ++i )
	{
		if( (xMatrix[i] = (double *)MALLOC(sizeof(double)*(unsigned long)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)(nFeIILevel*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)(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)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 < 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)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)nFeIILevel*sizeof(double) )) == NULL )
	{
		fprintf(ioQQQ," FeIICreate could not MALLOC Fe2LevelPop\n");
		puts( "[Stop in FeIICreate]" );
		cdEXIT(EXIT_FAILURE);
	}
	/*FeIIBoltzmann[NFE2LEVN];*/
	if( ( FeIIBoltzmann = (double*)MALLOC( (unsigned long)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)nFeIILevel)) ==NULL)
		BadMalloc();
	if( (ncs1=(int**)MALLOC(sizeof(int *)*(unsigned long)nFeIILevel)) ==NULL)
		BadMalloc();

	for( ipHi=1; ipHi < 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 < (nFeIILevel - 1); ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi < 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 < 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 < (nFeIILevel - 1); ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi < 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]));

			/* 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 (fcool) and output fecool from Cloudy 01.29.97god
 *********and cooling derivative (dfcool)
 ************ 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:FeIICollStrength;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<nFeIILevel; ++n)
	{
		for( ipLo=0; ipLo<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 */
	FeIICollStrength();

	/* pumping and spontantous decays */
	for( n=0; n<nFeIILevel; ++n)
	{
		for( ipHi=n+1; ipHi<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<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<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<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<nFeIILevel; ++n)
			{
				/*fprintf(ioQQQ,"\n");*/
				/* now print the matrix*/
				for( ipLo=0; ipLo<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 < 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_)*nFeIILevel+(J_)))

	/* copy current contents of xMatrix array over to special amat array*/
	for( ipHi=0; ipHi < nFeIILevel; ipHi++ )
	{
		for( i=0; i < 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(nFeIILevel, nFeIILevel, amat, nFeIILevel, ipiv, &info);
	getrs_wrapper('N', nFeIILevel, 1, amat, nFeIILevel, ipiv, yVector, nFeIILevel, &info);

	if( info != 0 )
	{
		fprintf( ioQQQ, " FeIILevelPops: dgetrs finds singular or ill-conditioned matrix\n" );
		puts( "[Stop in CO_step]" );
		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 < nFeIILevel; ipLo++ )
	{
		if(yVector[ipLo] <= 0. )
		{
			lgPopNeg = TRUE;
			fprintf(ioQQQ," 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," TROUBLE FeIILevelPops finds non-positive level population\n" );
	}

	/* now set line opacities, intensities, and level populations*/
	for( ipLo=0; ipLo < (nFeIILevel - 1); ipLo++ )
	{
		for( ipHi=ipLo+1; ipHi < 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] ) );
		}
	}

	/* the hydrogen Lya destruction rate, then probability */
	hydro.dstfe2lya = 0.;
	EnergyWN = 0.;
	/* count how many photons were removed-added */
	for( ipLo=0; ipLo < (nFeIILevel - 1); ipLo++ )
	{
		for( ipHi=ipLo+1; ipHi < 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.;
	}

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

	/* 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 < 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 */
	for( i=0; i < 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.fcool = 0.0f;
	/* this is be heating, not heating minus cooling */
	FeII.feheat = 0.f;
	/* derivative of cooling */
	FeII.dfcool = 0.0f;

	/* compute heating and cooling due to model atom */
	for( ipLo=0; ipLo < (nFeIILevel - 1); ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi < 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.fcool += (float)MAX2(0., OneLine);

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

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

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

/*
 *=====================================================================
 */
/*FeIICollStrength calculates collisional strenths for all lines that have no data */
void FeIICollStrength(void)
{
	/* this will be used to only reevaluate cs when temp changes */
	/* >>chng 00 jun 02, demoted next to float, PvH */
	static float OldTemp = -1.f;
	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( "<+>FeIICollStrength()\n", debug_fp );
#	endif

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

	for( ipLo=0; ipLo < (nFeIILevel - 1); ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi < 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,">>>impossible ncs1 in FeIILevelPops\n");
			}

			/*y = 1.438826f*Fe2LevN[ipHi][ipLo].EnergyWN/ phycon.te;*/
			y = Fe2LevN[ipHi][ipLo].EnergyWN/ 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 = (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 < 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 < nFeIILevel; ipHi++ )
	{
		/*FeIIBoltzmann[ipHi] = (float)sexp( 1.438826f*Fe2LevN[ipHi][0].EnergyWN/phycon.te );*/
		/* >>chng 99 may 21, from above to following slight different number from phyconst.h */
		FeIIBoltzmann[ipHi] = sexp( Fe2LevN[ipHi][0].EnergyWN/phycon.te_wn );
	}

	/* now possibly trim down atom if Boltzmann factors for upper levels are zero */
	ipHi = 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 */
	if( ipHi <= 1 )
	{
		/* this is basically a sanity check and can't possibly happen */
		fprintf( ioQQQ, " FeIILevelPops trimmed boltzmann factors too low.\n" );
		puts( "[Stop in FeIILevelPops]" );
		cdEXIT(EXIT_FAILURE);
	}
	else
	{
		/* is nearly all cases this does nothing since ipHi and nFeIILevel
		 * are equal . . . */
		nFeIILevel = ipHi+1;
		/*fprintf(ioQQQ," levels reset to %li\n",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<nFeIILevel; ++ipLo)
	{
		for( ipHi=ipLo+1; ipHi<nFeIILevel; ++ipHi )
		{
			/* collisional deexcitation rate coefficient from ipHi to lower ipLo
			 * note that it needs eden to become rate*/
			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 */
			Fe2Coll[ipLo][ipHi] = (float)(Fe2Coll[ipHi][ipLo]*FeIIBoltzmann[ipHi]/FeIIBoltzmann[ipLo]*
				Fe2LevN[ipHi][ipLo].gHi/Fe2LevN[ipHi][ipLo].gLo);
		}
	}

#	ifdef DEBUG_FUN
	fputs( " <->FeIICollStrength()\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 */
double FeIISumBand(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 < (nFeIILevel - 1); ipLo++ )
		{
			for( ipHi=ipLo + 1; ipHi < 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 < (nFeIILevel - 1); ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi < 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;
}

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

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

	/* called at end of iteration */
	for( ipLo=0; ipLo < nFeIILevel-1; ipLo++ )
	{
		for( ipHi=ipLo+1; ipHi < nFeIILevel; ipHi++ )
		{
			RT_line_one_tau_reset( &Fe2LevN[ipHi][ipLo],0.5 );
		}
	}

	/* call the line overlap routine */
	FeIIOvrLap();

#	ifdef DEBUG_FUN
	fputs( " <->FeIITauAver()\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 < nFeIILevel-1; ipLo++ )
	{
		for( ipHi=ipLo+1; ipHi < 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 < (nFeIILevel - 1); ipLo++ )
	{
		for( ipHi=ipLo+1; ipHi < 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");
	}

	if( !FeII.lgFeIION )
		return;

	/* this routine drives calls to make RT relations with large FeII atom */
	for( ipLo=0; ipLo < (nFeIILevel - 1); ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi < 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 < (nFeIILevel - 1); ipLo++ )
		{
			for( ipHi=ipLo + 1; ipHi < 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 < (nFeIILevel - 1); ipLo++ )
		{
			for( ipHi=ipLo + 1; ipHi < 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;
}

/*
 *=====================================================================
 */
/*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 < (nFeIILevel - 1); ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi < 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 < (nFeIILevel - 1); ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi < 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;
}

#if 0
/*
 *=====================================================================
 */
/*FeIIEmitOut add large FeII emission to outward beam - called once per zone in RT_diffuse */
void FeIIEmitOut(double VolFac, 
  double ref)
{
	long int ipHi, 
	  ipLo, 
	  ip;

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

	if( FeII.lgFeIION || dense.xIonDense[ipIRON][1] > 0. )
	{
#		ifdef DEBUG_FUN
		fputs( " <->FeIIEmitOut()\n", debug_fp );
#		endif
		return;
	}

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

			/* pointer to line energy in continuum array */
			ip = Fe2LevN[ipHi][ipLo].ipCont;
			rfield.outlin[ip-1] += 
			  (float)(Fe2LevN[ipHi][ipLo].phots*
			  VolFac*opac.tmn[ip-1]*Fe2LevN[ipHi][ipLo].ColOvTot);
			rfield.reflin[ip-1] += 
				(float)(Fe2LevN[ipHi][ipLo].phots* ref);
		}
	}

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

/*
 *=====================================================================
 */
/*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

	if( FeII.lgFeIION )
	{
		/* 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 < (nFeIILevel - 1); ipLo++ )
		{
			for( ipHi=ipLo + 1; ipHi < 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 in routine zero and it
	 * zero's out various elements of the FeII arrays
	 * it is called on every iteration
	 * */
	for( ipLo=0; ipLo < (nFeIILevel - 1); ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi < nFeIILevel; ipHi++ )
		{

			/* 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.Fe2L16Tot = 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( 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( 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;
		/* false, only punch subset in above init */
		/* first 64 do all lines */
		for( ipHi=1; ipHi<64; ++ipHi )
		{
			for( ipLo=0; ipLo<ipHi; ++ipLo )
			{
				fprintf(ioPUN,"%4li%4li ",ipLo,ipHi );
				Punch1LineData( &Fe2LevN[ipHi][ipLo] , ioPUN);
			}
		}
		fprintf( ioPUN , "\n");

		/* higher than 64 only do real transitions */
		for( ipHi=64; ipHi<nFeIILevel; ++ipHi )
		{
			for( ipLo=64; 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<=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 );
	ASSERT( ioPUN != NULL );

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

#	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 */
	nFeIILevel = 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.fcool = 0.;
	FeII.dfcool = 0.;
	FeII.feheat = 0.;

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

	/* will not compute large feii atom */
	FeII.lgFeIION = 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 */
		nFeIILevel = nFeIILevelAlloc;
	}
	else
	{
		/* space not allocated yet, set to largest possible number of levels */
		nFeIILevel = NFE2LEVN;
	}

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

/*
 *=====================================================================
 */
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<=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 );
	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 <=*/
	/*lint -e668 passing null pointer */
	if( nPUN <= nFeIILevel )
	{
		fprintf( ioPUN , "%e ",Fe2LevelPop[nPUN-1] );
	}
	else
	{
		fprintf( ioPUN , "%e ",0. );
	}
	/*lint +e668 passing null pointer */

#	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 */
	 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 , "Fe2Bands.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 ||  !FeII.lgFeIION)
	{

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

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

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

	/* find sum of deviation of departure coef from unity */
	*BigError = 0.;
	*pred = 0.;
	sum2 = 0;
	for( n=0; n<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)nFeIILevel ;
	ASSERT( (arg >= 0.) );
	*StdDev = sqrt( arg / (double)(nFeIILevel - 1.) );

	/* this is average value, should be unity */
	*pred /= (double)(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

	if( !FeII.lgFeIION )
	{
#		ifdef DEBUG_FUN
		fputs( " <->FeII_OTS()\n", debug_fp );
#		endif
		return;
	}

	for( ipLo=0; ipLo < (nFeIILevel - 1); ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi < 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 = (float)(
				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( FeII.lgFeIION && dense.xIonDense[ipIRON][1] > 0. )
	{
		/* outward line photons */
		for( ipLo=0; ipLo < (nFeIILevel - 1); ipLo++ )
		{
			for( ipHi=ipLo + 1; ipHi < 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 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.lgFeIION = TRUE;

	/* 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 */
			nFeIILevel = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				/* hit eol with no number - this is a problem */
				NoNumb(chCard);
			}
			if( 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( 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);
		}
	}

	/* print keyword print comment for each call to FeIILevelPops */
	else if( lgMatch("PRIN",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 ;
	for( n=0; n<nFeIILevel-1; ++n)
	{
		for( ipHi=n+1; ipHi<nFeIILevel; ++ipHi )
		{
			if( Fe2LevN[ipHi][n].ipCont > 0 ) 
				fprintf(io,"%li\t%li\t%.2e\n", n , ipHi , 
					Fe2LevN[ipHi][n].TauIn );
		}
	}
}

/* include FeII lines in punched optical depths, etc, called from PunchLineStuff */
void FeIIPunchLineStuff( FILE * io , float xLimit  , long index)
{
	long int n, ipHi ;
	for( n=0; n<nFeIILevel-1; ++n)
	{
		for( ipHi=n+1; ipHi<nFeIILevel; ++ipHi )
		{
			pun1Line( &Fe2LevN[ipHi][n] , io , xLimit  , index);
		}
	}
}

/* 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;
	long int n, ipHi ;

	press = 0.;
	if( FeII.lgFeIION )
	{
		for( n=0; n<nFeIILevel-1; ++n)
		{
			for( ipHi=n+1; ipHi<nFeIILevel; ++ipHi )
			{
				/* >>chng 04 aug 27, add test on ipCont for bogus line */
				if( Fe2LevN[ipHi][n].ipCont > 0  &&
					Fe2LevN[ipHi][n].PopHi > 1e-30 )
				{
					if( Fe2LevN[ipHi][n].PopHi > smallfloat &&
						Fe2LevN[ipHi][n].PopOpc > smallfloat )
					{
						double RadPres1 = 5.551e-2*(powi(Fe2LevN[ipHi][n].EnergyWN/
							1.e6,4))*(Fe2LevN[ipHi][n].PopHi/Fe2LevN[ipHi][n].gHi)/
							(Fe2LevN[ipHi][n].PopLo/Fe2LevN[ipHi][n].gLo)*
							RT_LineWidth(&Fe2LevN[ipHi][n]);
							press += RadPres1;
					}
				}
			}
		}
	}
	return press;
}

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

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

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

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

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

		/* gary change, this crashed with texclya = 0 at start */
		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 < (nFeIILevel - 1); ipLo++ )
		{
			for( ipHi=ipLo + 1; ipHi < 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] += (float)PumpRate;

					/* Lya pumping rate from n to upper ipHi */
					Fe2LPump[ipLo][ipHi] += (float)(PumpRate*Gup_ov_Glo);
				}
			}
		}
	}
}
#if defined(__HP_cc)
#pragma OPTIMIZE OFF
#pragma OPTIMIZE ON
#endif
