/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*PrtZone print out individual zone results, call by iter_end_check at very
 * end of zone calculations */
#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 "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"
#include "taulines.h"

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

	DEBUG_ENTRY( "PrtZone()" );

	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.nTrConvg )
	{
		/* 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_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1] ));
		fprintf(ioQQQ,"###\n");

		if( trace.nTrConvg )
		{
			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.nTrConvg )
	{ 
		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 */
	if( !wind.lgStatic() )
	{
		double fac;
		/* find denominator for fractional contributions */
		if( wind.AccelTotalOutward == 0. )
			fac = 1.;
		else
			fac = wind.AccelTotalOutward;
		fprintf( ioQQQ, 
			" Dynamics wind V:%.3e km/s a(grav):%.2e a(tot):%.2e Fr(cont):%6.3f "
			"Fr(line):%6.3f \n",
			wind.windv/1e5 ,
			-wind.AccelGravity,
			wind.AccelTotalOutward,
			wind.AccelCont/ fac, 
			wind.AccelLine/fac );

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

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

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

	hatmic = 0.;
	count_ptr<chem_atom> elHydrogen = unresolved_atom_list[ipHYDROGEN];
	for(mol = 0; mol < mole_global.num_calc; mol++) {
		if (mole_global.list[mol]->parentLabel.empty() && mole_global.list[mol]->nAtom.find(elHydrogen) !=
			 mole_global.list[mol]->nAtom.end())
			hatmic += mole.species[mol].den*mole_global.list[mol]->nAtom[elHydrogen];
	}
	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",findspecieslocal("H-")->den/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",findspecieslocal("H2+")->den/dense.gas_phase[ipHYDROGEN]));
	fprintf( ioQQQ, " H2+ HeH+");
	fprintf(ioQQQ,PrintEfmt("%9.2e",findspecieslocal("HeH+")->den/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_sp[ipH_LIKE][ipHYDROGEN].lgPrtDepartCoef )
	{
		fprintf( ioQQQ, " Hydrogen     " );
		fprintf(ioQQQ,PrintEfmt("%9.2e",  iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].DepartCoef));
		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" );
	}

	fixit(); // add a command line option to activate this
	// print departure coefficients for diatoms species
	for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
	{
		if( 0 )
			(*diatom)->H2_PrtDepartCoef();
	}

	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");
	}

	/* convert energy flux [erg cm-2 s-1] in attenuated incident continuum, 
	 * diffuse emission, into equivalent temperature */
	double TconTot = pow((rfield.EnergyIncidCont+rfield.EnergyDiffCont)/SPEEDLIGHT/7.56464e-15,0.25);
	double Tincid = pow(rfield.EnergyIncidCont/SPEEDLIGHT/7.56464e-15,0.25);
	double Tdiff = pow(rfield.EnergyDiffCont/SPEEDLIGHT/7.56464e-15,0.25);

	/* convert sum of flux into energy density, then equivalent pressure */
	double Pcon_nT = (Tincid + Tdiff) / SPEEDLIGHT;
	Pcon_nT /= BOLTZMANN;

	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, " HeI 2s3S");
		fprintf(ioQQQ,PrintEfmt("%9.2e", 
			iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe2s3S].Pop()/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_sp[ipH_LIKE][ipHELIUM].lgPrtDepartCoef )
		{
			fprintf( ioQQQ, " Helium       " );
			fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipHE_LIKE][ipHELIUM].fb[0].DepartCoef));
			fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].DepartCoef));
			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 relative to total helium gas phase density */
			double fac = 1./dense.gas_phase[ipHELIUM];
			fprintf( ioQQQ, " He singlet n " );
			fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe1s1S].Pop()*fac ));
			/* singlet n=2 complex */
			fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe2s1S].Pop()*fac ));
			fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe2p1P].Pop()*fac ));
			/* singlet n=3 complex */
			fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe3s1S].Pop()*fac ));
			fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe3p1P].Pop()*fac ));
			fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe3d1D].Pop()*fac ));

			fprintf( ioQQQ, " He tripl" );
			/* triplet n=2 complex */
			fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe2s3S].Pop()*fac ));
			fprintf(ioQQQ,PrintEfmt("%9.2e", 
				iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe2p3P0].Pop()*fac+
				iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe2p3P1].Pop()*fac+
				iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe2p3P2].Pop()*fac ));
			/* triplet n=3 complex */
			fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe3s3S].Pop()*fac ));
			fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe3p3P].Pop()*fac ));
			fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe3d3D].Pop()*fac ));
			fprintf( ioQQQ, "\n" );
		}
	}

	/* loop over iso sequences to see if any populations
	 * and/or departure coefficients need to be printed */
	for( long ipISO = ipH_LIKE; ipISO < NISO; ipISO++ )
	{
		for( nelem=ipISO; nelem<LIMELM; ++nelem )
		{
			if( dense.lgElmtOn[nelem] )
			{
				if( iso_sp[ipISO][nelem].lgPrtLevelPops )
				{
					iso_prt_pops(ipISO, nelem, false);
				}
				if( iso_sp[ipISO][nelem].lgPrtDepartCoef )
				{
					/* true says print departure coefficients
					 * instead of populations. */
					iso_prt_pops(ipISO, nelem, true);
				}
			}
		}
	}

	/* >>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.PresTotlCurr));
	fprintf( ioQQQ, " P( gas )");
	fprintf(ioQQQ,PrintEfmt("%9.2e", pressure.PresGasCurr));
	fprintf( ioQQQ, " P(Radtn)");
	fprintf(ioQQQ,PrintEfmt("%9.2e", pressure.pres_radiation_lines_curr));
	fprintf( ioQQQ, " Rad accl");
	fprintf(ioQQQ,PrintEfmt("%9.2e", wind.AccelTotalOutward));
	fprintf( ioQQQ, " ForceMul");
	fprintf(ioQQQ,PrintEfmt("%9.2e", wind.fmul));
	fprintf( ioQQQ, "\n" );

	fprintf( ioQQQ , " Texc(La)     ");
	fprintf(ioQQQ,PrintEfmt("%9.2e",  hydro.TexcLya ));
	/* attenuated incident continuum */
	fprintf( ioQQQ , " T(I con)");
	fprintf(ioQQQ,PrintEfmt("%9.2e",  Tincid ));
	fprintf( ioQQQ , " T(D con)");
	fprintf(ioQQQ,PrintEfmt("%9.2e",  Tdiff ));
	fprintf( ioQQQ , " T(U tot)");
	fprintf(ioQQQ,PrintEfmt("%9.2e",  TconTot ));
	/* print the total radiation density expressed as an equivalent gas pressure */
	fprintf( ioQQQ , " nT (c+d)");
	fprintf(ioQQQ,PrintEfmt("%9.2e", Pcon_nT ));
	/* 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( size_t nd=0; nd < gv.bin.size(); 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 save-out of quantum heating data to qheat(), by PvH */

	/* heavy element molecules */
	if( findspecieslocal("CO")->den > 0. )
	{
		fprintf( ioQQQ, " Molecules     CH/Ctot:");
		fprintf(ioQQQ,PrintEfmt("%9.2e", findspecieslocal("CH")->den/dense.gas_phase[ipCARBON]));
		fprintf( ioQQQ, " CH+/Ctot");
		fprintf(ioQQQ,PrintEfmt("%9.2e", findspecieslocal("CH+")->den/dense.gas_phase[ipCARBON]));
		fprintf( ioQQQ, " CO/Ctot:");
		fprintf(ioQQQ,PrintEfmt("%9.2e", findspecieslocal("CO")->den/dense.gas_phase[ipCARBON]));
		fprintf( ioQQQ, " CO+/Ctot");
		fprintf(ioQQQ,PrintEfmt("%9.2e", findspecieslocal("CO+")->den/dense.gas_phase[ipCARBON]));
		fprintf( ioQQQ, " H2O/Otot");
		fprintf(ioQQQ,PrintEfmt("%9.2e", findspecieslocal("H2O")->den/dense.gas_phase[ipOXYGEN]));
		fprintf( ioQQQ, " OH/Ototl");
		fprintf(ioQQQ,PrintEfmt("%9.2e", findspecieslocal("OH")->den/dense.gas_phase[ipOXYGEN]));
		fprintf( ioQQQ, "\n");
	}

	/* information about the large H2 molecule - this just returns if not turned on */
	for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
		(*diatom)->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", findspecieslocal("H2O+")->den/MAX2(1e-35,dense.gas_phase[ipOXYGEN]) ));
	fprintf( ioQQQ, " OH+/Otot" );
	fprintf(ioQQQ,PrintEfmt("%9.2e", findspecieslocal("OH+")->den/ 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", findspecieslocal("O2")->den/MAX2(1e-35,dense.gas_phase[ipOXYGEN])));
	fprintf( ioQQQ, " O2+/Otot" );
	fprintf(ioQQQ,PrintEfmt("%9.2e", findspecieslocal("O2+")->den/ 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();
	return;
}