/* 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 */
/*PrtFinal create PrtFinal pages of printout, emission line intensities, etc */
/*StuffComment routine to stuff comments into the stack of comments, def in lines.h */
/*gett2o3 analyze computed [OIII] spectrum to get t^2 */
/*gett2 analyze computed structure to get structural t^2 */
#include "cddefines.h"
#include "utilitymacros.h"
#include "iso.h"
#include "cddrive.h"
#include "physconst.h"
#include "lines.h"
#include "taulines.h"
#include "warnings.h"
#include "dense.h"
#include "grainvar.h"
#include "h2.h"
#include "hmi.h"
#include "thermal.h"
#include "hydrogenic.h"
#include "jeans_mass.h"
#include "rt.h"
#include "atmdat.h"
#include "timesc.h"
#include "opacity.h"
#include "struc.h"
#include "pressure.h"
#include "converge.h"
#include "geometry.h"
#include "called.h"
#include "iterations.h"
#include "version.h"
#include "colden.h"
#include "input.h"
#include "rfield.h"
#include "radius.h"
#include "peimbt.h"
#include "oxy.h"
#include "ipoint.h"
#include "lines_service.h"
#include "aver.h"
#include "wind.h"
#include "prt.h"

/*gett2o3 analyze computed [OIII] spectrum to get t^2 */
static void gett2o3(float *tsqr);

/*gett2 analyze computed structure to get structural t^2 */
static void gett2(float *tsqr);

/* return is TRUE if calculation ok, false otherwise */
void PrtFinal(void)
{
	short int *lgPrt;
	float *wavelength;
	float *sclsav , *scaled;
	long int *ipSortLines;
	double *xLog_line_lumin ;
	char **chSLab;
	char *chSTyp;
	char chCKey[5], 
	  chGeo[7], 
	  chPlaw[21] ;

	long int 
	  i, 
	  ipNegIntensity[33], 
	  ip2500, 
	  ip2kev, 
	  iprnt, 
	  j, 
	  nd, 
	  nline, 
	  nNegIntenLines;
	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 */
	/*lint -e666 expression with side effects */
	fprintf( ioQQQ, "\f\n");
	REPEAT_16( fprintf( ioQQQ, " "); )
	REPEAT_4( fprintf( ioQQQ, " "); )
	REPEAT_2( fprintf( ioQQQ, " ");)
	fprintf( ioQQQ, " "); /*)*/
	REPEAT_32( fprintf( ioQQQ, "*");)
	fprintf( ioQQQ, "*"); /*)*/
	fprintf( ioQQQ, "> Cloudy %7.7s <", version.chVersion);
	REPEAT_32( fprintf( ioQQQ, "*");)
	fprintf( ioQQQ, "\n" );
	/*lint +e666 */

	/*fprintf( ioQQQ, "\f\n                       *********************************> Cloudy ");
	fprintf( ioQQQ, "%7.7s <********************************\n", version.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]);

		/* copy input to CAPS to make sure hide not on it */
		strcpy( input.chCARDCAPS , input.chCardSav[i] );
		caps( input.chCARDCAPS );

		/* print if not continue or hide */
		/* >>chng 04 jan 21, add hide option */
		if( (strcmp(chCKey,"CONT")!= 0) && !lgMatch( "HIDE" , input.chCARDCAPS) )
		{
			/* 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" );
		}
	}

	/* print log of ionization parameter if greater than zero - U==0 for pdr calcs */
	if( rfield.uh > 0. )
	{
		a = log10(rfield.uh);
	}
	else
	{
		a = -37.;
	}

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

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

	if( warnings.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( warnings.lgCautns )
	{
		fprintf( ioQQQ, 
			"                       >>>>>>>>>> Cautions are present.\n" );
	}

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

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

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

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

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

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

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

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

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

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

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

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

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

	/* emission lines as logs of total luminosity
	 * either per unit inner area (lgPredLumin=.F.), or FOUR PI R SQRD (=.T.) */
	if(  radius.distance > 0. && radius.lgRadiusKnown && prt.lgPrintFluxEarth )
	{
		fprintf( ioQQQ, 
			"                                             Flux observed at the Earth (erg/s/cm^2)\n" );
	}

	else if(  prt.lgSurfaceBrightness )
	{
		fprintf( ioQQQ, 
			"                                             Surface brightness (erg/s/cm^2/" );
		if( prt.lgSurfaceBrightness_SR )
		{
			fprintf( ioQQQ, "sr)\n");
		}
		else
		{
			fprintf( ioQQQ, "srcsec^2)\n");
		}
	}

	else if( radius.lgPredLumin )
	{
		/* 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( geometry.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", 
			  geometry.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[LineSave.ipNormWavL].sumlin;

	/* check that this line has positive intensity */
	/* >>chng 01 auyg 21, final relative intensities are a float (since passed to sorting
	 * routine), so change test for zero to SMALLFLOAT */
	if( ((snorm <= SMALLDOUBLE ) || (LineSave.ipNormWavL < 0)) || (LineSave.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 *
	 ******************************************************************/

	/* >>chng 02 aug 29, from lgHCaseBOK to following - caught by Ryan Porter */
	if( !atmdat.lgHCaseBOK[1][ipHYDROGEN] )
	{
#		define NWLH	21
		/* list of all H case b lines */
		float wlh[NWLH] = {6563.e0f, 4861.e0f, 4340.e0f, 4102.e0f, 1.875e4f, 1.282e4f, 
						   1.094e4f, 1.005e4f, 2.625e4f, 2.166e4f, 1.945e4f, 7.458e4f, 
						   4.653e4f, 3.740e4f, 4.051e4f, 7.458e4f, 3.296e4f, 1216.e0f,
						   1026.e0f, 972.5e0f, 949.7e0f };

		/* table exceeded - kill all H case b predictions*/
		for( i=0; i < LineSave.nsum; i++ )
		{
			/* >>chng 04 jun 21, kill both case a and b at same time,
			 * actually lgHCaseBOK has separate A and B flags, but
			 * this is simpler */
			if( (strcmp( LineSv[i].chALab,"Ca B" )==0) || 
				(strcmp( LineSv[i].chALab,"Ca A" )==0) )
			{
				float errorwave;
				/* this logic must be kept parallel with which H lines are
				 * added as case B in lines_hydro - any automatic hosing of
				 * case b would kill both H I and He II */
				errorwave = WavlenErrorGet( LineSv[i].wavelength );
				for( j=0; j<NWLH; ++j )
				{
					if( fabs(LineSv[i].wavelength-wlh[j] ) <= errorwave )
					{
						LineSv[i].sumlin = 0.;
						break;
					}
				}
			}
		}
	}

	if( !atmdat.lgHCaseBOK[1][ipHELIUM] )
	{
		/* table exceeded - kill all He case b predictions*/
#		define NWLHE	20
		float wlhe[NWLHE] = {1640.e0f, 1215.e0f, 1085.e0f, 1025.e0f, 4686.e0f, 3203.e0f,
						     2733.e0f, 2511.e0f, 1.012e4f, 6560.e0f, 5412.e0f, 4860.e0f,
						     1.864e4f, 1.163e4f, 9345.e0f, 8237.e0f, 303.8e0f, 256.3e0f,
						     243.0e0f, 237.3e0f};
		for( i=0; i < LineSave.nsum; i++ )
		{
			if( (strcmp( LineSv[i].chALab,"Ca B" )==0) || 
				(strcmp( LineSv[i].chALab,"Ca A" )==0) )
			{
				float errorwave;
				/* this logic must be kept parallel with which H lines are
				 * added as case B in lines_hydro - any automatic hosing of
				 * case b would kill both H I and He II */
				errorwave = WavlenErrorGet( LineSv[i].wavelength );
				for( j=0; j<NWLHE; ++j )
				{
					if( fabs(LineSv[i].wavelength-wlhe[j] ) <= errorwave )
					{
						LineSv[i].sumlin = 0.;
						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)<=0 )
	{
		if( dense.lgElmtOn[ipHYDROGEN] )
		{
			/* this is a major logical error if hydrogen is turned on */
			fprintf( ioQQQ, " PrtFinal could not find TOTL 4861 with cdLine.\n" );
			puts( "[Stop in PrtFinal]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* get HeI 5876 */
	if( cdLine("He 1",5876.,&he5876,&absint)<=0 )
	{
		if( dense.lgElmtOn[ipHELIUM] && iso.numLevels[ipHE_LIKE][ipHELIUM] >= 14 )
		{
			/* this is a major logical error if hydrogen is turned on */
			fprintf( ioQQQ, " PrtFinal could not find He 1 5876 with cdLine.\n" );
			puts( "[Stop in PrtFinal]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* get HeII 4686 */
	if( cdLine("He 2",4686.,&he4686,&absint)<=0 )
	{
		if( dense.lgElmtOn[ipHELIUM] && iso.numLevels[ipH_LIKE][ipHELIUM] > 5 )
		{
			/* this is a major logical error if hydrogen is turned on 
			 * and the model atom has enough levels */
			fprintf( ioQQQ, " PrtFinal could not find He 2 4686 with cdLine.\n" );
			puts( "[Stop in PrtFinal]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

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

	if( dense.lgElmtOn[ipHELIUM] )
	{
		hescal = heabun/(dense.gas_phase[ipHELIUM]/dense.gas_phase[ipHYDROGEN]);
	}
	else
	{
		hescal = 0.;
	}

	/* get temperature from [OIII] spectrum, o may be turned off,
	 * but lines still dumped into stack */
	if( cdLine("O  3",5007.,&o5007,&absint)<=0 )
	{
		if( dense.lgElmtOn[ipOXYGEN] )
		{
			/* this is a major logical error if hydrogen is turned on */
			fprintf( ioQQQ, " PrtFinal could not find O  3 5007 with cdLine.\n" );
			puts( "[Stop in PrtFinal]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	if( cdLine("TOTL",4363.,&o4363,&absint)<=0 )
	{
		if( dense.lgElmtOn[ipOXYGEN] )
		{
			/* this is a major logical error if hydrogen is turned on */
			fprintf( ioQQQ, " PrtFinal could not find TOTL 4363 with cdLine.\n" );
			puts( "[Stop in PrtFinal]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* 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/oxy.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 07, reset cs to values at 10^4K,
			 * had been last temp in model */
			oxy.o3cs12 = 2.33239f;
			oxy.o3cs13 = 0.292628f;
			ratio = ratio/1.3333/(oxy.o3cs13/oxy.o3cs12);
			/* ratio of energies and branching ratio for 4363 */
			ratio = ratio/oxy.o3enro/oxy.o3br32;
			peimbt.toiiir = (float)(-oxy.o3ex23/log(ratio));
		}
	}

	else
	{
		peimbt.toiiir = 0.;
	}

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

	/* 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)<=0 )
	{
		fprintf( ioQQQ, " PrtFinal could not find thin 3646 with cdLine.\n" );
		puts( "[Stop in PrtFinal]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( cdLine("Ca B",4861.,&hbcab,&absint)<=0 )
	{
		fprintf( ioQQQ, " PrtFinal could not find Ca B 4861 with cdLine.\n" );
		puts( "[Stop in PrtFinal]" );
		cdEXIT(EXIT_FAILURE);
	}

	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 = (float *)MALLOC( sizeof(float)*(unsigned long)LineSave.nsum );
	if( scaled == NULL )
	{ 
		printf( " not enough memory to allocate scaled in PrtFinal\n" );
		puts( "[Stop in PrtFinal]" );
		cdEXIT(EXIT_FAILURE);
	}

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

	/* 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]" );
		cdEXIT(EXIT_FAILURE);
	}

	for( i=0; i < LineSave.nsum; i++ )
	{

		/* when normalizaton line is off-scale small (generally a model
		 * with mis-set parameters) the scaled intensity can be larger than
		 * a float - this is not actually a problem since the number will
		 * overflow the format and hence be unreadable */
		double scale = LineSv[i].sumlin/snorm*LineSave.ScaleNormLine;
		/* this will become a float, so limit dynamic range */
		scale = MIN2(BIGFLOAT , scale );
		scale = MAX2( -BIGFLOAT , scale );

		/* find logs of ALL line intensities/luminosities */
		scaled[i] = (float)scale;

		if( LineSv[i].sumlin > 0. )
		{
			xLog_line_lumin[i] = log10(LineSv[i].sumlin) + radius.Conv2PrtInten;
		}
		else
		{
			xLog_line_lumin[i] = -38.;
		}
	}

	/* now find which lines to print, which to ignore because they are the wrong type */
	for( i=0; i < LineSave.nsum; i++ )
	{
		/* never print unit normalization check, at least in main list */
		if( strcmp(LineSv[i].chALab,"Unit") == 0 )
		{
			lgPrt[i] = FALSE;
		}
		else if( strcmp(LineSv[i].chALab,"Coll") == 0 && !prt.lgPrnColl )
		{
			lgPrt[i] = FALSE;
		}
		else if( strcmp(LineSv[i].chALab,"Pump") == 0 && !prt.lgPrnPump )
		{
			lgPrt[i] = FALSE;
		}
		else if( strncmp(LineSv[i].chALab,"Inw",3) == 0 && !prt.lgPrnInwd )
		{
			lgPrt[i] = FALSE;
		}
		else if( strcmp(LineSv[i].chALab,"Heat") == 0 && !prt.lgPrnHeat )
		{
			lgPrt[i] = FALSE;
		}
		/* these are diffuse and transmitted continua - normally do not print
		 * nFnu or nInu */
		else if( !prt.lgPrnDiff && 
			( (strcmp(LineSv[i].chALab,"nFnu") == 0) || (strcmp(LineSv[i].chALab,"nInu") == 0) ) )
		{
			lgPrt[i] = FALSE;
		}
		else
		{
			lgPrt[i] = TRUE;
		}
	}

	/* do not print relatively faint lines unless requested */
	nNegIntenLines = 0;

	/* get space for wavelength */
	wavelength = (float *)MALLOC( sizeof(float)*(unsigned long)LineSave.nsum );
	if( wavelength == NULL )
	{ 
		printf( " not enough memory to allocate wavelength in PrtFinal\n" );
		puts( "[Stop in PrtFinal]" );
		cdEXIT(EXIT_FAILURE);
	}

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

	/* 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]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* 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]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* 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]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

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

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

	/* >>chng 01 aug 23, only print on second iteration when inward fractions known */
	if( (gv.lgDustOn && (!geometry.lgSphere)) && (geometry.covrt == 0.) && (iteration > 1))
	{
		long int nPrint=0;

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

		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 >  prt.TooFaint)
			{
				/* the label for the line */
				fprintf( ioQQQ, "%4.4s ",LineDSv[i].chSMDLab ); 

				/* print the wavelength for the line */
				prt_wl( ioQQQ , LineDSv[i].wavelength);

				/* print the log of the intensity/luminosity of the lines, the usual case */
				if( prt.lgPrtLineLog )
				{
					fprintf( ioQQQ, " %7.3f", log10(MAX2(LineDSv[i].smdlin,1e-35)) + radius.Conv2PrtInten );
				}
				else
				{
					/* print linear quantity instead */
					fprintf( ioQQQ, "  %.2e ", LineDSv[i].smdlin + radius.Conv2PrtInten );
				}

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

				/* now end the block with some spaces to set next one off */
				fprintf( ioQQQ, "      " );
				++nPrint;
			}
			if( nPrint==4 )
			{
				fprintf( ioQQQ, "\n " );
				nPrint = 0;
			}
		}
		if( nPrint !=0 )
		{
			fprintf( ioQQQ, "\n" );
		}
		fprintf( ioQQQ, 
			"\n                                                   Intrinsic line intensities\n" );
	}

	/* option to only print brighter lines */
	if( prt.lgFaintOn )
	{
		iprnt = 0;
		for( i=0; i < LineSave.nsum; i++ )
		{
			/* this insanity can happen when arrays overrun */
			ASSERT( iprnt <= i);
			if( scaled[i] >= prt.TooFaint && lgPrt[i] )
			{
				if( prt.lgPrtLineLog )
				{
					xLog_line_lumin[iprnt] = log10(LineSv[i].sumlin) + radius.Conv2PrtInten;
				}
				else
				{
					xLog_line_lumin[iprnt] = LineSv[i].sumlin * pow(10.,radius.Conv2PrtInten);
				}
				sclsav[iprnt] = scaled[i];
				chSTyp[iprnt] = LineSv[i].chSumTyp;
				strcpy( chSLab[iprnt], LineSv[i].chALab );
				wavelength[iprnt] = LineSv[i].wavelength;
				++iprnt;
			}
			else if( -scaled[i] > prt.TooFaint && lgPrt[i] )
			{
				/* negative intensities occur if line absorbs continuum */
				ipNegIntensity[nNegIntenLines] = i;
				nNegIntenLines = MIN2(32,nNegIntenLines+1);
			}
			/* special labels to give id for blocks of lines 
			 * do not add these labels when sorting by wavelength since invalid */
			else if( strcmp( LineSv[i].chALab, "####" ) == 0  &&!prt.lgSortLines )
			{
				strcpy( chSLab[iprnt], LineSv[i].chALab );
				xLog_line_lumin[iprnt] = 0.;
				sclsav[iprnt] = 0.;
				chSTyp[iprnt] = LineSv[i].chSumTyp;
				wavelength[iprnt] = LineSv[i].wavelength;
				++iprnt;
			}
		}
	}

	else
	{
		/* print everything */
		iprnt = LineSave.nsum;
		for( i=0; i < LineSave.nsum; i++ )
		{
			if( strcmp( LineSv[i].chALab, "####" ) == 0  )
			{
				strcpy( chSLab[i], LineSv[i].chALab );
				xLog_line_lumin[i] = 0.;
				sclsav[i] = 0.;
				chSTyp[i] = LineSv[i].chSumTyp;
				wavelength[i] = LineSv[i].wavelength;
			}
			else
			{
				sclsav[i] = scaled[i];
				chSTyp[i] = LineSv[i].chSumTyp;
				strcpy( chSLab[i], LineSv[i].chALab );
				wavelength[i] = LineSv[i].wavelength;
			}
			if( scaled[i] < 0. )
			{
				ipNegIntensity[nNegIntenLines] = i;
				nNegIntenLines = MIN2(32,nNegIntenLines+1);
			}
		}
	}

	/* the number of lines to print better be positive */
	ASSERT( iprnt > 0. );

	/* get space for array of indices for lines, for possible sorting */
	ipSortLines = (long *)MALLOC( sizeof(long)*(unsigned long)iprnt );
	if( ipSortLines == NULL )
	{ 
		printf( " not enough memory to allocate ipSortLines in ipSortLines\n" );
		puts( "[Stop in PrtFinal]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* reorder lines according to wavelength for comparison with spectrum
	 * including writing out the results */
	if( prt.lgSortLines )
	{
		int lgFlag;
		if( prt.lgSortLineWavelength )
		{
			/* first check if wavelength range specified */
			if( prt.wlSort1 >-0.1 )
			{
				j = 0;
				/* skip over those lines not desired */
				for( i=0; i<iprnt; ++i )
				{
					if( wavelength[i]>= prt.wlSort1 && wavelength[i]<= prt.wlSort2 )
					{
						if( j!=i )
						{
							sclsav[j] = sclsav[i];
							chSTyp[j] = chSTyp[i];
							strcpy( chSLab[j], chSLab[i] );
							wavelength[j] = wavelength[i];
						}
						++j;
					}
				}
				iprnt = j;
			}
			/*spsort netlib routine to sort array returning sorted indices */
			spsort(wavelength, 
				   iprnt, 
				  ipSortLines, 
				  /* flag saying what to do - 1 sorts into increasing order, not changing
				   * the original routine */
				  1, 
				  &lgFlag);
			if( lgFlag ) 
				TotalInsanity();
		}
		else if ( prt.lgSortLineIntensity )
		{
			/*spsort netlib routine to sort array returning sorted indices */
			spsort(sclsav, 
				   iprnt, 
				  ipSortLines, 
				  /* flag saying what to do - -1 sorts into decreasing order, not changing
				   * the original routine */
				  -1, 
				  &lgFlag);
			if( lgFlag ) 
				TotalInsanity();
		}
		else
		{
			/* only to keep lint happen, or in case vars clobbered */
			TotalInsanity();
		}
	}
	else
	{
		/* do not sort lines (usual case) simply print in incoming order */
		for( i=0; i<iprnt; ++i )
		{
			ipSortLines[i] = i;
		}
	}

	/*----------------------------------------------------------
	 *
	 * print out all lines which made it through the faint filter,
	 * there are iprnt lines to print */

	if( prt.lgPrtLineArray )
	{
		/* nline will be the number of horizontal lines - 
		 * the 4 represents the 4 columns */
		nline = (iprnt + 3)/4;
	}
	else
		/* this option print a single column of emission lines */
		nline = iprnt;
	{
	}
	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)
		{ 
			/* indix with possibly reorderd set of lines */
			long ip = ipSortLines[j];
			/* this is the actual print statement for the emission line spectrum*/
			if( strcmp( chSLab[ip], "####" ) == 0  )
			{
				/* special labels */
				/*fprintf( ioQQQ, "1111222223333333444444444      " ); */
				/* this string was set in */
				fprintf( ioQQQ, "%s",LineSave.chHoldComments[(int)wavelength[ip]] ); 
			}
			else
			{
				/* the label for the line */
				fprintf( ioQQQ, "%4.4s ",chSLab[ip] ); 

				/* print the wavelength for the line */
				prt_wl( ioQQQ , wavelength[ip]);

				/* print the log of the intensity/luminosity of the lines, the usual case */
				if( prt.lgPrtLineLog )
				{
					fprintf( ioQQQ, " %7.3f", xLog_line_lumin[ip] );
				}
				else
				{
					/* print linear quantity instead */
					fprintf( ioQQQ, "  %.2e ", xLog_line_lumin[ip] );
				}

				/* print scaled intensity with various formats,
				 * depending on order of magnitude.  value is
				 * always the same but the format changes. */
				if( sclsav[ip] < 9999.5 )
				{
					fprintf( ioQQQ, "%9.4f",sclsav[ip] );
				}
				else if( sclsav[ip] < 99999.5 )
				{
					fprintf( ioQQQ, "%9.3f",sclsav[ip] );
				}
				else if( sclsav[ip] < 999999.5 )
				{
					fprintf( ioQQQ, "%9.2f",sclsav[ip] );
				}
				else if( sclsav[ip] < 9999999.5 )
				{
					fprintf( ioQQQ, "%9.1f",sclsav[ip] );
				}
				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( nNegIntenLines > 0 )
	{
		fprintf( ioQQQ, " Lines with negative intensities -  Linear itensities relative to normalize line.\n" );
		fprintf( ioQQQ, "          " );

		for( i=0; i < nNegIntenLines; ++i )
		{
			fprintf( ioQQQ, "%ld %s %.0f %.1e, ", 
			  ipNegIntensity[i], 
			  LineSv[ipNegIntensity[i]].chALab, 
			  LineSv[ipNegIntensity[i]].wavelength, 
			  scaled[ipNegIntensity[i]] );
		}
		fprintf( ioQQQ, "\n" );
	}

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

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

	/* 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)thermal.totcol;
		if( (a >= 0.05) && LineSv[i].chSumTyp == 'c' )
		{
			ipNegIntensity[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 ", 
				LineSv[ipNegIntensity[i]].chALab);

			prt_wl( ioQQQ, LineSv[ipNegIntensity[i]].wavelength );

			fprintf( ioQQQ, ":%5.3f", 
				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)thermal.power;
		if( (a >= 0.05) && LineSv[i].chSumTyp == 'h' )
		{
			ipNegIntensity[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 ", 
				LineSv[ipNegIntensity[i]].chALab);

			prt_wl(ioQQQ, LineSv[ipNegIntensity[i]].wavelength);

			fprintf( ioQQQ, ":%5.3f", 
				d[i] );
		}
		fprintf( ioQQQ, "  \n" );
	}

	/* print out ionization parameters and radiation making it through */
	qlow = 0.;
	plow = 0.;
	for( i=0; i < (iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - 1); i++ )
	{
		/* N.B. in following en1ryd prevents overflow */
		plow += (rfield.flux[i] + rfield.ConInterOut[i]+ rfield.outlin[i] + rfield.outlin_noplot[i])*
		  EN1RYD*rfield.anu[i];
		qlow += rfield.flux[i] + rfield.ConInterOut[i]+ rfield.outlin[i] + rfield.outlin_noplot[i];
	}

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

	qion = 0.;
	pion = 0.;
	for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nflux; i++ )
	{
		/* N.B. in following en1ryd prevents overflow */
		pion += (rfield.flux[i] + rfield.ConInterOut[i]+ rfield.outlin[i] + rfield.outlin_noplot[i])*
		  EN1RYD*rfield.anu[i];
		qion += rfield.flux[i] + rfield.ConInterOut[i]+ rfield.outlin[i] + rfield.outlin_noplot[i];
	}

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

	if( pion > 0. )
	{
		qion = log10(qion) + radius.Conv2PrtInten;
		pion = log10(pion) + radius.Conv2PrtInten;
	}
	else
	{
		qion = 0.;
		pion = 0.;
	}

	/* derive ionization parameter for spherical geometry */
	if( rfield.qhtot > 0. )
	{
		if( rfield.lgUSphON )
		{
			/* RSTROM is stromgren radius */
			usp = rfield.qhtot/POW2(rfield.rstrom/radius.rinner)/
			  2.998e10/dense.gas_phase[ipHYDROGEN];
			usp = log10(usp);
		}
		else
		{
			/* no stromgren radius found, use outer radius */
			usp = rfield.qhtot/radius.r1r0sq/2.998e10/dense.gas_phase[ipHYDROGEN];
			usp = log10(usp);
		}
	}

	else
	{
		usp = -37.;
	}

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

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

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

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

	if( thermal.totcol > 0. )
	{
		thermal.totcol = log10(thermal.totcol) + radius.Conv2PrtInten;
	}
	else
	{
		thermal.totcol = 0.;
	}

	if( thermal.FreeFreeTotHeat > 0. )
	{
		thermal.FreeFreeTotHeat = log10(thermal.FreeFreeTotHeat) + radius.Conv2PrtInten;
	}
	else
	{
		thermal.FreeFreeTotHeat = 0.;
	}

	/* following is recombination line intensity */
	reclin = totlin('r');
	if( reclin > 0. )
	{
		reclin = log10(reclin) + radius.Conv2PrtInten;
	}
	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   P(rad/tot)mx:", 
	  powerl, 
	  thermal.totcol, 
	  a, 
	  reclin, 
	  thermal.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.TauAbsGeo[0][rt.ipxry-1] - rt.tauxry;
	coleff /= 2.14e-22;

	/* 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  Te-low: ");
	PrintE82(ioQQQ, thermal.tlowst);
	fprintf( ioQQQ, "  Te-hi: ");
	PrintE82(ioQQQ, thermal.thist);

	/* this is the intensity of the UV continuum at the illuminated face, relative to the Habing value, as
	 * defined by Tielens & Hollenbach 1985 */
	fprintf( ioQQQ, " G0TH85:");
	PrintE82( ioQQQ, hmi.UV_Cont_rel2_Habing_TH85_face );
	/* this is the intensity of the UV continuum at the illuminated face, relative to the Habing value, as
	 * defined by Draine & Bertoldi 1985 */
	fprintf( ioQQQ, "   G0DB96:");
	PrintE82( ioQQQ, hmi.UV_Cont_rel2_Draine_DB96_face );

	fprintf( ioQQQ, "\n");

	fprintf( ioQQQ, "  Emiss Measure    n(e)n(p) dl       ");
	PrintE93(ioQQQ, colden.dlnenp);
	fprintf( ioQQQ, "  n(e)n(He+)dl         ");
	PrintE93(ioQQQ, colden.dlnenHep);
	fprintf( ioQQQ, "  En(e)n(He++) dl      ");
	PrintE93(ioQQQ, colden.dlnenHepp);
	fprintf( ioQQQ, "\n");

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

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

	/* this is ratio of conv.nTotalIoniz, the number of times ConvBase 
	 * was called during the model, over the number of zones.
	 * This is 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( nzone > 0 )
	{
		znit = (double)(conv.nTotalIoniz)/(double)(nzone);
	}
	else
	{
		znit = 0.;
	}
	/* >>chng 02 jan 09, format from 5.3f to 5.2f */
	fprintf(ioQQQ, "  itr/zn:%5.2f",znit);

	/* sort-of same thing for large H2 molecule - number is number of level pop solns per zone */
	fprintf(ioQQQ, "  H2 itr/zn:%6.2f",H2_itrzn());

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

	/* log of total mass in grams if spherical, or gm/cm2 if plane parallel */
	{
		/* this is mass per unit inner area */
		double xmass = log10( SDIV(dense.xMassTotal) );
		xmass += (float)(1.0992099 + 2.*log10(radius.rinner));
		fprintf(ioQQQ," Mass %.3f",
			xmass);
	}
	fprintf(ioQQQ,"\n");

	/* this is the temperature derived from Lya - 21 cm optical depths */
	/* get harmonic mean HI temperature, as needed for 21 cm spin temperature */
	if( cdTemp( "opti",0,&THI,"volume" ) )
	{
		fprintf(ioQQQ,"\n>>>>>>>>>>>>>>>>> prtfinal, impossible problem getting 21cm opt.\n");
		TotalInsanity();
	}
	fprintf( ioQQQ, "   Temps(21 cm)   T(21cm/Ly a)  ");
	PrintE82(ioQQQ, THI );

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

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

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

	/* get mean 21 cm sping temperature */
	if( cdTemp( "spin",0,&THI,"volume" ) )
	{
		fprintf(ioQQQ,"\n>>>>>>>>>>>>>>>>> prtfinal, impossible problem getting spin temp.\n");
		TotalInsanity();
	}
	fprintf( ioQQQ, "\n        <Tspin>  ");
	PrintE82(ioQQQ, THI);
	fprintf(ioQQQ, "\n");

	/* timescale (sec here) for photoerosion of fe-ni */
	rate = timesc.TimeErode*2e-26;
	if( rate > SMALLFLOAT )
	{
		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 = jeans_mass.wmas/jeans_mass.TotMassColl;
	/* >>chng 02 aug 21, from radius.depth_x_fillfac to intergral of radius times fillfac */
	dmean = jeans_mass.TotMassColl/radius.depth_x_fillfac;
	tmean = jeans_mass.tmas/jeans_mass.TotMassColl;
	wmean = jeans_mass.wmas/jeans_mass.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(dense.wmole) - log10(dense.xMassDensity*
	  geometry.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.*dense.xMassDensity*
		  geometry.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 = jeans_mass.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.anu[ip2500-1]/rfield.widflx[ip2500-1];

		top = rfield.flux[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( jeans_mass.rjnmin < 37 )
	{
		jeans_mass.rjnmin = (float)pow(10.f,jeans_mass.rjnmin);
	}
	else
	{
		jeans_mass.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, jeans_mass.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, jeans_mass.rjnmin);
	fprintf( ioQQQ, "  M(sun):");
	PrintE82(ioQQQ, ajmin);
	fprintf( ioQQQ, " Alf(ox-tran): %10.4f\n", alfox );

	if( prt.lgPrintTime )
	{
		/* print execution time by default */
		alfox = cdExecTime();
	}
	else
	{
		/* flag set false with no time command, so that different runs can
		 * compare exactly */
		alfox = 0.;
	}

	/* some details about the hydrogen and helium ions */
	fprintf( ioQQQ, 
		" Hatom level%3ld  NHtopoff:%4ld  HatomType:%4.4s  HInducImp %2c"
		"  He tot level:%3ld  He2 level:  %3ld ExecTime %g\n", 
	  iso.numLevels[ipH_LIKE][ipHYDROGEN], 
	  iso.nTopOff[ipH_LIKE][ipHYDROGEN], 
	  hydro.chHTopType, 
	  TorF(hydro.lgHInducImp), 
	  iso.numLevels[ipHE_LIKE][ipHELIUM],
	  iso.numLevels[ipH_LIKE][ipHYDROGEN] , 
	  alfox );

	/* now give an indication of the convergence error budget */
	fprintf( ioQQQ, 
		" ConvrgError(%%)  <eden>%7.3f  MaxEden%7.3f  <H-C>%7.2f  Max(H-C)%8.2f  <Press>%8.3f MaxPres%7.3f\n", 
		conv.AverEdenError/nzone*100. , 
		conv.BigEdenError*100. ,
		conv.AverHeatCoolError/nzone*100. , 
		conv.BigHeatCoolError*100. ,
		conv.AverPressError/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( dense.gas_phase[ipHYDROGEN] < 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( gv.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,
	 * average quantitles are incremented in radius_increment, zeroed in IterStart */
	if( gv.lgDustOn && gv.lgGrainPhysicsOn )
	{
		long int i0,
		  i1;
		char chQHeat;
		double AV , AB;
		double total_dust2gas = 0.;

		fprintf( ioQQQ, "\n Average Grain Properties (over radius):\n" );

		for( i0=0; i0 < gv.nBin; i0 += 10 ) 
		{
			if( i0 > 0 )
				fprintf( ioQQQ, "\n" );

			/* this is upper limit to how many grain species we will print across lien */
			i1 = MIN2(i0+10,gv.nBin);

			fprintf( ioQQQ, "       " );
			for( nd=i0; nd < i1; nd++ )
			{
				chQHeat = (char)(gv.bin[nd]->lgEverQHeat ? '*' : ' ');
				fprintf( ioQQQ, "%-12.12s%c", gv.bin[nd]->chDstLab, chQHeat );
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "    nd:" );
			for( nd=i0; nd < i1; nd++ )
			{
				if( nd != i0 ) fprintf( ioQQQ,"   " );
				fprintf( ioQQQ, "%7ld   ", nd );
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, " <Tgr>:" );
			for( nd=i0; nd < i1; nd++ )
			{
				if( nd != i0 ) fprintf( ioQQQ,"   " );
				fprintf( ioQQQ,PrintEfmt("%10.3e", gv.bin[nd]->avdust/radius.depth_x_fillfac));
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, " <Vel>:" );
			for( nd=i0; nd < i1; nd++ )
			{
				if( nd != i0 ) fprintf( ioQQQ,"   " );
				fprintf( ioQQQ,PrintEfmt("%10.3e", gv.bin[nd]->avdft/radius.depth_x_fillfac));
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, " <Pot>:" );
			for( nd=i0; nd < i1; nd++ )
			{
				if( nd != i0 ) fprintf( ioQQQ,"   " );
				fprintf( ioQQQ,PrintEfmt("%10.3e", gv.bin[nd]->avdpot/radius.depth_x_fillfac));
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, " <D/G>:" );
			for( nd=i0; nd < i1; nd++ )
			{
				if( nd != i0 ) fprintf( ioQQQ,"   " );
				fprintf( ioQQQ,PrintEfmt("%10.3e", gv.bin[nd]->avDGRatio/radius.depth_x_fillfac));
				/* add up total dust to gas mass ratio */
				total_dust2gas += gv.bin[nd]->avDGRatio/radius.depth_x_fillfac;
			}
			fprintf( ioQQQ, "\n" );
		}

		fprintf(ioQQQ," Dust to gas ratio (by mass):");
		fprintf(ioQQQ,PrintEfmt("%10.3e", total_dust2gas));

		/* total extinction (conv to mags) at V and B per hydrogen, this includes
		 * forward scattering as an extinction process, so is what would be measured
		 * for a star, but not for an extended source where forward scattering
		 * should be discounted */
		AV = rfield.extin_mag_V_point/colden.colden[ipCOL_HTOT];
		AB = rfield.extin_mag_B_point/colden.colden[ipCOL_HTOT];
		/* print A_V/N(Htot) for point and extended sources */
		fprintf(ioQQQ,", A(V)/N(H)(pnt):%.3e, (ext):%.3e", 
			AV,
			rfield.extin_mag_V_extended/colden.colden[ipCOL_HTOT] );

		/* ratio of total to selective extinction */
		fprintf(ioQQQ,", R:");

		/* grey grains have AB - AV == 0 */
		if( fabs(AB-AV)>SMALLFLOAT )
		{
			fprintf(ioQQQ,"%.3e", AV/(AB-AV) );
		}
		else
		{
			fprintf(ioQQQ,"%.3e", 0. );
		}
		fprintf(ioQQQ," AV(ext):%.3e (pnt):%.3e\n",
			rfield.extin_mag_V_extended, 
			rfield.extin_mag_V_point);
	}

	/* now release saved arrays */
	free( wavelength );
	free( ipSortLines );
	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( xLog_line_lumin );

	/* option to make short printout */
	if( !prt.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( iterations.lgLastIt )
		{
			/* option to print column densities, set with print column density command */
			if( prt.lgPrintColumns )
				PrtColumns(ioQQQ);
			/* print mean ionization fractions for all elements */
			PrtMeanIon('i', FALSE, ioQQQ);
			/* print mean ionization fractions for all elements with density weighting*/
			PrtMeanIon('i', TRUE , ioQQQ);
			/* print mean temperature fractions for all elements */
			PrtMeanIon('t', FALSE , ioQQQ);
			/* print mean temperature fractions for all elements with density weighting */
			PrtMeanIon('t', TRUE , ioQQQ);
			/* print emergent continuum */
			PrtContinuum();
		}
	}

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

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

/* routine to stuff comments into the stack of comments,
 * return is index to this comment */
long int StuffComment( char * chComment )
{
	long int n , i;

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

	/* only do this one time per coreload */
	if( LineSave.ipass == 0 )
	{
		if( LineSave.nComment >= NHOLDCOMMENTS )
		{
			fprintf( ioQQQ, " Too many comments have been entered into line array; increase the value of NHOLDCOMMENTS.\n" );
			puts( "[Stop in StuffComment]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* want this to finally be 33 char long to match format */
#		define NWIDTH 33
		strcpy( LineSave.chHoldComments[LineSave.nComment], chComment );

		/* current length of this string */
		n = (long)strlen( LineSave.chHoldComments[LineSave.nComment] );
		for( i=0; i<NWIDTH-n-1-6; ++i )
		{
			strcat( LineSave.chHoldComments[LineSave.nComment], ".");
		}

		strcat( LineSave.chHoldComments[LineSave.nComment], "..");

		for( i=0; i<6; ++i )
		{
			strcat( LineSave.chHoldComments[LineSave.nComment], " ");
		}
	}

	++LineSave.nComment;

#	ifdef DEBUG_FUN
	fputs( " <->StuffComment()\n", debug_fp );
#	endif
	return( LineSave.nComment-1 ) ;

}

/*gett2 analyze computed structure to get structural t^2 */
static void gett2(float *tsqr)
{
	long int i;

	double tmean;
	double a, 
	  as, 
	  b;

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

	/* get T, t^2 */
	a = 0.;
	b = 0.;

	ASSERT( nzone < struc.nzlim);
	for( i=0; i < nzone; i++ )
	{
		as = (double)(struc.volstr[i])*(double)(struc.hiistr[i])*
		  (double)(struc.ednstr[i]);
		a += (double)(struc.testr[i])*as;
		/* B is used twice below */
		b += as;
	}

	if( b <= 0. )
	{
		*tsqr = 0.;
	}
	else
	{
		/* following is H+ weighted mean temp over vol */
		tmean = a/b;
		a = 0.;

		ASSERT( nzone < struc.nzlim );
		for( i=0; i < nzone; i++ )
		{
			as = (double)(struc.volstr[i])*(double)(struc.hiistr[i])*
			  struc.ednstr[i];
			a += (POW2((double)(struc.testr[i]-tmean)))*as;
		}
		*tsqr = (float)(a/(b*(POW2(tmean))));
	}


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

/*gett2o3 analyze computed [OIII] spectrum to get t^2 */
static void gett2o3(float *tsqr)
{
	long int i;
	double tmean;
	double a, 
	  as, 
	  b;

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

	/* get T, t^2 */	a = 0.;
	b = 0.;
	ASSERT(nzone<struc.nzlim);
	for( i=0; i < nzone; i++ )
	{
		as = (double)(struc.volstr[i])*(double)(struc.o3str[i])*
		  (double)(struc.ednstr[i]);
		a += (double)(struc.testr[i])*as;

		/* B is used twice below */
		b += as;
	}

	if( b <= 0. )
	{
		*tsqr = 0.;
	}

	else
	{
		/* following is H+ weighted mean temp over vol */
		tmean = a/b;
		a = 0.;
		ASSERT( nzone < struc.nzlim );
		for( i=0; i < nzone; i++ )
		{
			as = (double)(struc.volstr[i])*(double)(struc.o3str[i])*
			  struc.ednstr[i];
			a += (POW2((double)(struc.testr[i]-tmean)))*as;
		}
		*tsqr = (float)(a/(b*(POW2(tmean))));
	}


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

