/*CoolIron compute iron cooling */
#include <string.h>
#include "cddefines.h"
#include "coolheavy.h"
#include "poplevls.h"
#include "resion.h"
#include "taulines.h"
#include "tepowers.h"
#include "phycon.h"
#include "h.h"
#include "logte.h"
#include "fe2neg.h"
#include "heating.h"
#include "converge.h"
#include "holod.h"
#include "trace.h"
#include "tcool.h"
#include "hydrogenic.h"
#include "ionfracs.h"
#include "leveln.h"
#include "ligbar.h"
#include "sexp.h"
#include "makecs.h"
#include "coladd.h"
#include "level2.h"
#include "putcs.h"
#include "beseq.h"
#include "zonecnt.h"
#include "level3.h"
#include "pop3.h"
#include "coolmetals.h"
#include "fe14cs.h"
#include "fe4lev12.h"
#include "popexc.h"
#include "fe2lev16.h"
#include "pop371.h"
#define	NLFE2	6

void CoolIron()
{
	int lgNegPop;
	long int i, 
	  ipdest[NLFE2][NLFE2], 
	  j;
	double cs ,
	  dFe2dt, 
	  pump[NLFE2][NLFE2],
	  data[NLFE2][NLFE2], 
	  dest[NLFE2][NLFE2], 
	  popFe2[NLFE2], 
	  tused,
	  cs2s2p, 
	  cs2s3p;
	float 
	  p2,
	  PumpLyaFeII;

	/* stat weights for Fred's 6 level model FeII atom */
	static double gFe2[NLFE2]={1.,1.,1.,1.,1.,1.};
	/* excitation energies (Kelvin) for Fred's 6 level model FeII atom */
	static double ex[NLFE2]={0.,3.32e4,5.68e4,6.95e4,1.15e5,1.31e5};

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

	/* series of FeI lines from Dima Verner's list, each 2-lev atom
	 *
	 * Fe I 3884 */
	MakeCS(&TauLines[ipFeI3884]);
	level2(&TauLines[ipFeI3884]);

	/* Fe I 3729 */
	MakeCS(&TauLines[ipFeI3729]);
	level2(&TauLines[ipFeI3729]);

	/* Fe I 3457 */
	MakeCS(&TauLines[ipFeI3457]);
	level2(&TauLines[ipFeI3457]);

	/* Fe I 3021 */
	MakeCS(&TauLines[ipFeI3021]);
	level2(&TauLines[ipFeI3021]);

	/* Fe I 2966 */
	MakeCS(&TauLines[ipFeI2966]);
	level2(&TauLines[ipFeI2966]);


	if( FeII.lgFeIION )
	{
		/* this is the Verner's 371 level FeII atom
		 *666    error! evaluate this destruction probability!! */
		hydro.dstfe2lya = 0.;
		/* >>chng 97 jan 17, evaluate large FeII atom cooling every time temp changes
		 * >>chng 97 jan 31, added check on zone since must reeval even const temp
		 * >>chng 99 May 21, reeval when zone or temperature changes, but not when
		 * abundance changes, since we can simply rescale cooling */

		/* only redo large feii when redoing everything */
		if( xIonFracs[2][25] == 0. )
		{
			/* zero abundance, do nothing */
			holod.dfcool = 0.;
			holod.fcool = 0.;
			holod.feheat = 0.;
		}
		else
		{
			/* these are used to only evaluated feii one time per temperature, zone, and abundance */
			static float TUsed = 0. , 
				AbunUsed = 0.;
			/* remember which zone last evaluated with */
			static long int nZUsed=-1,
				/* make sure at least two calls per zone */
				nCall=0;;

			/* totally reevaluate large atom if new zone, or cooling is significant
			 * and temperture changed, we are in search phase,
			 * lgSlow option set true with atom feii slow, forces constant
			 * evaluation of atom */
			if( FeII.lgSlow || 
				conv.lgSearch ||
				ZoneCnt.nzone !=nZUsed || 
				/* on new call, nCall is now set at previous zone's number of calls.
				 * it is set to zero below, so on second call, nCall is 0.  On 
				 * third call nCall is 1.  Check for <1 is to insure at least two calls */
				nCall < 1 ||
				( phycon.te != TUsed &&  holod.fcool/tcool.ctot> 0.002 &&  
				fabs(xIonFracs[2][25]-AbunUsed)/MAX2(AbunUsed,SMALLFLOAT)> 0.002 ) )
			{

				if( nZUsed == ZoneCnt.nzone )
				{
					/* not first call, increment, check above to make sure at least
					 * two evaluations */
					++nCall;
				}
				else
				{
					/* first call this zone set nCall to zero*/
					nCall = 0;
				}

				/* option to trace convergence and feii calls */
				if( trace.lgTrConvg>4 )
				{
					fprintf( ioQQQ, "       CoolIron calling pop371 since ");
					if( phycon.te!= TUsed )
					{
						fprintf( ioQQQ, 
							"temperature changed, old new are %g %g, nCall %li\n", 
							TUsed, phycon.te , nCall);
					}
					else if( ZoneCnt.nzone !=nZUsed )
					{
						fprintf( ioQQQ, 
							"new zone, nCall %li\n", nCall );
					}
					else if( nCall < 2 )
					{
						fprintf( ioQQQ, 
							"not second nCall %li\n" , nCall );
					}
					else if( phycon.te!= TUsed &&  holod.fcool/tcool.ctot> 0.001 )
					{
						fprintf( ioQQQ, 
							"temp or cooling changed, new are %g %g nCall %li\n", 
							phycon.te, holod.fcool, nCall );
					}
				}

				/* remember parameters for current conditions */
				TUsed = phycon.te;
				AbunUsed = xIonFracs[2][25];
				nZUsed = ZoneCnt.nzone;

				/* this print turned on with atom feii print command */
				if( FeII.lgPrint)
				{
					fprintf(ioQQQ,
						" pop371 called zone %4li te %5f abun %10e c(fe/tot):%6f nCall %li\n", 
						ZoneCnt.nzone,phycon.te,AbunUsed,holod.fcool/tcool.ctot,nCall);
				}

				/* FeII.lgSimulate was set true with simulate flag on atom feii command,
				 * for bebugging without actually calling the routine */
				if( !FeII.lgSimulate)
				{
					pop371();
				}
				/*fprintf(ioQQQ," pop371 called\n");*/

			}
			else if( xIonFracs[2][25]!=AbunUsed )
			{
				float ratio;
				/* this branch, same zone and temperature, but small change in abundance, so just
				 * rescale cooling and deriv by this change.  assumption is that very small changes
				 * in abundance occurs as ots rates damp out */
				if( trace.lgTrConvg>4 )
				{
					fprintf( ioQQQ, "       CoolIron rescaling pop371 since small abund change\n");
				}
				ratio = xIonFracs[2][25]/AbunUsed;
				holod.fcool *= ratio;
				holod.dfcool *= ratio;
				holod.feheat *= ratio;
				AbunUsed = xIonFracs[2][25];
			}
			else
			{
				/* this is case where temp is unchanged, so heating and cooling same too */
				if( trace.lgTrConvg>4 )
				{
					fprintf( ioQQQ, "       CoolIron NOT calling pop371\n");
				}
			}

			/* evaluate some strong lines that would have been evaluated by the 16 level atom */
			FeIIFillLow16();

			/* now update heating/cooling */
			/* total cooling */
			CoolHeavy.Fecool = holod.fcool;

			/* counts as heating derivative if negative cooling */
			if( CoolHeavy.Fecool > 0. )
			{
				tcool.dCooldT += holod.dfcool;
			}

			coladd("Fe 2",0,MAX2(0.,holod.fcool));

			/* add negative cooling to heating stack */
			HeatingCom.heating[27][25] = MAX2(0.,-holod.fcool);

			if( trace.lgTrace && trace.lgCoolTr )
			{
				fprintf( ioQQQ, " Large FeII returns te, cooling, dc=%11.3e%11.3e%11.3e\n", 
				  phycon.te, CoolHeavy.Fecool, holod.dfcool );
			}
		}
	}
	else
	{
		holod.feheat = 0.;

		/* following treatment of Fe II follows
		 * >>refer Wills, B.J., Wills, D., Netzer, H. 1985, ApJ, 288, 143
		 * all elements are used, and must be set to zero if zero */

		/*zero out all arrays, then check that upper diag remains zero below */
		for( i=0; i < NLFE2; i++ )
		{
			for( j=0; j < NLFE2; j++ )
			{
				data[j][i] = 0.;
				dest[j][i] = 0.;
				pump[j][i] = 0.;
				ipdest[j][i] = 0;
			}
		}

		/* now put in real data for lines */
		data[0][1] = 1.;
		data[0][2] = TauLines[ipTuv3].Pesc*TauLines[ipTuv3].Aul;
		dest[0][2] = TauLines[ipTuv3].Pdest*TauLines[ipTuv3].Aul;
		pump[0][2] = TauLines[ipTuv3].pump;
		ipdest[0][2] = TauLines[ipTuv3].ipCont;

		data[0][5] = TauLines[ipTFe16].Pesc*TauLines[ipTFe16].Aul;
		dest[0][5] = TauLines[ipTFe16].Pdest*TauLines[ipTFe16].Aul;
		/* continuum pumping of n=6 */
		pump[0][5] = TauLines[ipTFe16].pump;
		/* Ly-alpha pumping */
		if( xIonFracs[2][25] > 0. )
		{
			PumpLyaFeII = hydro.hn[0][IP2P]*h.hii*
				HydroLines[0][IP2P][IP1S].Aul*
			  hydro.dstfe2lya/xIonFracs[2][25];
			pump[0][5] += PumpLyaFeII;
		}
		else
		{
			PumpLyaFeII = 0.;
		}

		ipdest[0][5] = TauLines[ipTFe16].ipCont;

		data[1][2] = TauLines[ipTr48].Pesc*TauLines[ipTr48].Aul;
		dest[1][2] = TauLines[ipTr48].Pdest*TauLines[ipTr48].Aul;
		pump[1][2] = TauLines[ipTr48].pump;
		ipdest[1][2] = TauLines[ipTr48].ipCont;

		data[1][5] = TauLines[ipTFe26].Pesc*TauLines[ipTFe26].Aul;
		dest[1][5] = TauLines[ipTFe26].Pdest*TauLines[ipTFe26].Aul;
		pump[1][5] = TauLines[ipTFe26].pump;
		ipdest[1][5] = TauLines[ipTFe26].ipCont;

		data[2][3] = TauLines[ipTFe34].Pesc*TauLines[ipTFe34].Aul;
		dest[2][3] = TauLines[ipTFe34].Pdest*TauLines[ipTFe34].Aul;
		pump[2][3] = TauLines[ipTFe34].pump;
		ipdest[2][3] = TauLines[ipTFe34].ipCont;

		data[2][4] = TauLines[ipTFe35].Pesc*TauLines[ipTFe35].Aul;
		dest[2][4] = TauLines[ipTFe35].Pdest*TauLines[ipTFe35].Aul;
		pump[2][4] = TauLines[ipTFe35].pump;
		ipdest[2][4] = TauLines[ipTFe35].ipCont;

		data[3][5] = TauLines[ipTFe46].Pesc*TauLines[ipTFe46].Aul;
		dest[3][5] = TauLines[ipTFe46].Pdest*TauLines[ipTFe46].Aul;
		pump[3][5] = TauLines[ipTFe46].pump;
		ipdest[3][5] = TauLines[ipTFe46].ipCont;

		data[4][5] = TauLines[ipTFe56].Pesc*TauLines[ipTFe56].Aul;
		dest[4][5] = TauLines[ipTFe56].Pdest*TauLines[ipTFe56].Aul;
		pump[4][5] = TauLines[ipTFe56].pump;
		ipdest[4][5] = TauLines[ipTFe56].ipCont;

		/* these are collision strengths */
		data[1][0] = 1.;
		data[2][0] = 12.;
		data[3][0] = 1.;
		data[4][0] = 1.;
		data[5][0] = 12.;
		data[2][1] = 6.;
		data[3][1] = 1.;
		data[4][1] = 1.;
		data[5][1] = 12.;
		data[3][2] = 6.;
		data[4][2] = 12.;
		data[5][2] = 1.;
		data[4][3] = 1.;
		data[5][3] = 12.;
		data[5][4] = 6.;

		/*void LevelN(long,long,float,double[],double[],double[],double*,
		double*,double*,long*,float*,float*,STRING,int*);*/
		LevelN(NLFE2,
			NLFE2,
			xIonFracs[2][25],
			gFe2,
		  ex,
		  popFe2,
		  (double*)data,
		  (double*)dest,
		  (double*)pump,
		  (long*)ipdest,
		  &CoolHeavy.Fecool,
		  &dFe2dt,
		  "FeII",
		  &lgNegPop);

		if( lgNegPop )
			Fe2Neg.nFe2Neg += 1;

		coladd("Fe 2",0,MAX2(0.,CoolHeavy.Fecool));
		HeatingCom.heating[27][25] = MAX2(0.,-CoolHeavy.Fecool);
		tcool.dCooldT += dFe2dt;
		TauLines[ipTuv3].PopLo = (float)popFe2[0];
		TauLines[ipTuv3].PopOpc = (float)(popFe2[0] - popFe2[2]);
		resion.emuv3 = (float)(popFe2[2]*data[0][2]*TauLines[ipTuv3].EnergyErg);
		TauLines[ipTuv3].phots = (float)(popFe2[2]*data[0][2]);
		resion.otsfe2 = (float)(popFe2[2]*TauLines[ipTuv3].Pdest*TauLines[ipTuv3].Aul);

		TauLines[ipTr48].PopLo = (float)popFe2[1];
		TauLines[ipTr48].PopOpc = (float)(popFe2[1] - popFe2[2]);
		CoolHeavy.rmt48 = popFe2[2]*data[1][2]*TauLines[ipTr48].EnergyErg;
		TauLines[ipTr48].phots = (float)(popFe2[2]*data[1][2]);

		CoolHeavy.for7 = popFe2[1]*data[0][1]*4.65e-12;

		TauLines[ipTFe16].PopLo = (float)popFe2[0];
		TauLines[ipTFe16].PopOpc = (float)(popFe2[0] - popFe2[5]);
		TauLines[ipTFe16].phots = (float)(popFe2[5]*data[0][5]);
		CoolHeavy.cFe16 = TauLines[ipTFe16].phots*TauLines[ipTFe16].EnergyErg;

		TauLines[ipTFe26].PopLo = (float)(popFe2[1]);
		TauLines[ipTFe26].PopOpc = (float)(popFe2[1] - popFe2[5]);
		TauLines[ipTFe26].phots = (float)(popFe2[5]*data[1][5]);
		CoolHeavy.cFe26 = TauLines[ipTFe26].phots*TauLines[ipTFe26].EnergyErg;

		TauLines[ipTFe34].PopLo = (float)popFe2[2];
		TauLines[ipTFe34].PopOpc = (float)(popFe2[2] - popFe2[3]);
		TauLines[ipTFe34].phots = (float)(popFe2[3]*data[2][3]);
		CoolHeavy.cFe34 = TauLines[ipTFe34].phots*TauLines[ipTFe34].EnergyErg;

		TauLines[ipTFe35].PopLo = (float)popFe2[2];
		TauLines[ipTFe35].PopOpc = (float)(popFe2[2] - popFe2[4]);
		TauLines[ipTFe35].phots = (float)(popFe2[4]*data[2][4]);
		CoolHeavy.cFe35 = TauLines[ipTFe35].phots*TauLines[ipTFe35].EnergyErg;

		TauLines[ipTFe46].PopLo = (float)popFe2[3];
		TauLines[ipTFe46].PopOpc = (float)(popFe2[3] - popFe2[5]);
		TauLines[ipTFe46].phots = (float)(popFe2[5]*data[3][5]);
		CoolHeavy.cFe46 = TauLines[ipTFe46].phots*TauLines[ipTFe46].EnergyErg;

		TauLines[ipTFe56].PopLo = (float)popFe2[4];
		/* tFe56(ipLnPopOpc) = popFe2(5) - popFe2(6)
		 * >>chng 96 apr 23 this levels mases due to Lya pump */
		TauLines[ipTFe56].PopOpc = (float)(popFe2[4] - popFe2[5]*0.);
		TauLines[ipTFe56].phots = (float)(popFe2[5]*data[4][5]);
		CoolHeavy.cFe56 = TauLines[ipTFe56].phots*TauLines[ipTFe56].EnergyErg;

		/* Jack's funny FeII lines, data from 
		 * >>refer Johansson, S., Brage, T., Leckrone, D.S., Nave, G. &
		 * >>refer Wahlgren, G.M. 1995, ApJ 446, 361 */
		PutCS(10.,&TauLines[ipT191]);
		level2(&TauLines[ipT191]);

		/* do 16 level ground + first two terms */
		Fe2Lev16();
	}

	/* lump 3p and 3f together; cs=
	 * >>refer Garstang, R.H., Robb, W.D., Rountree, S.P. 1978, ApJ, 222, 384
	 * A from
	 * >>refer Garstang, R.H., 1957, Vistas in Astronomy, 1, 268
	 * FE III 5270, is 20.9% of total */
	CoolHeavy.c5270 = popexc(5.53,25.,30.,0.435,2.948e4,xIonFracs[3][25])*
	  3.78e-12;
	coladd("Fe 3",5270,CoolHeavy.c5270);

	/* FeIII 1122 entire multiplet - atomic data=A's Dima, CS = guess */
	PutCS(25.,&TauLines[ipT1122]);
	level2(&TauLines[ipT1122]);

	/* call 12 level atom */
	Fe4Lev12();

	/* FE V 3892 + 3839, data from Shields */
	CoolHeavy.c3892 = popexc(7.4,25.,5.,0.6,3.7e4,xIonFracs[5][25])*
	  5.11e-12;
	coladd("Fe 5",3892,CoolHeavy.c3892);

	/* FE VI 5177+5146+4972+4967
	 * data from 
	 * >>refer Garstang, R.H., Robb, W.D., Rountree, S.P. 1978, ApJ, 222, 384 */
	CoolHeavy.c5177 = popexc(1.9,28.,18.,0.52,2.78e4,xIonFracs[6][25])*
	  3.84e-12;
	coladd("Fe 6",5177,CoolHeavy.c5177);

	/* >>chng 97 jan 6 added fir lines, cs from 
	 * >>refer Keenan, F.P., & Norrington, P.H. 1987 A&A, 181, 370
	 * [FeVII] 9.51, 7.81 micron lines */
	if( phycon.te < 1.5e4 )
	{
		cs = MIN2(1.49,374.8/(TePowers.sqrte*TePowers.te10));
		PutCS(cs,&TauLines[ipFe0795]);
		cs = MIN2(1.57,188.5/(TePowers.sqrte*TePowers.te02));
		PutCS(cs,&TauLines[ipFe0778]);
		cs = MIN2(0.389,51.32/(TePowers.sqrte*TePowers.te03));
		PutCS(cs,&TauDummy);
	}
	else
	{
		cs = MIN2(1.17,7.929/(TePowers.te20/TePowers.te001));
		PutCS(cs,&TauLines[ipFe0795]);
		cs = MIN2(1.27,5.638/(TePowers.te10*TePowers.te05*TePowers.te005));
		PutCS(cs,&TauLines[ipFe0778]);
		cs = MIN2(0.314,1.491/(TePowers.te10*TePowers.te05*TePowers.te01*
		  TePowers.te001*TePowers.te001));
		PutCS(cs,&TauDummy);
	}
	level3(&TauLines[ipFe0795],&TauLines[ipFe0778],&TauDummy);

	/* FE VII 6087, A from 
	 * >>refer Nussbaumer, H., & Storey, P.J. 1982, A&A, 113, 21
	 * cs from 
	 * >>refer Keenan, F.P., & Norrington, P.H. 1987 A&A, 181, 370 */
	tused = MAX2(4.0,logte.alogte);
	tused = MIN2(5.1,tused);
	cs = 2.2684 - tused*(0.60893 - 0.052004*tused);
	/* >>chng 97 jan 27, should have been total collision strength of multiplet
	 * thanks for Tino Oliva for finding this problem */
	cs *= 2.96;
	CoolHeavy.c6087 = popexc(cs,21.,5.,0.94,2.32e4,xIonFracs[7][25])*
	  3.27e-12*0.577/0.94;
	CoolHeavy.Fe5722 = CoolHeavy.c6087*0.629*6087./5722.;
	tcool.dCooldT += CoolHeavy.c6087*2.37e4*tcool.tsq1;
	coladd("Fe 7",6087,CoolHeavy.c6087);
	coladd("Fe 7",5722,CoolHeavy.Fe5722);

	/* Fe IX 245,242
	 * all atomic data from 
	 * >>refer Flower, D.R. 1976, A&A, 56, 451 */
	PutCS(0.123,&TauLines[ipT245]);
	/* beseq(cs23,cs24,cs34,tarray,a41)
	 * C245 = BESEQ( .087 ,.038 ,.188 , t245 ,71. )  * 8.14E-11 */
	beseq(.087,.038,.188,&TauLines[ipT245],71.);

	CoolHeavy.c242 = PopLevls.PopLevels[3]*8.22e-11*71.;

	/* Fe X 6374, A from 
	 * >>refer Mason, H. 1975, MNRAS 170, 651
	 * cs 
	 * >>refer Mohan, M., Hibbert, A., & Kingston, A.E. 1994, ApJ, 434, 389
	 * not used newer cs from 
	 * >>refer Pelan, J., & Berrington, K.A. 1995, A&A Suppl, 110, 209
	 * does not agree with Mohan et al. by factor of 4
	 * 20x larger than Mason numbers used pre 1994 */
	cs = 13.3-2.*MIN2(7.0,logte.alogte);
	cs = MAX2(0.1,cs);
	/* 666 error! following to stop insane FeX strengths
	 * >>chng 96 july 11, removed 1/10 factor, so using real data, 90.01
	 * cs = cs * 0.1
	 * >>chng 96 june 3, transferred following
	 * >>chng 97 jan 31, give up on this mess, use cs of unity */
	cs = 1.;
	PutCS(cs,&TauLines[ipFe106375]);
	level2(&TauLines[ipFe106375]);
	/* c6374 = popexc(cs ,4.,2.,69.5,2.26e4,xIonFracs(26,10))*3.12e-12
	 * dCooldT = dCooldT + c6374 * 2.26e4 * tsq1
	 * call coladd( 'Fe10' , 6374 , c6374 )
	 *
	 * this is the E1 line that can pump [FeX] */
	cs = 0.85*sexp(0.045*1.259e6/phycon.te);
	cs = MAX2(0.05,cs);
	PutCS(cs,&TauLines[ipT352]);
	level2(&TauLines[ipT352]);

	/* Fe XI 7892, 6.08 mic, CS from 
	 * >>refer Mason, H. 1975, MNRAS 170, 651
	 * A from 
	 * >>refer Mendoza, C., & Zeippen, C.J. 1983, MNRAS, 202, 981 */
	cs = 0.27;
	/* >>chng 97 jan 31, set cs to unity, above ignored resonances */
	cs = 1.;
	PutCS(cs,&TauLines[ipTFe07]);
	PutCS(0.036,&TauLines[ipTFe61]);
	PutCS(0.11,&TauDummy);
	level3(&TauLines[ipTFe07],&TauLines[ipTFe61],&TauDummy);

	/* Kafatos, M., & Lynch, J.P. 1980, ApJS, 42, 611 */
	CoolHeavy.c1467 = pop3(9.,5.,1.,.24,.028,.342,100.,1012.,9.3,5.43e4,
	  6.19e4,&p2,xIonFracs[11][25],0.)*1012.*1.36e-11;
	coladd("Fe11",1467,CoolHeavy.c1467);
	CoolHeavy.c2649 = p2*91.0*7.512e-12;
	coladd("Fe11",2469,CoolHeavy.c2649);

	/* A's for 2-1 and 3-1 from 
	 * >>refer Tayal, S.S., & Henry, R.J.W. 1986, ApJ, 302, 200
	 * CS fro 
	 * >>refer Tayal, S.S., Henry, R.J.W., Pradhan, A.K. 1987, ApJ, 319, 951
	 * a(3-2) scaled from both ref */
	CoolHeavy.c2568 = pop3(4.,10.,6.,0.72,0.69,2.18,8.1e4,1.84e6,1.33e6,
	  6.37e4,4.91e4,&p2,xIonFracs[12][25],0.)*1.33e6*6.79e-12;
	coladd("Fe12",2568,CoolHeavy.c2568);
	CoolHeavy.c1242 = CoolHeavy.c2568*2.30*1.38;
	coladd("Fe12",1242,CoolHeavy.c1242);
	CoolHeavy.c2170 = p2*8.09e4*8.82e-12;
	coladd("Fe12",2170,CoolHeavy.c2170);

	/* [Fe XIII] 1.07, 1.08 microns, no collision data available */
	PutCS(1.,&TauLines[ipFe1310]);
	PutCS(1.,&TauLines[ipFe1311]);
	PutCS(1.,&TauDummy);
	level3(&TauLines[ipFe1310],&TauLines[ipFe1311],&TauDummy);

	/* Fe XIV 5303, cs from 
	 * >>refer Storey, P.J., Mason, H.E., Saraph, H.E., 1996, A&A, 309, 677 */
	fe14cs(logte.alogte,&cs);
	/* >>chng 97 jan 31, set cs to unity, above is VERY large, >10 */
	cs = 1.;
	CoolHeavy.c5303 = popexc(cs,2.,4.,60.3,2.71e4,xIonFracs[14][25])*
	  3.75e-12;
	tcool.dCooldT += CoolHeavy.c5303*2.71e4*tcool.tsq1;
	coladd("Fe14",5303,CoolHeavy.c5303);

	/* Fe 18 974.8A;, cs from 
	 * >>refer Saraph, H.E. & Tully, J.A. 1994, A&AS, 107, 29 */
	cs = MIN2(0.143,0.804/(TePowers.te10*TePowers.te05));
	PutCS(cs,&TauLines[ipFe18975]);
	level2(&TauLines[ipFe18975]);

	CoolHeavy.c7082 = pop3(5.,1.,3.,0.0132,0.0404,0.0137,0.505,1.46e4,
	  41.2,1.083e5,2.03e4,&p2,xIonFracs[19][25],0.)*41.2*
	  2.81e-12;
	CoolHeavy.c1118 = CoolHeavy.c7082*354.4*6.335;
	CoolHeavy.c1328 = p2*0.505*1.50e-11;
	coladd("Fe19",7082,CoolHeavy.c7082);
	coladd("Fe19",1118,CoolHeavy.c1118);
	coladd("Fe19",1328,CoolHeavy.c1328);

	CoolHeavy.c592 = popexc(0.0913,9.,5.,1.64e4,2.428e5,xIonFracs[19][25])*
	  3.36e-11;
	coladd("Fe19",592,CoolHeavy.c592);

	/* following all from same ref */
	PutCS(0.0151,&TauLines[ipTFe13]);
	PutCS(0.0424,&TauLines[ipTFe23]);
	PutCS(0.0147,&TauDummy);
	level3(&TauLines[ipTFe13],&TauLines[ipTFe23],&TauDummy);

	/* FE 22 845.6, C.S.=guess, A from
	 * >>refer Froese Fischer, C. 1983, J.Phys. B, 16, 157 */
	CoolHeavy.fs846 = popexc(0.2,2.,4.,1.5e4,1.701e5,xIonFracs[22][25])*
	  2.35e-11;
	coladd("Fe22",846,CoolHeavy.fs846);

	/* FE 23 262.6, 287A, 1909-LIKE, 
	 * cs from 
	 * >>refer Bhatia, A.K., & Mason, H.E. 1986, A&A, 155, 413 */
	CoolHeavy.c263 = popexc(0.04,1.,9.,1.6e7,5.484e5,xIonFracs[23][25])*
	  7.58e-11;
	coladd("Fe23",262,CoolHeavy.c263);

	/* FE 24, 255, 192 Li seq doublet
	 * >>refer Cochrane, D.M., & McWhirter, R.W.P. 1983, PhyS, 28, 25 */
	ligbar(26,&TauLines[ipT192],&TauLines[ipT11],&cs2s2p,&cs2s3p);

	PutCS(cs2s2p,&TauLines[ipT192]);
	level2(&TauLines[ipT192]);

	/* funny factor (should have been 0.5) due to energy change */
	PutCS(cs2s2p*0.376,&TauLines[ipT255]);
	level2(&TauLines[ipT255]);

	PutCS(cs2s3p,&TauLines[ipT11]);
	level2(&TauLines[ipT11]);

	/* 666 error!
	 * following not in cooling function */
	TauLines[ipT353].PopOpc = xIonFracs[11][25];
	TauLines[ipT353].PopLo = xIonFracs[11][25];
	TauLines[ipT353].PopHi = 0.;
	TauLines[ipT347].PopOpc = xIonFracs[14][25];
	TauLines[ipT347].PopLo = xIonFracs[14][25];
	TauLines[ipT347].PopHi = 0.;


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