/*PrtFinal create PrtFinal pages of printout, emission line intensities, etc */
#include <string.h>
#include <float.h>
#include <limits.h>
#include "cddefines.h"
#include "cddrive.h"
#include "physconst.h"
#include "linesave.h"
#include "nhe1lvl.h"
#include "taulines.h"
#include "grainvar.h"
#include "hydrogenic.h"
#include "colmas.h"
#include "jmin.h"
#include "hsrate.h"
#include "densty.h"
#include "timesc.h"
#include "opac.h"
#include "tehigh.h"
#include "pressure.h"
#include "badstp.h"
#include "converge.h"
#include "itagan.h"
#include "warngs.h"
#include "uspher.h"
#include "qs.h"
#include "sphere.h"
#include "called.h"
#include "pshort.h"
#include "fourpi.h"
#include "faint.h"
#include "printit.h"
#include "itercnt.h"
#include "zonecnt.h"
#include "norm.h"
#include "tff.h"
#include "abrems.h"
#include "date.h"
#include "he.h"
#include "u.h"
#include "phycon.h"
#include "heat.h"
#include "colden.h"
#include "tcool.h"
#include "input.h"
#include "rfield.h"
#include "nh.h"
#include "radius.h"
#include "xry.h"
#include "peimbt.h"
#include "o3data.h"
#include "filfac.h"
#include "prnline.h"
#include "betaver.h"
#include "htopoff.h"
#include "hphoto.h"
#include "ipoint.h"
#include "cap4.h"
#include "totlin.h"
#include "prtalltau.h"
#include "aver.h"
#include "molcol.h"
#include "sexp.h"
#include "gett2.h"
#include "gett2o3.h"
#include "pmeani.h"
#include "spect.h"
#include "pcontn.h"
#include "wind.h"
#include "printe82.h"
#include "prtfinal.h"

/* return is TRUE if calculation ok, false otherwise */
void PrtFinal(void)
{
	short int *lgPrt;
	long int *linsav;
	double *sclsav , *scaled , *asumli ;
	char **chSLab;
	char *chSTyp;
	char chCKey[5], 
	  chGeo[7], 
	  chLab[NDUST][11], 
	  chPlaw[21] ;

	long int 
	  i, 
	  iindex[32], 
	  ip2500, 
	  ip2kev, 
	  iprnt, 
	  j, 
	  nd, 
	  nline, 
	  nneg;
	double o4363, 
	  bacobs, 
	  absint, 
	  bacthn, 
	  hbcab, 
	  hbeta, 
	  o5007;

	double a, 
	  ajmass, 
	  ajmin, 
	  alfox, 
	  bot, 
	  bottom, 
	  bremtk, 
	  coleff, 
	  /* nb 8 is used for following preset in loop */
	  d[8], 
	  dmean, 
	  ferode, 
	  he4686, 
	  he5876, 
	  heabun, 
	  hescal, 
	  pion, 
	  plow, 
	  powerl, 
	  qion, 
	  qlow, 
	  rate, 
	  ratio, 
	  reclin, 
	  rjeans, 
	  snorm,
	  tmean, 
	  top, 
	  THI,/* HI 21 cm spin temperature */
	  uhel, 
	  uhl, 
	  usp, 
	  wmean, 
	  znit;
	 
	  double bac, 
	  tel, 
	  x;

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

	/* return if not talking */
	if( !called.lgTalk )
	{ 
#		ifdef DEBUG_FUN
		fputs( " <->PrtFinal()\n", debug_fp );
#		endif
		return;
	}

	/* print out header, or parts that were saved */
	
	/* this would be a major logical error */
	assert( LineSave.nsum > 1 );

	/* print name and version number */
	fprintf( ioQQQ, "\f\n                       *********************************> Cloudy ");
	fprintf( ioQQQ, "%7.7s <********************************\n", date.chVersion );

	for( i=0; i <= input.nSave; i++ )
	{
		/* print command line, unless it is a continue line */

		/* copy start of command to key, making it into capitals */
		cap4(chCKey,input.chCardSav[i]);

		/* print if not continue */
		if( strcmp(chCKey,"CONT") != 0 )
		{
			/* print leading space and first * */
			fprintf( ioQQQ, "                       * ");

			j = 0;
			/* print actual command */
			while( input.chCardSav[i][j] !='\0' )
			{
				fprintf( ioQQQ, "%c", input.chCardSav[i][j] );
				++j;
			}

			/* flush out rest of space then * */
			while( j<80 )
			{
				fprintf( ioQQQ, "%c", ' ' );
				++j;
			}
			/* end the line */
			fprintf( ioQQQ, "*\n" );
		}
	}

	if( u.uh > 0. )
	{
		a = log10(u.uh);
	}
	else
	{
		a = -37.;
	}

	fprintf( ioQQQ, 
		"                       *********************************> Log(U):%6.2f <*********************************\n", 
	  a );

	if( BetaVer.nBetaVer > 0 )
	{
		fprintf( ioQQQ, 
			"\n                        This is a beta test version of the code, and is intended for testing only.\n\n" );
	}

	if( warngs.lgWarngs )
	{
		fprintf( ioQQQ, "  \n" );
		fprintf( ioQQQ, "                       >>>>>>>>>> Warning!\n" );
		fprintf( ioQQQ, "                       >>>>>>>>>> Warning!\n" );
		fprintf( ioQQQ, "                       >>>>>>>>>> Warning!  Warnings exist, this calculation has serious problems.\n" );
		fprintf( ioQQQ, "                       >>>>>>>>>> Warning!\n" );
		fprintf( ioQQQ, "                       >>>>>>>>>> Warning!\n" );
		fprintf( ioQQQ, "  \n" );
	}
	else if( warngs.lgCautns )
	{
		fprintf( ioQQQ, 
			"                       >>>>>>>>>> Cautions are present.\n" );
	}

	if( badstp.lgBadStop )
	{
		fprintf( ioQQQ, 
			"                       >>>>>>>>>> The calculation stopped for unintended reasons.\n" );
	}

	if( itagan.lgIterAgain )
	{
		fprintf( ioQQQ, 
			"                       >>>>>>>>>> Another iteration is needed.\n" );
	}

	LineSave.arealg = fourpi.pirsq;

	/* open or closed geometry?? */
	if( sphere.lgSphere )
	{
		strcpy( chGeo, "Closed" );
	}
	else
	{
		strcpy( chGeo, "  Open" );
	}

	/* now give description of pressure law */
	if( strcmp(pressure.chCPres,"CPRE") == 0 )
	{
		strcpy( chPlaw, " Constant  Pressure " );
	}

	else if( strcmp(pressure.chCPres,"CDEN") == 0 )
	{
		strcpy( chPlaw, "  Constant  Density " );
	}

	else if( (strcmp(pressure.chCPres,"POWD") == 0 || strcmp(pressure.chCPres
	  ,"POWR") == 0) || strcmp(pressure.chCPres,"POWC") == 0 )
	{
		strcpy( chPlaw, "  Power Law Density " );
	}

	else if( strcmp(pressure.chCPres,"SINE") == 0 )
	{
		strcpy( chPlaw, " Rapid Fluctuations " );
	}

	else if( strncmp(pressure.chCPres , "DLW" , 3) == 0 )
	{
		strcpy( chPlaw, " Special Density Lw " );
	}

	else if( strcmp(pressure.chCPres,"TIME") == 0 )
	{
		strcpy( chPlaw, "   Time  Dependent  " );
	}

	else if( strcmp(pressure.chCPres,"HYDR") == 0 )
	{
		strcpy( chPlaw, " Hydrostatic Equlib " );
	}

	else if( strcmp(pressure.chCPres,"WIND") == 0 )
	{
		strcpy( chPlaw, "  Radia Driven Wind " );
	}

	else if( strcmp(pressure.chCPres,"GLOB") == 0 )
	{
		strcpy( chPlaw, " Globule            " );
	}

	else
	{
		strcpy( chPlaw, " UNRECOGNIZED CPRES " );
	}

	fprintf( ioQQQ, "\n                     Emission Line Spectrum. %20.20sModel.  %6.6s geometry.  Iteration%2ld of%2ld.\n", 
	  chPlaw, chGeo, IterCnt.iter, IterCnt.itermx + 1 );

	/* emission lines as logs of total luminosity
	 * either per unit inner area (lgFourPi=.F.), or FOUR PI R SQRD (=.T.) */
	if( fourpi.lgFourPi )
	{
		/* four pi r^2, prog works in sqr cm, multiply by this */
		if( radius.lgCylnOn )
		{
			fprintf( ioQQQ, 
				"                                         Luminosity (erg/s) emitted by partial  cylinder.\n" );
		}

		else if( sphere.covgeo == 1. )
		{
			fprintf( ioQQQ, 
				"                                     Luminosity (erg/s) emitted by shell with full coverage.\n" );
		}

		else
		{
			fprintf( ioQQQ, 
				"                               Luminosity (erg/s) emitted by shell with a covering factor of%6.1f%%.\n", 
			  sphere.covgeo*100. );
		}
	}

	else
	{
		fprintf( ioQQQ, "                                                    Intensity (erg/s/cm^2)\n" );
	}

	/* throw a blank line */
	fprintf( ioQQQ, " %c\n", ' ' );

	/* this is the intensity of the line spectrum will be normalized to */
	snorm = LineSv[norm.ipNormWavL].sumlin;

	/* check that this line has positive intensity */
	if( ((snorm == 0.) || (norm.ipNormWavL < 0)) || (norm.ipNormWavL > 
	  LineSave.nsum) )
	{
		fprintf( ioQQQ, " >>Normalization line has zero intensity, its intensity was set to 1.  \n >>Please consider using another normalization line (this is set with the NORMALIZE command).\n" );
		fprintf( ioQQQ, " >>The relative intensities will be meaningless, and many lines may appear too faint.\n" );
		snorm = 1.;
	}

	/******************************************************************
	 * kill Hummer & Storey case b predictions if outside their table *
	 ******************************************************************/

	if( !CaseBHS.lgHCaseBOK )
	{
		/* table exceeded - kill all H case b predictions*/
		for( i=0; i < LineSave.nsum; i++ )
		{
			if( strcmp( LineSv[i].chALab,"Ca B" )==0 )
			{
				if( LineSv[i].lin==40   || LineSv[i].lin==6563 ||
					 LineSv[i].lin==4861 || LineSv[i].lin==18   ||
					 LineSv[i].lin==13   || LineSv[i].lin==12   ||
					 LineSv[i].lin==4340 || LineSv[i].lin==10   ||
					 LineSv[i].lin==21   || LineSv[i].lin==74   ||
					 LineSv[i].lin==46   || LineSv[i].lin==26
					 )
				{
					LineSv[i].sumlin = 0.;
				}
				else if( LineSv[i].lin==37 )
				{
					LineSv[i].sumlin = 0.;
					/* this is the last H line - bail */
					break;
				}
			}
		}
	}

	if( !CaseBHS.lgHeCaseBOK )
	{
		/* table exceeded - kill all He case b predictions*/
		for( i=0; i < LineSave.nsum; i++ )
		{
			if( strcmp( LineSv[i].chALab,"Ca B" )==0 )
			{
				if( LineSv[i].lin==1640 || LineSv[i].lin==1215 ||
					 LineSv[i].lin==1085 || LineSv[i].lin==4686 ||
					 LineSv[i].lin==3203 || LineSv[i].lin==2733 ||
					 LineSv[i].lin==10   || LineSv[i].lin==6560 ||
					 LineSv[i].lin==5412 || LineSv[i].lin==18   ||
					 LineSv[i].lin==11   
					 )
				{
					LineSv[i].sumlin = 0.;
				}
				/* this is the last HeII line, bail */
				else if( LineSv[i].lin==9345 )
				{
					LineSv[i].sumlin = 0.;
					/* this is the last He line - bail */
					break;
				}
			}
		}
	}

	/**********************************************************
	 *analyse line arrays for abundances, temp, etc           *
	 **********************************************************/

	/* find apparent helium abundance, will not find these if helium is turned off */

	if( !cdLine("TOTL",4861,&hbeta,&absint) )
	{
		fprintf( ioQQQ, " PrtFinal could not find TOTL 4861 with cdLine.\n" );
		puts( "[Stop in PrtFinal]" );
		exit(1);
	}

	/* get HeI 5876 */
	if( !cdLine("TOTL",5876,&he5876,&absint) )
	{
		he5876 = 0.;
	}

	/* get HeII 4686 */
	if( !cdLine("He 2",4685,&he4686,&absint) )
	{
		he4686 = 0.;
	}

	if( hbeta > 0. )
	{
		heabun = (he4686*0.078 + he5876*0.739)/hbeta;
	}
	else
	{
		heabun = 0.;
	}

	hescal = heabun/(he.ahe/phycon.hden);

	/* get temperature from [OIII] spectrum, o may be turned off,
	 * but lines still dumped into stack */
	if( !cdLine("O  3",5007,&o5007,&absint) )
	{
		fprintf( ioQQQ, " PrtFinal could not find O  3 5007 with cdLine.\n" );
		puts( "[Stop in PrtFinal]" );
		exit(1);
	}

	if( !cdLine("TOTL",4363,&o4363,&absint) )
	{
		fprintf( ioQQQ, " PrtFinal could not find TOTL 4363 with cdLine.\n" );
		puts( "[Stop in PrtFinal]" );
		exit(1);
	}

	/* first find low density limit OIII zone temp */
	if( (o4363 != 0.) && (o5007 != 0.) )
	{
		/* following will assume coll excit only, so correct
		 * for 4363's that cascade as 5007 */
		bot = o5007 - o4363/o3data.o3enro;

		if( bot > 0. )
		{
			ratio = o4363/bot;
		}
		else
		{
			ratio = 0.178;
		}

		if( ratio > 0.177 )
		{
			/* ratio was above infinite temperature limit */
			peimbt.toiiir = 1e10;
		}
		else
		{
			/* data for following set in oiii cooling routine
			 * ratio of collision strengths, fac of 4/3 makes 5007+4959
			 * >>chng 96 sep 7, reset cs to values at 10^4K,
			 * had been last temp in model */
			o3data.o3cs12 = 2.33239f;
			o3data.o3cs13 = 0.292628f;
			ratio = ratio/1.3333/(o3data.o3cs13/o3data.o3cs12);
			/* ratio of energies and branching ratio for 4363 */
			ratio = ratio/o3data.o3enro/o3data.o3br32;
			peimbt.toiiir = (float)(-o3data.o3ex23/log(ratio));
		}
	}

	else
	{
		peimbt.toiiir = 0.;
	}

	/* find temperature from Balmer continuum */
	if( !cdLine("Bac ",3646,&bacobs,&absint) )
	{
		fprintf( ioQQQ, " PrtFinal could not find Bac 3546 with cdLine.\n" );
		puts( "[Stop in PrtFinal]" );
		exit(1);
	}

	/* we pulled hbeta out of stack with cdLine above */
	if( hbeta > 0. )
	{
		bac = bacobs/hbeta;
	}
	else
	{
		bac = 0.;
	}

	if( bac > 0. )
	{
		/*----------------------------------------------------------*
		 ***** TableCurve c:\tcwin2\balcon.for Sep 6, 1994 11:13:38 AM
		 ***** log bal/Hbet
		 ***** X= log temp
		 ***** Y= 
		 ***** Eqn# 6503  y=a+b/x+c/x^2+d/x^3+e/x^4+f/x^5
		 ***** r2=0.9999987190883581
		 ***** r2adj=0.9999910336185065
		 ***** StdErr=0.001705886369042427
		 ***** Fval=312277.1895753243
		 ***** a= -4.82796940090671
		 ***** b= 33.08493347410885
		 ***** c= -56.08886262205931
		 ***** d= 52.44759454803217
		 ***** e= -25.07958990094209
		 ***** f= 4.815046760060006
		 *----------------------------------------------------------*
		 * bac is double precision!!!! */
		x = 1.e0/log10(bac);
		tel = -4.827969400906710 + x*(33.08493347410885 + x*(-56.08886262205931 + 
		  x*(52.44759454803217 + x*(-25.07958990094209 + x*(4.815046760060006)))));

		if( tel > 1. && tel < 5. )
		{
			peimbt.tbac = (float)pow(10.,tel);
		}
		else
		{
			peimbt.tbac = 0.;
		}
	}

	else
	{
		peimbt.tbac = 0.;
	}

	/* find temperature from optically thin Balmer continuum and case B H-beta */
	if( !cdLine("thin",3646,&bacthn,&absint) )
	{
		fprintf( ioQQQ, " PrtFinal could not find thin 3646 with cdLine.\n" );
		puts( "[Stop in PrtFinal]" );
		exit(1);
	}

	if( !cdLine("Ca B",4861,&hbcab,&absint) )
	{
		fprintf( ioQQQ, " PrtFinal could not find Ca B 4861 with cdLine.\n" );
		puts( "[Stop in PrtFinal]" );
		exit(1);
	}

	if( hbcab > 0. )
	{
		bacthn /= hbcab;
	}
	else
	{
		bacthn = 0.;
	}

	if( bacthn > 0. )
	{
		/*----------------------------------------------------------*
		 ***** TableCurve c:\tcwin2\balcon.for Sep 6, 1994 11:13:38 AM
		 ***** log bal/Hbet
		 ***** X= log temp
		 ***** Y= 
		 ***** Eqn# 6503  y=a+b/x+c/x^2+d/x^3+e/x^4+f/x^5
		 ***** r2=0.9999987190883581
		 ***** r2adj=0.9999910336185065
		 ***** StdErr=0.001705886369042427
		 ***** Fval=312277.1895753243
		 ***** a= -4.82796940090671
		 ***** b= 33.08493347410885
		 ***** c= -56.08886262205931
		 ***** d= 52.44759454803217
		 ***** e= -25.07958990094209
		 ***** f= 4.815046760060006
		 *----------------------------------------------------------*
		 * bac is double precision!!!! */
		x = 1.e0/log10(bacthn);
		tel = -4.827969400906710 + x*(33.08493347410885 + x*(-56.08886262205931 + 
		  x*(52.44759454803217 + x*(-25.07958990094209 + x*(4.815046760060006)))));

		if( tel > 1. && tel < 5. )
		{
			peimbt.tbcthn = (float)pow(10.,tel);
		}
		else
		{
			peimbt.tbcthn = 0.;
		}
	}
	else
	{
		peimbt.tbcthn = 0.;
	}

	/* we now have temps from OIII ratio and BAC ratio, now
	 * do Peimbert analysis, getting To and t^2 */

	peimbt.tohyox = (float)((8.5*peimbt.toiiir - 7.5*peimbt.tbcthn - 228200. + 
	  sqrt(POW2(8.5*peimbt.toiiir-7.5*peimbt.tbcthn-228200.)+9.128e5*
	  peimbt.tbcthn))/2.);

	if( peimbt.tohyox > 0. )
	{
		peimbt.t2hyox = (float)((peimbt.tohyox - peimbt.tbcthn)/(1.7*peimbt.tohyox));
	}
	else
	{
		peimbt.t2hyox = 0.;
	}

	/*----------------------------------------------------------
	 *
	 * first set scaled lines */

	/* get space for scaled */
	scaled = (double *)malloc( sizeof(double)*(unsigned long)LineSave.nsum );
	if( scaled == NULL )
	{ 
		printf( " not enough memory to allocate scaled in PrtFinal\n" );
		puts( "[Stop in PrtFinal]" );
		exit(1);
	}

	/* get space for asumli */
	asumli = (double *)malloc( sizeof(double)*(unsigned long)LineSave.nsum );
	if( asumli == NULL )
	{ 
		printf( " not enough memory to allocate asumli in PrtFinal\n" );
		puts( "[Stop in PrtFinal]" );
		exit(1);
	}

	/* this is option to not print certain contributions */
	/* gjf 98 june 10, get space for array lgPrt */
	lgPrt = (short int *)malloc( sizeof(short int)*(unsigned long)LineSave.nsum );
	if( lgPrt == NULL )
	{ 
		printf( " not enough memory to allocate lgPrt in PrtFinal\n" );
		puts( "[Stop in PrtFinal]" );
		exit(1);
	}

	for( i=0; i < LineSave.nsum; i++ )
	{
		scaled[i] = 
			LineSv[i].sumlin/snorm*norm.ScaleNormLine;

		if( LineSv[i].sumlin > 0. )
		{
			asumli[i] = log10(LineSv[i].sumlin) + fourpi.pirsq;
		}
		else
		{
			asumli[i] = -38.;
		}
	}

	for( i=0; i < LineSave.nsum; i++ )
	{
		if( strcmp(LineSv[i].chALab,"Coll") == 0 && !PrnLine.lgPrnColl )
		{
			lgPrt[i] = FALSE;
		}
		else if( strcmp(LineSv[i].chALab,"Pump") == 0 && !PrnLine.lgPrnPump )
		{
			lgPrt[i] = FALSE;
		}
		else if( strncmp(LineSv[i].chALab,"Inw",3) == 0 && !PrnLine.lgPrnInwd )
		{
			lgPrt[i] = FALSE;
		}
		else if( strcmp(LineSv[i].chALab,"Heat") == 0 && !PrnLine.lgPrnHeat )
		{
			lgPrt[i] = FALSE;
		}
		else if( strcmp(LineSv[i].chALab,"nFnu") == 0 && !PrnLine.lgPrnDiff )
		{
			lgPrt[i] = FALSE;
		}
		else
		{
			lgPrt[i] = TRUE;
		}
	}

	/* do not print relatively faint lines unless requested */
	nneg = 1;

	/* get space for linsav */
	linsav = (long int *)malloc( sizeof(long int)*(unsigned long)LineSave.nsum );
	if( linsav == NULL )
	{ 
		printf( " not enough memory to allocate linsav in PrtFinal\n" );
		puts( "[Stop in PrtFinal]" );
		exit(1);
	}

	/* get space for sclsav */
	sclsav = (double *)malloc( sizeof(double)*(unsigned long)LineSave.nsum );
	if( sclsav == NULL )
	{ 
		printf( " not enough memory to allocate sclsav in PrtFinal\n" );
		puts( "[Stop in PrtFinal]" );
		exit(1);
	}

	/* get space for chSTyp */
	chSTyp = (char *)malloc( sizeof(char)*(unsigned long)LineSave.nsum );
	if( chSTyp == NULL )
	{ 
		printf( " not enough memory to allocate chSTyp in PrtFinal\n" );
		puts( "[Stop in PrtFinal]" );
		exit(1);
	}

	/* get space for chSLab,
	 * first define array of pointers*/
	/* char chSLab[NLINES][5];*/
	chSLab = ((char**)malloc((unsigned long)LineSave.nsum*sizeof(char*)));
	if( chSLab == NULL )
	{ 
		printf( " not enough memory to allocate 1st chSTyp in PrtFinal\n" );
		puts( "[Stop in PrtFinal]" );
		exit(1);
	}

	/* now allocate all the labels for each of the above lines */
	for( i=0; i<LineSave.nsum; ++i)
	{
		chSLab[i] = (char*)malloc(5*sizeof(char));
		if( chSLab[i] == NULL )
		{ 
			printf( " not enough memory to allocate 2st chSTyp in PrtFinal\n" );
			puts( "[Stop in PrtFinal]" );
			exit(1);
		}
	}


	/* set iindex to bomb to make sure set in following */
	for(i=0; i< 32; i++ )
	{
		iindex[i] = LONG_MAX;
	}

	for(i=0;i<8;++i)
	{
		d[i] = -DBL_MAX;
	}

	if( (GrainVar.lgDustOn && (!sphere.lgSphere)) && (sphere.covrt == 0.) )
	{
		long int nPrint=0;

		/* print special lines with continuum transfer in */
		fprintf( ioQQQ, 
			"                                                   Emergent line intensities\n" );
		bottom = MAX2(1e-30,LineDSv[0].smdlin);
		/* problem with too many dust lines here */

		fprintf( ioQQQ, "      " );
		for( i=0; i < LineSave.ndsum; ++i )
		{
			/*limit = MIN2(nlow+4,LineSave.ndsum);*/

			/* this is the actual print statement for the massive set of emission lines 
			 * when grains are present*/
			if( LineDSv[i].smdlin/bottom >  faint.TooFaint)
			{
				fprintf( ioQQQ, "%4.4s%5ld%7.3f%9.4f", 
				   LineDSv[i].chSMDLab ,
				   LineDSv[i].lind, 
				   log10(MAX2(LineDSv[i].smdlin,1e-35)) + fourpi.pirsq, 
				   LineDSv[i].smdlin/bottom );
				fprintf( ioQQQ, "      " );
				++nPrint;
			}
			if( nPrint==4 )
			{
				fprintf( ioQQQ, "\n" );
				nPrint = 0;
				fprintf( ioQQQ, "      " );
			}
		}
		if( nPrint !=0 )
		{
			fprintf( ioQQQ, "\n" );
		}
		fprintf( ioQQQ, 
			"\n                                                   Intrinsic line intensities\n" );
	}

	/* option to only print brighter lines */
	if( faint.lgFaintOn )
	{
		iprnt = 0;
		for( i=0; i < LineSave.nsum; i++ )
		{
			if( scaled[i] >= faint.TooFaint && lgPrt[i] )
			{
				asumli[iprnt] = log10(LineSv[i].sumlin) + fourpi.pirsq;
				sclsav[iprnt] = scaled[i];
				chSTyp[iprnt] = LineSv[i].chSumTyp;
				strcpy( chSLab[iprnt], LineSv[i].chALab );
				linsav[iprnt] = LineSv[i].lin;
				iprnt += 1;
			}
			else if( -scaled[i] > faint.TooFaint && lgPrt[i] )
			{
				/* negative intensities occur if line absorbs continuum */
				iindex[nneg-1] = i+1;
				iindex[nneg] = LineSv[i].lin;
				nneg = MIN2(32,nneg+2);
			}
		}
	}

	else
	{
		/* print everything */
		iprnt = LineSave.nsum;
		for( i=0; i < LineSave.nsum; i++ )
		{
			sclsav[i] = scaled[i];
			chSTyp[i] = LineSv[i].chSumTyp;
			strcpy( chSLab[i], LineSv[i].chALab );
			linsav[i] = LineSv[i].lin;
			if( scaled[i] < 0. )
			{
				iindex[nneg-1] = i+1;
				iindex[nneg] = LineSv[i].lin;
				nneg = MIN2(32,nneg+2);
			}
		}
	}

	/*----------------------------------------------------------
	 *
	 * print out all lines which made it through the faint filter */
	/* nline will be the number of horizontal lines - 
	 * the 4 represents the 4 columns */
	nline = (iprnt + 3)/4;
	for( i=0; i < nline; i++ )
	{
		fprintf( ioQQQ, "      " );

		/* this skips over nline per step, to go to the next colum in the output */
		for( j=i; j<iprnt; j = j + nline)
		{ 
			/* this is the actual print statement for the emission line spectrum*/

			/* the label for the line */
			fprintf( ioQQQ, "%4.4s",chSLab[j] ); 

			/* the integer wavelength */
			fprintf( ioQQQ, "%5ld",	linsav[j]);

			/* log of the intensity of the lines */
			fprintf( ioQQQ, "%7.3f", asumli[j] );

			/* print scaled intensity with various formats,
			 * depending on order of magnitude.  value is
			 * always the same but the format changes. */
			if( sclsav[j] < 9999.5 )
			{
				fprintf( ioQQQ, "%9.4f",sclsav[j] );
			}
			else if( sclsav[j] < 99999.5 )
			{
				fprintf( ioQQQ, "%9.3f",sclsav[j] );
			}
			else if( sclsav[j] < 999999.5 )
			{
				fprintf( ioQQQ, "%9.2f",sclsav[j] );
			}
			else if( sclsav[j] < 9999999.5 )
			{
				fprintf( ioQQQ, "%9.1f",sclsav[j] );
			}
			else
			{
				fprintf( ioQQQ, "*********" );
			}

			/* now end the block with some spaces to set next one off */
			fprintf( ioQQQ, "      " );
		}
		/* now end the lines */
		fprintf( ioQQQ, "      \n" );
	}

	if( nneg > 1 )
	{
		fprintf( ioQQQ, " Lines with negative intensities -  Linear itensities relative to normalize line.\n" );
		fprintf( ioQQQ, "          " );

		for( i=0; i < (nneg - 1); i += 2 )
		{
			fprintf( ioQQQ, "%5ld%5.5s%5ld%9.1e", iindex[i], 
			  LineSv[iindex[i]-1].chALab, iindex[i+1], 
			  scaled[iindex[i]-1] );
		}
		fprintf( ioQQQ, "\n" );
	}

	/* add lines to truncated list, using known line ratios, with printou */
	LineSave.npxdd = iprnt;

	/* with area since sumlin does */
	snorm = LineSv[norm.ipNormWavL].sumlin;

	/* reorder lines according to wavelength for comparison with spectrum
	 * including writing out the results */
	if( printit.psort == 1 )
		spect();

	/* now find which were the very strongest, more that 5% of cooling */
	j = 0;
	for( i=0; i < LineSave.nsum; i++ )
	{
		a = (double)LineSv[i].sumlin/(double)tcool.totcol;
		if( (a >= 0.05) && LineSv[i].chSumTyp == 'c' )
		{
			iindex[j] = i;
			d[j] = a;
			j = MIN2(j+1,7);
		}
	}

	fprintf( ioQQQ, "\n\n\n %s\n", input.chTitle );
	if( j != 0 )
	{
		fprintf( ioQQQ, " Cooling: " );
		for( i=0; i < j; i++ )
		{
			fprintf( ioQQQ, " %4.4s%5ld:%5.3f", 
				LineSv[iindex[i]].chALab, LineSv[iindex[i]].lin, d[i] );
		}
		fprintf( ioQQQ, "  \n" );
	}

	/* now find strongest heating, more that 5% of total */
	j = 0;
	for( i=0; i < LineSave.nsum; i++ )
	{
		a = (double)LineSv[i].sumlin/(double)heat.power;
		if( (a >= 0.05) && LineSv[i].chSumTyp == 'h' )
		{
			iindex[j] = i;
			d[j] = a;
			j = MIN2(j+1,7);
		}
	}

	if( j != 0 )
	{
		fprintf( ioQQQ, " Heating: " );
		for( i=0; i < j; i++ )
		{
			fprintf( ioQQQ, " %4.4s%5ld:%5.3f", 
				LineSv[iindex[i]].chALab, LineSv[iindex[i]].lin, d[i] );
		}
		fprintf( ioQQQ, "  \n" );
	}

	/* print out ionization parameters and radiation making it through */
	qlow = 0.;
	plow = 0.;
	for( i=0; i < (nh.ipHn[0][IP1S] - 1); i++ )
	{
		/* N.B. in following en1ryd prevents overflow */
		plow += (rfield.flux[i] + rfield.outcon[i] + rfield.outlin[i])*
		  EN1RYD*rfield.anu[i];
		qlow += rfield.flux[i] + rfield.outcon[i] + rfield.outlin[i];
	}

	qlow *= radius.r1r0sq;
	plow *= radius.r1r0sq;
	if( plow > 0. )
	{
		qlow = log10(qlow) + fourpi.pirsq;
		plow = log10(plow) + fourpi.pirsq;
	}
	else
	{
		qlow = 0.;
		plow = 0.;
	}

	qion = 0.;
	pion = 0.;
	for( i=nh.ipHn[0][IP1S]-1; i < rfield.nflux; i++ )
	{
		/* N.B. in following en1ryd prevents overflow */
		pion += (rfield.flux[i] + rfield.outcon[i] + rfield.outlin[i])*
		  EN1RYD*rfield.anu[i];
		qion += rfield.flux[i] + rfield.outcon[i] + rfield.outlin[i];
	}

	qion *= radius.r1r0sq;
	pion *= radius.r1r0sq;

	if( pion > 0. )
	{
		qion = log10(qion) + fourpi.pirsq;
		pion = log10(pion) + fourpi.pirsq;
	}
	else
	{
		qion = 0.;
		pion = 0.;
	}

	/* derive ionization parameter for spherical geometry */
	if( qs.qhtot > 0. )
	{
		if( uspher.lgUSphON )
		{
			/* RSTROM is stromgren radius */
			usp = qs.qhtot/POW2(uspher.rstrom/radius.rinner)/
			  2.998e10/phycon.hden;
			usp = log10(usp);
		}
		else
		{
			/* no stromgren radius found, use outer radius */
			usp = qs.qhtot/radius.r1r0sq/2.998e10/phycon.hden;
			usp = log10(usp);
		}
	}

	else
	{
		usp = -37.;
	}

	if( u.uh > 0. )
	{
		uhl = log10(u.uh);
	}
	else
	{
		uhl = -37.;
	}

	if( u.uheii > 0. )
	{
		uhel = log10(u.uheii);
	}
	else
	{
		uhel = -37.;
	}

	fprintf( ioQQQ, "\n IONIZE PARMET:  U(1-)%8.4f  U(4-):%8.4f  U(sp):%6.2f  Q(ion):   %6.3f  L(ion):%8.3f Q(low):%6.2f    P(low)%9.3f\n", 
	  uhl, uhel, usp, qion, pion, qlow, plow );

	a = fabs((heat.power-tcool.totcol)*100.)/heat.power;
	/* output power and total cooling; can be neg for shocks, collisional ioniz */
	if( heat.power > 0. )
	{
		powerl = log10(heat.power) + fourpi.pirsq;
	}
	else
	{
		powerl = 0.;
	}

	if( tcool.totcol > 0. )
	{
		tcool.totcol = log10(tcool.totcol) + fourpi.pirsq;
	}
	else
	{
		tcool.totcol = 0.;
	}

	if( hydro.FreeFreeTotHeat > 0. )
	{
		hydro.FreeFreeTotHeat = log10(hydro.FreeFreeTotHeat) + fourpi.pirsq;
	}
	else
	{
		hydro.FreeFreeTotHeat = 0.;
	}

	/* following is recombination line intensity */
	reclin = totlin('r');
	if( reclin > 0. )
	{
		reclin = log10(reclin) + fourpi.pirsq;
	}
	else
	{
		reclin = 0.;
	}

	fprintf( ioQQQ, 
		" ENERGY BUDGET:  Heat:%8.3f  Coolg:%8.3f  Error:%5.1f%%  Rec Lin:%8.3f  F-F  H%7.3f    RadBetaMax:", 
	  powerl, tcool.totcol, a, reclin, hydro.FreeFreeTotHeat );
	PrintE82( ioQQQ, pressure.RadBetaMax );
	fprintf( ioQQQ, "\n");

	/* effective x-ray column density, from 0.5keV attenuation, no scat
	 * IPXRY is pointer to 73.5 Ryd */
	coleff = opac.tauabs[0][xry.ipxry-1] - xry.tauxry;
	coleff /= 2.14e-22;

	fprintf( ioQQQ, "\n Column density  H12:");
	PrintE93(ioQQQ, coldenCom.colden[IPCOLUMNDENSITY-1]);
	fprintf( ioQQQ, "  H II:");
	PrintE93(ioQQQ, coldenCom.colden[IPCHII-1]);
	fprintf( ioQQQ, "  HI:");
	PrintE93(ioQQQ, coldenCom.colden[IPCHI-1]);
	fprintf( ioQQQ, "  H-: ");
	PrintE93(ioQQQ, coldenCom.colden[IPCHMIN-1]);
	fprintf( ioQQQ, "  H2: ");
	PrintE93(ioQQQ, coldenCom.colden[IPCOLH2-1]);
	fprintf( ioQQQ, "  H2+:");
	PrintE93(ioQQQ, coldenCom.colden[IPCH2PLS-1]);
	fprintf( ioQQQ, " HeH+:");
	PrintE93(ioQQQ, coldenCom.colden[IPCHEHP-1] );
	fprintf( ioQQQ, "\n");

	/* print molecular element column densities */
	molcol("PRIN");

	/* find t^2 from H II structure */
	gett2(&peimbt.t2hstr);

	/* find t^2 from OIII structure */
	gett2o3(&peimbt.t2o3str);

	fprintf( ioQQQ, "\n     Col(Heff):      ");
	PrintE93(ioQQQ, coleff);
	fprintf( ioQQQ, "  snd travl time  ");
	PrintE82(ioQQQ, timesc.sound);
	fprintf( ioQQQ, " sec  NeN+dl: ");
	PrintE82(ioQQQ, abrems.dlnenp);
	fprintf( ioQQQ, "  Te-low:");
	PrintE82(ioQQQ, tehigh.tlowst);
	fprintf( ioQQQ, " Te-hi:");
	PrintE82(ioQQQ, tehigh.thist);
	fprintf( ioQQQ, "\n");

	/* following is wl where gas goes thick to brems, in cm */
	if( tffCom.tff > 0. )
	{
		bremtk = RYDLAM*1e-8/tffCom.tff;
	}
	else
	{
		bremtk = 1e30;
	}

	/* this is number of times ionize was called during the model, over number of zones,
	 * a measure of the convergence of the model - includes inital search so worse for
	 * short calculations 
	 * It is a measure of how hard the model was to converge */
	if( ZoneCnt.nzone > 0 )
	{
		znit = (double)(conv.nTotalIoniz)/(double)(ZoneCnt.nzone);
	}
	else
	{
		znit = 0.;
	}

	/* apparent helium abundance */
	fprintf( ioQQQ, " He/Ha:");
	PrintE82( ioQQQ, heabun);

	/* he/h relative to correct helium abundance */
	fprintf(ioQQQ, "  =%7.2f*true  Lthin:",hescal);

	/* wavelength were structure is optically thick to brems absorption */
	PrintE82( ioQQQ, bremtk);
	fprintf(ioQQQ, "  itr/zn:%5.3f",znit);

	/* say whether we used stored opacities (T) or derived them from scratch (F) */
	fprintf(ioQQQ, "  File Opacity:  %1c ",TorF(opac.lgOpacExist));

	/* get harmonic mean HI temperature, as needed for 21 cm spin temperature */
	fprintf(ioQQQ, " <nH/T>:");
	if( !cdTemp( "21cm",0,&THI,"volume" ) )
	{
		fprintf(ioQQQ,"\n>>>>>>>>>>>>>>>>> prtfinal, impossible problem getting 21cm temp.\n");
	}
	PrintE82(ioQQQ, THI);

	/* now convert this HI spin temperature into a brightness temperature */
	THI *= (1. - sexp( TauLines[ipH21cm].TauIn ) );
	fprintf( ioQQQ, " TB21cm");
	PrintE82(ioQQQ, THI);
	fprintf(ioQQQ, "\n");

	/* timescale (sec here) for photoerosion of fe-ni */
	rate = timesc.TimeErode*2e-26;
	if( rate > 1e-30 )
	{
		ferode = 1./rate;
	}
	else
	{
		ferode = 0.;
	}

	/* mean acceleration */
	if( wind.acldr > 0. )
	{
		wind.AccelAver /= wind.acldr;
	}
	else
	{
		wind.AccelAver = 0.;
	}

	/* DMEAN is mean density (gm per cc); mean temp is weighted by mass density */
	wmean = colmas.wmas/colmas.TotMassColl;
	dmean = colmas.TotMassColl/GrainVar.avhden;
	tmean = colmas.tmas/colmas.TotMassColl;
	wmean = colmas.wmas/colmas.TotMassColl;

	/*fprintf( ioQQQ, "   <a>:%8.2e  erdeFe%7.1e  Tcompt%8.2e  Tthr%8.2e  <Tden>: %8.2e  <dens>:%8.2e <Mol>:%8.2e\n", 
	  wind.AccelAver, ferode, timesc.tcmptn, timesc.ttherm, tmean, 
	  dmean, wmean );*/
	fprintf( ioQQQ, "   <a>:");
	PrintE82(ioQQQ , wind.AccelAver);
	fprintf( ioQQQ, "  erdeFe");
	PrintE71(ioQQQ , ferode);
	fprintf( ioQQQ, "  Tcompt");
	PrintE82(ioQQQ , timesc.tcmptn);
	fprintf( ioQQQ, "  Tthr");
	PrintE82(ioQQQ , timesc.ttherm);
	fprintf( ioQQQ, "  <Tden>: ");
	PrintE82(ioQQQ , tmean);
	fprintf( ioQQQ, "  <dens>:");
	PrintE82(ioQQQ , dmean);
	fprintf( ioQQQ, " <Mol>:");
	PrintE82(ioQQQ , wmean);
	fprintf(ioQQQ,"\n");

	/* Jeans length and mass - this is mean over model */
	rjeans = 7.79637 + (log10(tmean) - log10(densty.wmole) - log10(densty.xMassDensity*
	  filfac.FillFac))/2.;

	if( rjeans < 36. )
	{
		rjeans = (double)pow(10.,rjeans);
		/* AJMASS is Jeans mass in solar units */
		ajmass = 3.*log10(rjeans/2.) + log10(4.*PI/3.*densty.xMassDensity*
		  filfac.FillFac) - log10(SOLAR_MASS);
		if( ajmass < 37 )
		{
			ajmass = pow(10.,ajmass);
		}
		else
		{
			ajmass = 0.;
		}
	}
	else
	{
		rjeans = 0.;
		ajmass = 0.;
	}

	/* Jeans length and mass - this is smallest over model */
	ajmin = jmin.ajmmin - log10(SOLAR_MASS);
	if( ajmin < 37 )
	{
		ajmin = pow(10.,ajmin);
	}
	else
	{
		ajmin = 0.;
	}

	/* estimate alpha (o-x) */
	if( rfield.anu[rfield.nflux-1] > 150. )
	{
		/* generate pointers to energies where alpha ox will be evaluated */
		ip2kev = ipoint(147.);
		ip2500 = ipoint(0.3648);

		/* now get fluxes there */
		bottom = (rfield.flux[ip2500-1] + rfield.ConOutNoInter[ip2500-1])*
		  rfield.anu[ip2500-1]/rfield.widflx[ip2500-1];
		top = (rfield.flux[ip2kev-1]+ rfield.ConOutNoInter[ip2kev-1] )*
		  rfield.anu[ip2kev-1]/rfield.widflx[ip2kev-1];

		/* generate alpha ox if denominator is non-zero */
		if( bottom > 1e-30 && top > 1e-30 )
		{
			ratio = log10(top) - log10(bottom);
			if( ratio < 36. && ratio > -36. )
			{
				ratio = pow(10.,ratio);
			}
			else
			{
				ratio = 0.;
			}
		}

		else
		{
			ratio = 0.;
		}

		if( ratio > 0. )
		{
			alfox = log(ratio)/log(rfield.anu[ip2kev-1]/rfield.anu[ip2500-1]);
		}
		else
		{
			alfox = 0.;
		}
	}
	else
	{
		alfox = 0.;
	}

	if( jmin.rjnmin < 37 )
	{
		jmin.rjnmin = (float)pow(10.,jmin.rjnmin);
	}
	else
	{
		jmin.rjnmin = FLT_MAX;
	}

	/*fprintf( ioQQQ, "     Mean Jeans  l(cm)%8.2e  M(sun)%8.2e  smallest - -  len(cm):%8.2e  M(sun):%8.2e Alf(ox-tran): %10.4f\n", 
	  rjeans, ajmass, jmin.rjnmin, ajmin, alfox );*/
	fprintf( ioQQQ, "     Mean Jeans  l(cm)");
	PrintE82(ioQQQ, rjeans);
	fprintf( ioQQQ, "  M(sun)");
	PrintE82(ioQQQ, ajmass);
	fprintf( ioQQQ, "  smallest:     len(cm):");
	PrintE82(ioQQQ, jmin.rjnmin);
	fprintf( ioQQQ, "  M(sun):");
	PrintE82(ioQQQ, ajmin);
	fprintf( ioQQQ, " Alf(ox-tran): %10.4f\n", alfox );

	/* some details about the hydrogen and helium ions */
	fprintf( ioQQQ, 
		" Hatom level%3ld  NHtopoff:%4ld  HatomType:%4.4s  HInducImp %2c"
		"  He sin level:%3ld  He2 level:  %3ld ExecTime %g\n", 
	  nhlevel, HTopOff.nHTopOff, HTopOff.chHTopType, TorF(HPhoto.lgHInducImp), 
	  NHE1LVL, nhlevel , cdExecTime() );

	/* now give an indication of the convergence error budget */
	fprintf( ioQQQ, 
		" ConvrgError(%%)  <eden>%7.3f  MaxEden%7.3f  <H-C>%7.3f  Max(H-C)   %5.2f  <Press>%8.3f MaxPres%7.3f\n", 
		conv.AverEdenError/ZoneCnt.nzone*100. , 
		conv.BigEdenError*100. ,
		conv.AverHeatCoolError/ZoneCnt.nzone*100. , 
		conv.BigHeatCoolError*100. ,
		conv.AverPressError/ZoneCnt.nzone*100. , 
		conv.BigPressError*100. );

	/* print out some average quantities */
	aver("prin",0.,0.,"          ");

	/* print out Peimbert analysis, tsqden default 1e7, changed
	 * with set tsqden command */
	if( phycon.hden < peimbt.tsqden )
	{
		fprintf(  ioQQQ, " \n" );

		/* temperature from the [OIII] 5007/4363 ratio */
		fprintf(  ioQQQ, " Peimbert T(OIIIr)");
		PrintE82( ioQQQ , peimbt.toiiir);

		/* temperature from predicted Balmer jump relative to Hbeta */
		fprintf(  ioQQQ, " T(Bac)");
		PrintE82( ioQQQ , peimbt.tbac);

		/* temperature predicted from optically thin Balmer jump rel to Hbeta */
		fprintf(  ioQQQ, " T(Hth)");
		PrintE82( ioQQQ , peimbt.tbcthn);

		/* t^2 predicted from the structure, weighted by H */
		fprintf(  ioQQQ, " t2(Hstrc)");
		fprintf( ioQQQ,PrintEfmt("%9.2e", peimbt.t2hstr));

		/* temperature from both [OIII] and the balmer jump rel to Hbeta */
		fprintf(  ioQQQ, " T(O3-BAC)");
		PrintE82( ioQQQ , peimbt.tohyox);

		/* t2 from both [OIII] and the balmer jump rel to Hbeta */
		fprintf(  ioQQQ, " t2(O3-BC)");
		fprintf( ioQQQ,PrintEfmt("%9.2e", peimbt.t2hyox));

		/* sructural t2 from the O+2 predicted radial dependence */
		fprintf(  ioQQQ, " t2(O3str)");
		fprintf( ioQQQ,PrintEfmt("%9.2e", peimbt.t2o3str));

		fprintf(  ioQQQ, "\n");

		if( GrainVar.lgDustOn )
		{
			fprintf( ioQQQ, " Be careful: grains exist.  This spectrum was not corrected for reddening before analysis.\n" );
		}
	}

	/* print average (over radius) grain dust parameters if lgDustOn is true,
	 * avergaeg queantitles are incremented in tauinc, zeroed in startiter */
	if( GrainVar.lgDustOn )
	{
		/*lint -e785 */
		float avdust[NDUST]={FLT_MAX} , 
			avdft[NDUST]={FLT_MAX} , 
			avdpot[NDUST]={FLT_MAX};
		/*lint +e785 */
		nd = 0;
		for( i=0; i < NDUST; i++ )
		{
			if( GrainVar.lgDustOn1[i] )
			{
				/* only save numbers for species that are turned on */
				if( GrainVar.lgQHeat[i] )
				{
					chSTyp[nd] = '*';
				}
				else
				{
					chSTyp[nd] = ' ';
				}

				/* >>chng 97 feb 28, only print out needed ones */
				/* >>chng 99 may 02, use local variable to keep integrety of
				 * of original vars, so that they can later be used in asserts */
				avdust[nd] = GrainVar.avdust[i]/GrainVar.reftot;
				avdft[nd] = GrainVar.avdft[i]/GrainVar.reftot;
				avdpot[nd] = GrainVar.avdpot[i]/GrainVar.reftot;
				strcpy( chLab[nd], chGrainVar.chDstLab[i] );
				++nd;
			}
		}

		fprintf( ioQQQ, "\n   Average grain Properties: \n" );

		fprintf( ioQQQ, "            " );
		for( i=0; i < nd; i++ )
		{
			fprintf( ioQQQ, "%10.10s%c", chLab[i], chSTyp[i] );
		}
		fprintf( ioQQQ, " \n" );

		fprintf( ioQQQ, "   <Tdust>:" );
		for( i=0; i < nd; i++ )
		{
			fprintf( ioQQQ, "  ");
			fprintf( ioQQQ,PrintEfmt("%9.2e", avdust[i]));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "   <Vel D>:" );
		for( i=0; i < nd; i++ )
		{
			fprintf( ioQQQ, "  ");
			fprintf( ioQQQ,PrintEfmt("%9.2e", avdft[i]));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "   <Pot D>:" );
		for( i=0; i < nd; i++ )
		{
			fprintf( ioQQQ, "  ");
			fprintf( ioQQQ,PrintEfmt("%9.2e", avdpot[i]));
		}
		fprintf( ioQQQ, "\n" );
	}

	/* now release saved arrays */
	free( linsav );
	free( sclsav );
	free( lgPrt );
	free( chSTyp) ;

	/* now return the space claimed for the chSLab array*/
	for( i=0; i<LineSave.nsum; ++i)
	{
		free( chSLab[i] );
	}
	free(chSLab );

	free( scaled );
	free( asumli );

	/* option to make short printout */
	if( !pshort.lgPrtShort && called.lgTalk )
	{
		/* print log of optical depths, 
		 * calls prtmet if print line optical depths command entered */
		PrtAllTau();
		/* only print mean ionization and emergent continuum on last iteration */
		if( IterCnt.lgLastIt )
		{
			/* print mean ionization fractions for all elements */
			pmeani('i', ioQQQ);
			/* print mean temperature fractions for all elements */
			pmeani('t', ioQQQ);
			/* print emergent continuum */
			pcontn();
		}
	}

	/* print input title for model */
	fprintf( ioQQQ, "%s\n\n", input.chTitle );

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