/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * For conditions of distribution and use see copyright notice in license.txt */
/*PrtZone print out individual zone results */
/*prt_H_like_Pops print out hydrogenic populations */
/*prt_H_like_DeparCoef print out hydrogenic populations */
#include "cddefines.h"
#include "physconst.h"
#include "iso.h"
#include "grainvar.h"
#include "pressure.h"
#include "wind.h"
#include "conv.h"
#include "trace.h"
#include "helike.h"
#include "magnetic.h"
#include "dense.h"
#include "called.h"
#include "dynamics.h"
#include "h2.h"
#include "mole.h"
#include "secondaries.h"
#include "opacity.h"
#include "colden.h"
#include "geometry.h"
#include "hmi.h"
#include "rfield.h"
#include "thermal.h"
#include "radius.h"
#include "phycon.h"
#include "abund.h"
#include "hydrogenic.h"
#include "ionbal.h"
#include "elementnames.h"
#include "atomfeii.h"
#include "prt.h"
/*prt_H_like_Pops print out hydrogenic populations */
static void prt_H_like_Pops( long nelem );

/*prt_H_like_DeparCoef print out hydrogenic populations */
static void prt_H_like_DeparCoef( long nelem );

void PrtZone(void)
{
	char chField7[8];
	char chLet, 
	  chQHMark;
	long int i, 
	  ishift, 
	  nelem ,
	  nd,
		mol;
	double cdif, 
	  coninc,
	  con_density,
	  fac, 
	  hatmic;

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

	if( thermal.lgUnstable )
	{
		chLet = 'u';
	}
	else
	{
		chLet = ' ';
	}

	/* middle of zone for printing
	rmidle = radius.Radius - radius.drad*0.5*radius.dRadSign;
	dmidle = radius.depth - radius.drad*0.5; */

	/* option to print single line when quiet but tracing convergence
	 * with "trace convergence" command */
	if( called.lgTalk || trace.lgTrConvg )
	{
		/* print either ####123 or ###1234 */
		if( nzone <= 999 )
		{
			sprintf( chField7, "####%3ld", nzone );
		}
		else
		{
			sprintf( chField7, "###%4ld", nzone );
		}

		fprintf(ioQQQ, " %7.7s %cTe:",chField7, chLet);
		PrintE93(ioQQQ,phycon.te);
		fprintf(ioQQQ," Hden:");
		PrintE93(ioQQQ,dense.gas_phase[ipHYDROGEN]);
		fprintf(ioQQQ," Ne:");
		PrintE93(ioQQQ,dense.eden);
		fprintf(ioQQQ," R:");
		PrintE93(ioQQQ,radius.Radius_mid_zone );
		fprintf(ioQQQ," R-R0:");
		PrintE93(ioQQQ,radius.depth_mid_zone);
		fprintf(ioQQQ," dR:");
		PrintE93(ioQQQ,radius.drad);
		fprintf(ioQQQ," NTR:%3ld Htot:",conv.nPres2Ioniz);
		PrintE93(ioQQQ,thermal.htot);
		fprintf(ioQQQ," T912:");
		fprintf(ioQQQ,PrintEfmt("%9.2e",opac.TauAbsGeo[0][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1] ));
		fprintf(ioQQQ,"###\n");
		  
		if( trace.lgTrConvg )
		{
			fprintf( ioQQQ, " H:%.2e %.2e 2H2/H: %.2e He: %.2e %.2e %.2e\n", 
			  dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN], 
			  dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN], 
			  2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN],
			  dense.xIonDense[ipHELIUM][0]/SDIV(dense.gas_phase[ipHELIUM]), 
			  dense.xIonDense[ipHELIUM][1]/SDIV(dense.gas_phase[ipHELIUM]), 
			  dense.xIonDense[ipHELIUM][2]/SDIV(dense.gas_phase[ipHELIUM])
			  );
		}
	}

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

	/* lgDenFlucOn set to true in zero, only false when variable abundances are on,
	 * lgAbTaON set true when element table used */
	if( !dense.lgDenFlucOn || abund.lgAbTaON )
	{
		fprintf( ioQQQ, " Abun:" );
		for( i=0; i < LIMELM; i++ )
		{
			fprintf( ioQQQ,PrintEfmt("%8.1e", dense.gas_phase[i] ));
		}
		fprintf( ioQQQ, "\n" );
	}

	/*-------------------------------------------------
	 * print wind parameters if windy model */
	fac = wind.windv*dense.gas_phase[ipHYDROGEN]*radius.r1r0sq;
	if( wind.windv != 0. )
	{
		if( wind.AccelTot == 0. )
		{
			fac = 1.;
		}
		else
		{
			fac = wind.AccelTot;
		}
		fprintf( ioQQQ, " WIND; V:");
		fprintf(ioQQQ,PrintEfmt("%10.3e",wind.windv/1e5));
		fprintf( ioQQQ, "km/s G:");
		fprintf(ioQQQ,PrintEfmt("%9.2e",-wind.agrav));
		fprintf( ioQQQ, " Accel:");
		fprintf(ioQQQ,PrintEfmt("%9.2e",wind.AccelTot));
		fprintf( ioQQQ, " Fr(cont):%6.3f Fr(line):%6.3f Fr(dP):%6.3f\n", 
		  wind.AccelCont/ fac, wind.AccelLine/fac, wind.AccelPres/fac );

		/* print line with radiation pressure if significant */
		if( pressure.pbeta > .05 )
			PrtLinePres();

		/* print advection information when velocity nenabledegative */
		if( dynamics.lgAdvection )
			DynaPrtZone();
	}

	/*---------------------------------------------------- */

	hatmic = 0.;
	for (mol = 0; mol < N_H_MOLEC; mol++) {
		hatmic += hmi.Hmolec[mol]*hmi.nProton[mol];
	}
	assert(hatmic > 0.);
	hatmic = (dense.xIonDense[ipHYDROGEN][0] + dense.xIonDense[ipHYDROGEN][1])/hatmic;

	fprintf( ioQQQ, " Hydrogen     ");
	fprintf(ioQQQ,PrintEfmt("%9.2e",dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN]));
	fprintf(ioQQQ,PrintEfmt("%9.2e",dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN]));
	fprintf( ioQQQ, " H+o/Hden");
	fprintf(ioQQQ,PrintEfmt("%9.2e",hatmic ));
	fprintf(ioQQQ,PrintEfmt("%9.2e",hmi.Hmolec[ipMHm]/dense.gas_phase[ipHYDROGEN] ));
	fprintf( ioQQQ, " H-    H2");
	/* this is total H2, the sum of "ground" and excited */
	fprintf(ioQQQ,PrintEfmt("%9.2e",hmi.H2_total/dense.gas_phase[ipHYDROGEN]));
	fprintf(ioQQQ,PrintEfmt("%9.2e",hmi.Hmolec[ipMH2p]/dense.gas_phase[ipHYDROGEN]));
	fprintf( ioQQQ, " H2+ HeH+");
	fprintf(ioQQQ,PrintEfmt("%9.2e",hmi.Hmolec[ipMHeHp]/dense.gas_phase[ipHYDROGEN]));
	fprintf( ioQQQ, " Ho+ ColD");
	fprintf(ioQQQ,PrintEfmt("%9.2e",colden.colden[ipCOL_H0]));
	fprintf(ioQQQ,PrintEfmt("%9.2e",colden.colden[ipCOL_Hp]));
	fprintf( ioQQQ, "\n");

	/* print departure coef if desired */
	if( iso.lgPrtDepartCoef[ipH_LIKE][ipHYDROGEN] )
	{
		fprintf( ioQQQ, " Hydrogen     " );
		fprintf(ioQQQ,PrintEfmt("%9.2e",  iso.DepartCoef[ipH_LIKE][ipHYDROGEN][ipH1s]));
		fprintf(ioQQQ,PrintEfmt("%9.2e",  1.));
		fprintf( ioQQQ, " H+o/Hden");
		fprintf(ioQQQ,PrintEfmt("%9.2e", (dense.xIonDense[ipHYDROGEN][0] + dense.xIonDense[ipHYDROGEN][1])/dense.gas_phase[ipHYDROGEN]));
		fprintf(ioQQQ,PrintEfmt("%9.2e", hmi.hmidep));
		fprintf( ioQQQ, " H-    H2");
		fprintf(ioQQQ,PrintEfmt("%9.2e", hmi.h2dep));
		fprintf( ioQQQ, "      H2+");
		fprintf(ioQQQ,PrintEfmt("%9.2e", hmi.h2pdep));
		fprintf( ioQQQ, "      H3+");
		fprintf(ioQQQ,PrintEfmt("%9.2e",hmi.h3pdep));
		fprintf( ioQQQ, "\n" );
	}

	if( prt.lgPrintHeating )
	{
		fprintf( ioQQQ, "              ");
		fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[0][0]/thermal.htot));
		fprintf( ioQQQ,"         ");
		fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[0][15]/thermal.htot));
		fprintf( ioQQQ,"         ");
		fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[0][16]/thermal.htot));
		fprintf( ioQQQ,"\n");
	}

	/* temperature correspoding to radiation fields */
	coninc = 0.;
	cdif = 0.;
	for( i=0; i < rfield.nflux; i++ )
	{
		/* integrated energy flux, ergs s^-1 cm^-2 */
		coninc += rfield.flux[i]*(rfield.anu[i]*EN1RYD);
		cdif += (rfield.outlin[i] + rfield.outlin_noplot[i] +
			rfield.ConInterOut[i])* (rfield.anu[i]*EN1RYD);
	}
	/* convert flux in attenuated incident continuum, diffuse emission, into equivalent temperature */
	coninc = pow(coninc/SPEEDLIGHT/7.56464e-15,0.25);
	cdif = pow(cdif/SPEEDLIGHT/7.56464e-15,0.25);

	/* convert sum of flux into energy density, then equivalent pressure */
	con_density = (coninc + cdif) / SPEEDLIGHT;
	con_density /= BOLTZMANN;

	/* print info on H-like level populations */
	for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
	{
		if( dense.lgElmtOn[nelem] )
		{
			if( iso.lgPrtLevelPops[ipH_LIKE][nelem] )
			{
				prt_H_like_Pops(nelem);
			}

			if( iso.lgPrtDepartCoef[ipH_LIKE][nelem] )
			{
				prt_H_like_DeparCoef(nelem);
			}
		}
	}

	if( prt.lgPrintHeating )
	{
		fprintf( ioQQQ, "              ");
		fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[0][1]/thermal.htot ));
		fprintf( ioQQQ, "         ");
		fprintf(ioQQQ,PrintEfmt("%9.2e", 0. ));
		fprintf( ioQQQ, " BoundCom");
		fprintf(ioQQQ,PrintEfmt("%9.2e", ionbal.CompRecoilHeatLocal/ thermal.htot));
		fprintf( ioQQQ, "   Extra:");
		fprintf(ioQQQ,PrintEfmt("%9.2e",thermal.heating[0][20]/thermal.htot));
		fprintf( ioQQQ, "   Pairs:");
		fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[0][21]/ thermal.htot ));
		fprintf( ioQQQ,"  H-lines\n");
	}

	/* Helium */
	if( dense.lgElmtOn[ipHELIUM] )
	{
		fprintf( ioQQQ, " Helium       " );
		for( i=0; i < 3; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", dense.xIonDense[ipHELIUM][i]/dense.gas_phase[ipHELIUM]) );
		}

		fprintf( ioQQQ, " He I2SP3");
		fprintf(ioQQQ,PrintEfmt("%9.2e", 
			iso.Pop2Ion[ipHE_LIKE][ipHELIUM][ipHe2s3S]*dense.xIonDense[ipHELIUM][1]/dense.gas_phase[ipHELIUM] ));
		fprintf( ioQQQ, " Comp H,C");
		fprintf(ioQQQ,PrintEfmt("%9.2e", rfield.cmheat ));
		fprintf(ioQQQ,PrintEfmt("%9.2e",  rfield.cmcool*phycon.te));
		fprintf( ioQQQ , " Fill Fac");
		fprintf(ioQQQ,PrintEfmt("%9.2e", geometry.FillFac));
		fprintf( ioQQQ , " Gam1/tot");
		fprintf(ioQQQ,PrintEfmt("%9.2e", hydro.H_ion_frac_photo));
		fprintf( ioQQQ, "\n");

		/* option to print departure coef */
		if( iso.lgPrtDepartCoef[ipH_LIKE][ipHELIUM] )
		{
			fprintf( ioQQQ, " Helium       " );
			fprintf(ioQQQ,PrintEfmt("%9.2e", iso.DepartCoef[ipHE_LIKE][ipHELIUM][0]));
			fprintf(ioQQQ,PrintEfmt("%9.2e", iso.DepartCoef[ipH_LIKE][ipHELIUM][ipH1s]));
			fprintf(ioQQQ,PrintEfmt("%9.2e", 1.));

			fprintf( ioQQQ, " Comp H,C");
			fprintf(ioQQQ,PrintEfmt("%9.2e", rfield.cmheat ));
			fprintf(ioQQQ,PrintEfmt("%9.2e", rfield.cmcool*phycon.te ));
			fprintf( ioQQQ , " Fill Fac");
			fprintf(ioQQQ,PrintEfmt("%9.2e", geometry.FillFac ));
			fprintf( ioQQQ , " Gam1/tot");
			fprintf(ioQQQ,PrintEfmt("%9.2e", hydro.H_ion_frac_photo));
			fprintf( ioQQQ, "\n");
		}

		/* print heating from He (and others) if desired
		 * entry "lines" is induced line heating
		 * 1,12 ffheat:  2,3 he triplets, 1,20 compton */
		if( prt.lgPrintHeating )
		{
			/*fprintf( ioQQQ, "            %10.3e%10.3e    Lines:%10.2e%10.2e  Compton:%10.3e FF Heatig%10.3e\n", 
			  thermal.heating[1][0]/thermal.htot, thermal.heating[1][1]/
			  thermal.htot, thermal.heating[0][22]/thermal.htot, thermal.heating[1][2]/
			  thermal.htot, thermal.heating[0][19]/thermal.htot, thermal.heating[0][11]/
			  thermal.htot );*/
			fprintf( ioQQQ, "              ");
			fprintf(ioQQQ,PrintEfmt("%9.2e",thermal.heating[1][0]/thermal.htot));
			fprintf(ioQQQ,PrintEfmt("%9.2e",thermal.heating[1][1]/thermal.htot));
			fprintf( ioQQQ, "   Lines:");
			fprintf(ioQQQ,PrintEfmt("%9.2e",thermal.heating[0][22]/thermal.htot));
			fprintf(ioQQQ,PrintEfmt("%9.2e",thermal.heating[1][2]/thermal.htot));
			fprintf( ioQQQ, " Compton:");
			fprintf(ioQQQ,PrintEfmt("%9.2e",thermal.heating[0][19]/thermal.htot));
			fprintf( ioQQQ, " FFHeatig");
			fprintf(ioQQQ,PrintEfmt("%9.2e",thermal.heating[0][11]/thermal.htot));
			fprintf( ioQQQ, "\n");
		}

		if( dense.lgElmtOn[ipHELIUM] )
		{
			/* helium singlets and triplets */
			fac = dense.xIonDense[ipHELIUM][1]/dense.gas_phase[ipHELIUM];
			fprintf( ioQQQ, " He singlet n " );
			fprintf(ioQQQ,PrintEfmt("%9.2e", iso.Pop2Ion[ipHE_LIKE][1][ipHe1s1S]*fac ));
			/* singlet n=2 complex */
			if( iso.numPrintLevels[ipHE_LIKE][ipHELIUM]>= ipHe2p1P )
			{
				fprintf(ioQQQ,PrintEfmt("%9.2e", iso.Pop2Ion[ipHE_LIKE][1][ipHe2s1S]*fac ));
				fprintf(ioQQQ,PrintEfmt("%9.2e", iso.Pop2Ion[ipHE_LIKE][1][ipHe2p1P]*fac ));
			}
			else
			{
				fprintf(ioQQQ,PrintEfmt("%9.2e", 0. ));
				fprintf(ioQQQ,PrintEfmt("%9.2e", 0. ));
			}
			/* singlet n=3 complex */
			if( iso.numPrintLevels[ipHE_LIKE][ipHELIUM]>= ipHe3p1P )
			{
				fprintf(ioQQQ,PrintEfmt("%9.2e", iso.Pop2Ion[ipHE_LIKE][ipHELIUM][ipHe3s1S]*fac ));
				fprintf(ioQQQ,PrintEfmt("%9.2e", iso.Pop2Ion[ipHE_LIKE][ipHELIUM][ipHe3p1P]*fac ));
				fprintf(ioQQQ,PrintEfmt("%9.2e", iso.Pop2Ion[ipHE_LIKE][ipHELIUM][ipHe3d1D]*fac ));
			}
			else
			{
				fprintf(ioQQQ,PrintEfmt("%9.2e", 0. ));
				fprintf(ioQQQ,PrintEfmt("%9.2e", 0. ));
				fprintf(ioQQQ,PrintEfmt("%9.2e", 0. ));
			}

			fprintf( ioQQQ, " He tripl" );
			/* triplet n=2 complex */
			if( iso.numPrintLevels[ipHE_LIKE][ipHELIUM]>= ipHe2p3P2 )
			{
				fprintf(ioQQQ,PrintEfmt("%9.2e", iso.Pop2Ion[ipHE_LIKE][1][ipHe2s3S]*fac ));
				fprintf(ioQQQ,PrintEfmt("%9.2e", 
					iso.Pop2Ion[ipHE_LIKE][1][ipHe2p3P0]*fac+
					iso.Pop2Ion[ipHE_LIKE][1][ipHe2p3P1]*fac+
					iso.Pop2Ion[ipHE_LIKE][1][ipHe2p3P2]*fac ));
			}
			else
			{
				fprintf(ioQQQ,PrintEfmt("%9.2e", 0. ));
				fprintf(ioQQQ,PrintEfmt("%9.2e", 0. ));
			}
			/* triplet n=3 complex */
			if( iso.numPrintLevels[ipHE_LIKE][ipHELIUM]> ipHe3d3D )
			{
				fprintf(ioQQQ,PrintEfmt("%9.2e", iso.Pop2Ion[ipHE_LIKE][1][ipHe3s3S]*fac ));
				fprintf(ioQQQ,PrintEfmt("%9.2e", iso.Pop2Ion[ipHE_LIKE][1][ipHe3p3P]*fac ));
				fprintf(ioQQQ,PrintEfmt("%9.2e", iso.Pop2Ion[ipHE_LIKE][1][ipHe3d3D]*fac ));
			}
			else
			{
				fprintf(ioQQQ,PrintEfmt("%9.2e", 0. ));
				fprintf(ioQQQ,PrintEfmt("%9.2e", 0. ));
				fprintf(ioQQQ,PrintEfmt("%9.2e", 0. ));
			}
			fprintf( ioQQQ, "\n" );
		}
	}

	/* loop over all possible he-like elements to see if any need to be printed */
	for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
	{
		if( dense.lgElmtOn[nelem] )
		{
			if( iso.lgPrtLevelPops[ipHE_LIKE][nelem] )
			{
				prt_He_like_Pops(nelem);
			}
			if( iso.lgPrtDepartCoef[ipHE_LIKE][nelem] )
			{
				prt_He_like_DeparCoef(nelem);
			}
		}
	}

	/* >>chng 01 dec 08, move pressure to line before grains, after radiation properties */
	/* gas pressure, pressure due to incident radiation field, rad accel */
	fprintf( ioQQQ, " Pressure      NgasTgas");
	fprintf(ioQQQ,PrintEfmt("%9.2e", pressure.PresGasCurr/BOLTZMANN));
	fprintf( ioQQQ, " P(total)");
	fprintf(ioQQQ,PrintEfmt("%9.2e", pressure.PresRadCurr + pressure.PresGasCurr));
	fprintf( ioQQQ, " P( gas )");
	fprintf(ioQQQ,PrintEfmt("%9.2e", pressure.PresGasCurr));
	fprintf( ioQQQ, " P(Radtn)");
	fprintf(ioQQQ,PrintEfmt("%9.2e", pressure.PresRadCurr));
	fprintf( ioQQQ, " Rad accl");
	fprintf(ioQQQ,PrintEfmt("%9.2e", wind.AccelTot));
	fprintf( ioQQQ, " ForceMul");
	fprintf(ioQQQ,PrintEfmt("%9.2e", wind.fmul));
	fprintf( ioQQQ, "\n" );

	fprintf( ioQQQ , "               Texc(La)");
	fprintf(ioQQQ,PrintEfmt("%9.2e",  hydro.TexcLya ));
	fprintf( ioQQQ , " T(contn)");
	fprintf(ioQQQ,PrintEfmt("%9.2e",  coninc ));
	fprintf( ioQQQ , " T(diffs)");
	fprintf(ioQQQ,PrintEfmt("%9.2e",  cdif ));
	/* print the total radiation density expressed as an equivalent gas pressure */
	fprintf( ioQQQ , " nT (c+d)");
	fprintf(ioQQQ,PrintEfmt("%9.2e", con_density ));
	/* print the radiation to gas pressure */
	fprintf( ioQQQ , " Prad/Gas");
	fprintf(ioQQQ,PrintEfmt("%9.2e", pressure.pbeta ));
	/* magnetic to gas pressure ratio */
	fprintf( ioQQQ , " Pmag/Gas");
	fprintf(ioQQQ,PrintEfmt("%9.2e",  magnetic.pressure / pressure.PresGasCurr) );
	fprintf( ioQQQ, "\n" );

	if( gv.lgGrainPhysicsOn )
	{
		for( nd=0; nd < gv.nBin; nd++ )
		{
			/*  Change things so the quantum heated dust species are marked with an
			*  asterisk just after the name (K Volk)
			*  added QHMARK here and in the write statement */
			chQHMark = (char)(( gv.bin[nd]->lgQHeat && gv.bin[nd]->lgUseQHeat ) ? '*' : ' ');
			fprintf( ioQQQ, "%-12.12s%c  DustTemp",gv.bin[nd]->chDstLab, chQHMark);
			fprintf(ioQQQ,PrintEfmt("%9.2e", gv.bin[nd]->tedust));
			fprintf( ioQQQ, " Pot Volt");
			fprintf(ioQQQ,PrintEfmt("%9.2e", gv.bin[nd]->dstpot*EVRYD));
			fprintf( ioQQQ, " Chrg (e)");
			fprintf(ioQQQ,PrintEfmt("%9.2e", gv.bin[nd]->AveDustZ));
			fprintf( ioQQQ, " drf cm/s");
			fprintf(ioQQQ,PrintEfmt("%9.2e", gv.bin[nd]->DustDftVel));
			fprintf( ioQQQ, " Heating:");
			fprintf(ioQQQ,PrintEfmt("%9.2e", gv.bin[nd]->GasHeatPhotoEl));
			fprintf( ioQQQ, " Frac tot");
			fprintf(ioQQQ,PrintEfmt("%9.2e", gv.bin[nd]->GasHeatPhotoEl/thermal.htot));
			fprintf( ioQQQ, "\n" );
		}
	}
	/* >>chng 00 apr 20, moved punch-out of quantum heating data to qheat(), by PvH */

	/* heavy element molecules */
	if( co.hevmol[ipCO] > 0. )
	{
		fprintf( ioQQQ, " Molecules     CH/Ctot:");
		fprintf(ioQQQ,PrintEfmt("%9.2e", co.hevmol[ipCH]/dense.gas_phase[ipCARBON]));
		fprintf( ioQQQ, " CH+/Ctot");
		fprintf(ioQQQ,PrintEfmt("%9.2e", co.hevmol[ipCHP]/dense.gas_phase[ipCARBON]));
		fprintf( ioQQQ, " CO/Ctot:");
		fprintf(ioQQQ,PrintEfmt("%9.2e", co.hevmol[ipCO]/dense.gas_phase[ipCARBON]));
		fprintf( ioQQQ, " CO+/Ctot");
		fprintf(ioQQQ,PrintEfmt("%9.2e", co.hevmol[ipCOP]/dense.gas_phase[ipCARBON]));
		fprintf( ioQQQ, " H2O/Otot");
		fprintf(ioQQQ,PrintEfmt("%9.2e", co.hevmol[ipH2O]/dense.gas_phase[ipOXYGEN]));
		fprintf( ioQQQ, " OH/Ototl");
		fprintf(ioQQQ,PrintEfmt("%9.2e", co.hevmol[ipOH]/dense.gas_phase[ipOXYGEN]));
		fprintf( ioQQQ, "\n");
	}

	/* information about the large H2 molecule - this just returns if not turned on */
	H2_Prt_Zone();

	/* Lithium, Beryllium */
	if( dense.lgElmtOn[ipLITHIUM] || dense.lgElmtOn[ipBERYLLIUM] || 
		(secondaries.csupra[ipHYDROGEN][0]>0.) )
	{
		fprintf( ioQQQ, " Lithium      " );
		for( i=0; i < 4; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", dense.xIonDense[ipLITHIUM][i]/MAX2(1e-35,dense.gas_phase[ipLITHIUM]) ));
		}
		fprintf( ioQQQ, " Berylliu" );
		for( i=0; i < 5; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e",  dense.xIonDense[ipBERYLLIUM][i]/MAX2(1e-35,dense.gas_phase[ipBERYLLIUM])) );
		}

		/* print secondary ionization rate for atomic hydrogen */
		fprintf( ioQQQ, " sec ion:" );
		fprintf(ioQQQ,PrintEfmt("%9.2e", secondaries.csupra[ipHYDROGEN][0]) );
		fprintf( ioQQQ, "\n" );

		/* option to print heating due to these stages*/
		if( prt.lgPrintHeating )
		{
			fprintf( ioQQQ, "              " );
			for( i=0; i < 3; i++ )
			{
				fprintf(ioQQQ,PrintEfmt("%9.2e",  thermal.heating[ipLITHIUM][i]/ thermal.htot) );
			}
			fprintf( ioQQQ, "                    " );

			for( i=0; i < 4; i++ )
			{
				fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[ipBERYLLIUM][i]/thermal.htot ));
			}
			fprintf( ioQQQ, "\n" );
		}
	}

	/* Boron */
	if( dense.lgElmtOn[ipBORON] )
	{
		fprintf( ioQQQ, " Boron        " );
		for( i=0; i < 6; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", dense.xIonDense[ipBORON][i]/MAX2(1e-35,dense.gas_phase[ipBORON]) ));
		}
		fprintf( ioQQQ, "\n" );

		/* option to print heating*/
		if( prt.lgPrintHeating )
		{ 
			fprintf( ioQQQ, "              " );
			for( i=0; i < 5; i++ )
			{
				fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[ipBORON][i]/thermal.htot ));
			}
			fprintf( ioQQQ, "\n" );
		}
	}

	/* Carbon */
	fprintf( ioQQQ, " Carbon       " );
	for( i=0; i < 7; i++ )
	{
		fprintf(ioQQQ,PrintEfmt("%9.2e", dense.xIonDense[ipCARBON][i]/SDIV(dense.gas_phase[ipCARBON])) );
	}
	/* some molecules trail the line */
	fprintf( ioQQQ, " H2O+/O  " );
	fprintf(ioQQQ,PrintEfmt("%9.2e", co.hevmol[ipH2OP]/MAX2(1e-35,dense.gas_phase[ipOXYGEN]) ));
	fprintf( ioQQQ, " OH+/Otot" );
	fprintf(ioQQQ,PrintEfmt("%9.2e", co.hevmol[ipOHP]/ MAX2(1e-35,dense.gas_phase[ipOXYGEN]) ));
	/* print extra heating, normally zero */
	fprintf( ioQQQ, " Hex(tot)" );
	fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[0][20] ));
	fprintf( ioQQQ, "\n" );

	/* option to print heating*/
	if( prt.lgPrintHeating )
	{
		fprintf( ioQQQ, "              " );
		for( i=0; i < ipCARBON+1; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[ipCARBON][i]/ thermal.htot) );
		}
		fprintf( ioQQQ, "\n" );
	}

	/* Nitrogen */
	fprintf( ioQQQ, " Nitrogen     " );
	for( i=1; i <= 8; i++ )
	{
		fprintf(ioQQQ,PrintEfmt("%9.2e",dense.xIonDense[ipNITROGEN][i-1]/ SDIV(dense.gas_phase[ipNITROGEN]) ));
	}
	fprintf( ioQQQ, " O2/Ototl" );
	fprintf(ioQQQ,PrintEfmt("%9.2e", co.hevmol[ipO2]/MAX2(1e-35,dense.gas_phase[ipOXYGEN])));
	fprintf( ioQQQ, " O2+/Otot" );
	fprintf(ioQQQ,PrintEfmt("%9.2e", co.hevmol[ipO2P]/ MAX2(1e-35,dense.gas_phase[ipOXYGEN]) ));
	fprintf( ioQQQ, "\n" );

	/* option to print heating*/
	if( prt.lgPrintHeating )
	{
		fprintf( ioQQQ, "              " );
		for( i=0; i < ipNITROGEN+1; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[ipNITROGEN][i]/ thermal.htot ));
		}
		fprintf( ioQQQ, "\n" );
	}

#	if 0
	/* Oxygen */
	fprintf( ioQQQ, " Oxygen       " );
	for( i=1; i <= 9; i++ )
	{
		fprintf(ioQQQ,PrintEfmt("%9.2e",dense.xIonDense[ipOXYGEN][i-1]/ SDIV(dense.gas_phase[ipOXYGEN]) ));
	}
	fprintf( ioQQQ, "\n" );

	/* option to print heating*/
	if( prt.lgPrintHeating )
	{
		fprintf( ioQQQ, "              " );
		for( i=0; i < ipOXYGEN+1; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[ipOXYGEN][i]/ thermal.htot ));
		}
		fprintf( ioQQQ, "\n" );
	}
#	endif

	/* now print rest of elements inside loops */
	/* fluorine through Magnesium */
	for( nelem=ipOXYGEN; nelem < ipALUMINIUM; ++nelem )
	{
		if( dense.lgElmtOn[nelem] )
		{
			/* print the element name and amount of shift */
			fprintf( ioQQQ, " %10.10s   ", elementnames.chElementName[nelem]);

			for( i=0; i < nelem+2; i++ )
			{
				fprintf(ioQQQ,PrintEfmt("%9.2e", dense.xIonDense[nelem][i]/dense.gas_phase[nelem] ));
			}
			fprintf( ioQQQ, "\n" );

			/* print heating but only if needed */
			if( prt.lgPrintHeating )
			{
				fprintf( ioQQQ, "              " );
				for( i=0; i < nelem+1; i++ )
				{
					fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[nelem][i]/thermal.htot ));
				}
				fprintf( ioQQQ, "\n" );
			}
		}
	}

	/* Aluminium through Zinc */
	for( nelem=ipALUMINIUM; nelem < LIMELM; ++nelem )
	{
		if( dense.lgElmtOn[nelem] )
		{
			/* number of ionization stages to print across the page */
			/*@-redef@*/
			enum {LINE=13};
			/*@+redef@*/
			ishift = MAX2(0,dense.IonHigh[nelem]-LINE+1);

			/* print the element name and amount of shift */
			fprintf( ioQQQ, " %10.10s%2ld ", elementnames.chElementName[nelem],ishift );

			for( i=0; i < LINE; i++ )
			{
				fprintf(ioQQQ,PrintEfmt("%9.2e", dense.xIonDense[nelem][i+ishift]/dense.gas_phase[nelem]) );
			}
			fprintf( ioQQQ, "\n" );

			/* print heating but only if needed */
			if( prt.lgPrintHeating )
			{
				fprintf( ioQQQ, "              " );
				for( i=0; i < LINE; i++ )
				{
					fprintf(ioQQQ,
						PrintEfmt("%9.2e", thermal.heating[nelem][i+ishift]/thermal.htot ));
				}
				fprintf( ioQQQ, "\n" );
			}
		}
	}

	/* call FeII print routine if large atom is turned on */
	FeIIPrint();

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

/*prt_H_like_Pops print out hydrogenic populations */
static void prt_H_like_Pops( long nelem )
{
	long i ,
		limit ;

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

	/* lgPrintHLevPops is set true with print h-like command, usually false
	 * when set true, also parse in optional number of levels to print, nPrintHLevPops
	 * will be smaller than current value of iso.numLevels_max[ipH_LIKE][nelem] if 
	 * nothing entered, so set levels 
	 * before print command */
	if( dense.IonHigh[nelem] != nelem + 1  )
	{
#		ifdef DEBUG_FUN
		fputs( " <->prt_H_like_Pops()\n", debug_fp );
#		endif
		/* return if no hydrogenic species present */
		return;
	}

	/* it is impossible for this to fail since ParsePrint set to at least 3 */
	ASSERT( iso.numPrintLevels[ipH_LIKE][nelem] > 2);

	/* print straight level pop per ionized hydrogen */
	fprintf( ioQQQ, " %2s%2li 1S-12   ",
		elementnames.chElementSym[nelem],(nelem+1) );

	limit = MIN2(iso.numPrintLevels[ipH_LIKE][nelem],13);
	for( i=ipH1s; i < limit; i++ )
	{
		fprintf( ioQQQ,PrintEfmt("%9.2e",  iso.Pop2Ion[ipH_LIKE][nelem][i] ));
	}
	fprintf(ioQQQ,"\n");

	if( iso.numPrintLevels[ipH_LIKE][nelem] >= limit )
	{
		long int nprint = 0;

		fprintf( ioQQQ, " %2s%2li  rest   ",
		elementnames.chElementSym[nelem], (nelem+1) );

		for( i=limit; i < iso.numPrintLevels[ipH_LIKE][nelem]; i++ )
		{
			fprintf( ioQQQ,PrintEfmt("%9.2e", iso.Pop2Ion[ipH_LIKE][nelem][i] ));
			++nprint;

			if( nprint == 13 )
			{
				nprint = 0;
				/* throw a newline every 14th column */
				fprintf( ioQQQ, "\n" );
				if( i!=iso.numPrintLevels[ipH_LIKE][nelem]-1 )
				{
					/* space out if not first line with label */
					fprintf( ioQQQ, "              " );
				}
			}
		}
		/* end of loop, check whether we need another eol */
		if( nprint != 0 )
		{
			/* do not throw another if iso.numPrintLevels[ipH_LIKE][ipHYDROGEN] is multiple of 13 */
			fprintf( ioQQQ, "\n" );
		}
	}

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

/*prt_H_like_Pops print out hydrogenic populations */
static void prt_H_like_DeparCoef( long nelem )
{
	long i ,
		limit ;

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

	/* lgPrintHLevPops is set true with print h-like command, usually false
	 * when set true, also parse in optional number of levels to print, nPrintHLevPops
	 * will be current value of iso.numLevels_max[ipH_LIKE][nelem] if nothing entered, so set levels
	 * before print command */
	if( dense.IonHigh[nelem] != nelem + 1  )
	{
#		ifdef DEBUG_FUN
		fputs( " <->prt_H_like_DeparCoef()\n", debug_fp );
#		endif
		/* return if no hydrogenic species present */
		return;
	}

	/* it is impossible for this to fail since ParsePrint set to at least 3 */
	ASSERT( iso.numPrintLevels[ipH_LIKE][nelem] > 2);

	limit = MIN2(iso.numPrintLevels[ipH_LIKE][nelem],13);
	/* print this line so that it agrees with level populations above */
	fprintf( ioQQQ, " %2s%2li 1S-12   ",
		elementnames.chElementSym[nelem],(nelem+1) );

	for( i=ipH1s; i < limit; i++ )
	{
		fprintf( ioQQQ,PrintEfmt("%9.2e", iso.DepartCoef[ipH_LIKE][nelem][i]) );
	}
	fprintf( ioQQQ, "\n" );

	if( iso.numPrintLevels[ipH_LIKE][nelem] >= limit )
	{
		long int nprint = 0;

		fprintf( ioQQQ, " %2s%2li  rest   ",
		elementnames.chElementSym[nelem], (nelem+1) );

		for( i=limit; i < iso.numPrintLevels[ipH_LIKE][nelem]; i++ )
		{
			fprintf( ioQQQ,PrintEfmt("%9.2e", iso.DepartCoef[ipH_LIKE][nelem][i] ));
			++nprint;
			if( nprint == 13 )
			{
				/* throw a newline every 13th column */
				fprintf( ioQQQ, "\n" );
				nprint = 0;
				if( i!=iso.numPrintLevels[ipH_LIKE][nelem]-1 )
				{
					/* throw a newline every 13th column but not on the last one */
					fprintf( ioQQQ, "              " );
				}
			}
		}
		if( nprint != 0 )
		{
			/* do not throw another if iso.numPrintLevels[ipH_LIKE][ipHYDROGEN] is multiple of 13 */
			fprintf( ioQQQ, "\n" );
		}
	}

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