/*PrtHeader print program's header, including luminosities and ionization parameters */
#include "cddefines.h"
#include "physconst.h"
#include "phycon.h"
#include "bounds.h"
#include "prhd.h"
#include "totlum.h"
#include "u.h"
#include "he.h"
#include "rfield.h"
#include "qs.h"
#include "fourpi.h"
#include "called.h"
#include "nh.h"
#include "compton.h"
#include "only.h"
#include "tenden.h"
#include "edgtcm.h"
#include "ffun.h"
#include "ipoint.h"
#include "printe82.h"
#include "prtheader.h"

void PrtHeader(void)
{
	long int i, 
	  ip2500, 
	  ip2kev;
	double absbol, 
	  absv, 
	  alfox, 
	  avpow, 
	  bolc, 
	  fbeta, 
	  fluxv, 
	  gpowl, 
	  pballog, 
	  pionl, 
	  qballog, 
	  qgaml, 
	  qheiil, 
	  qhel, 
	  ql, 
	  qxl, 
	  radpwl, 
	  ratio, 
	  solar, 
	  tcomp, 
	  tradio;

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

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

	fprintf( ioQQQ, "           %4ldCellPeak",rfield.nflux );
	PrintE82(ioQQQ, rfield.anu[prhd.ipeak-1] );
	fprintf( ioQQQ, "   Lo");
	fprintf( ioQQQ,PrintEfmt("%9.2e", rfield.anu[0] - rfield.widflx[0]/2. ));
	fprintf( ioQQQ, "=%6.2fcm   Hi-Con:",
	  9.117e-6/(rfield.anu[0] - rfield.widflx[0]/2.) );
	PrintE82(ioQQQ,rfield.anu[rfield.nflux-1] + rfield.widflx[rfield.nflux-1]/2.);
	fprintf(ioQQQ," Ryd   E(hi):");
	PrintE82(ioQQQ,bounds.egamry);
	fprintf(ioQQQ,"Ryd     E(hi):  %9.2f MeV\n", bounds.egamry*0.0000136 );

	if( prhd.xpow > 0. )
	{
		prhd.xpow = (float)(log10(prhd.xpow) + fourpi.pirsq);
		qxl = log10(prhd.qx) + fourpi.pirsq;
	}
	else
	{
		prhd.xpow = 0.;
		qxl = 0.;
	}

	if( prhd.powion > 0. )
	{
		pionl = log10(prhd.powion) + fourpi.pirsq;
		avpow = prhd.powion/qs.qhtot/EN1RYD;
	}
	else
	{
		pionl = 0.;
		avpow = 0.;
	}

	/* >>chng 97 mar 18, break these two out here, so that returns zero
	 *when no radiation - had been done in the print statement so pirsq was printed */
	if( prhd.pbal > 0. )
	{
		pballog = log10(MAX2(1e-30,prhd.pbal)) + fourpi.pirsq;
		qballog = log10(MAX2(1e-30,qs.qbal)) + fourpi.pirsq;
	}
	else
	{
		pballog = 0.;
		qballog = 0.;
	}

	if( fourpi.pirsq > 0. )
	{
		fprintf( ioQQQ, "           L(nu>1ryd):%9.4f   Average nu:",pionl);
		PrintE93(ioQQQ, avpow);
		fprintf( ioQQQ,"   L( X-ray):%9.4f   L(BalC):%9.4f     Q(Balmer C):%9.4f\n", 
		  prhd.xpow, pballog, qballog );
	}
	else
	{
		fprintf( ioQQQ, "           I(nu>1ryd):%9.4f   Average nu:",pionl);
		PrintE93(ioQQQ, avpow);
		fprintf( ioQQQ,"   I( X-ray):%9.4f   I(BalC):%9.4f     Phi(BalmrC):%9.4f\n", 
		  prhd.xpow, 
		  log10(MAX2(1e-30,prhd.pbal)), 
		  log10(MAX2(1e-30,qs.qbal)) );
	}

	if( qs.qhe > 0. )
	{
		qhel = log10(qs.qhe) + fourpi.pirsq;
	}
	else
	{
		qhel = 0.;
	}

	if( qs.qheii > 0. )
	{
		qheiil = log10(qs.qheii) + fourpi.pirsq;
	}
	else
	{
		qheiil = 0.;
	}

	if( prhd.q > 0. )
	{
		ql = log10(prhd.q) + fourpi.pirsq;
	}
	else
	{
		ql = 0.;
	}

	if( fourpi.pirsq != 0. )
	{
		fprintf( ioQQQ, 
			"           Q(1.0-1.8):%9.4f   Q(1.8-4.0):%9.4f   Q(4.0-20):"
			"%9.4f   Q(20--):%9.4f     Ion pht flx:", 
		  ql, 
		  qhel, 
		  qheiil, 
		  qxl);
		PrintE93(ioQQQ, qs.qhtot );
	}
	else
	{
		fprintf( ioQQQ, 
			"           phi(1.0-1.8):%7.4f   phi(1.8-4.0):%7.3f   phi(4.0-20):"
			"%7.3f   phi(20--):%7.3f     Ion pht flx:", 
		  ql, 
		  qhel, 
		  qheiil, 
		  qxl );
		PrintE93(ioQQQ, qs.qhtot );
	}
	fprintf( ioQQQ,"\n");

	/* estimate alpha (o-x) */
	if( rfield.anu[rfield.nflux-1] > 150. )
	{
		/* the ratio of fluxes is nearly 403.3^alpha ox */
		ip2kev = ipoint(147.);
		ip2500 = ipoint(0.3645);
		if( rfield.flux[ip2500-1] > 1e-28 )
		{
			ratio = (rfield.flux[ip2kev-1]*rfield.anu[ip2kev-1]/rfield.widflx[ip2kev-1])/
			  (rfield.flux[ip2500-1]*rfield.anu[ip2500-1]/rfield.widflx[ip2500-1]);
		}
		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( prhd.GammaLumin > 0. )
	{
		gpowl = log10(prhd.GammaLumin) + fourpi.pirsq;
		qgaml = log10(prhd.qgam) + fourpi.pirsq;
	}
	else
	{
		gpowl = 0.;
		qgaml = 0.;
	}

	if( prhd.pradio > 0. )
	{
		radpwl = log10(prhd.pradio) + fourpi.pirsq;
	}
	else
	{
		radpwl = 0.;
	}

	if( fourpi.pirsq > 0. )
	{
		fprintf( ioQQQ, 
			"           L(gam ray):%9.4f   Q(gam ray):%9.4f   L(Infred):%9.4f   Alf(ox):%9.4f     Total lumin:%9.4f\n", 
		  gpowl, qgaml, radpwl, alfox, log10(totlum.TotalLumin) + 
		  fourpi.pirsq );
	}
	else
	{
		fprintf( ioQQQ, 
			"           I(gam ray):%9.4f   phi(gam r):%9.4f   I(Infred):%9.4f   Alf(ox):%9.4f     Total inten:%9.4f\n", 
		  gpowl, qgaml, radpwl, alfox, log10(totlum.TotalLumin) + 
		  fourpi.pirsq );
	}

	/* magnitudes */
	if( fourpi.pirsq > 0. )
	{
		solar = log10(totlum.TotalLumin) + fourpi.pirsq - 33.5828;
		absbol = 4.75 - 2.5*solar;

		/* flux density in V, erg / s / cm2 / hz */
		fluxv = ffun(0.1643)*HPLANCK*0.1643;

		/* absv = 4.79 - 2.5 * (LOG10(MAX(1e-30,fluxv)) + pirsq - 18.742 -
		 *  1  0.016)
		 * allen 76, page 197 */
		absv = -2.5*(log10(MAX2(1e-30,fluxv)) + fourpi.pirsq - 20.654202);

		fbeta = ffun(0.1875)*HPLANCK*0.1875*6.167e14;
		/* >>chng 97 mar 18, following branch so zero returned when no radiation at all */
		if( fbeta > 0. )
		{
			fbeta = log10(MAX2(1e-37,fbeta)) + fourpi.pirsq;
		}
		else
		{
			fbeta = 0.;
		}

		bolc = absbol - absv;
		fprintf( ioQQQ, 
			"           log L/Lsun:%9.4f   Abs bol mg:%9.4f   Abs V mag:%9.4f   Bol cor:%9.4f     nuFnu(Bbet):%9.4f\n", 
		  solar, 
		  absbol, 
		  absv, 
		  bolc, 
		  fbeta );
	}

	Compton.cmcool = 0.;
	Compton.cmheat = 0.;

	for( i=0; i < rfield.nflux; i++ )
	{
		/* CSIGC is Tarter expression times ANU(I)*3.858E-25 */
		Compton.cmcool += (rfield.flux[i] + rfield.outlin[i] + rfield.outcon[i])*
		  Compton.csigc[i];

		/* Compton heating with correction for induced Compton scattering
		 * CSIGH is Tarter expression times ANU(I)**2 * 3.858E-25 */
		Compton.cmheat += (rfield.flux[i] + rfield.outlin[i] + rfield.outcon[i])*
		  Compton.csigh[i]*(1. + rfield.occnum[i]);
	}

	Compton.cmheat *= phycon.eden*1e-15;

	/* 6.338E-6 is k in inf mass Rydbergs */
	Compton.cmcool *= phycon.eden*4.*6.338e-6*1e-15;

	/* all of following need factor of 10^-15 to be true rates */
	if( Compton.cmcool > 0. )
	{
		Compton.lgComUndr = FALSE;
		tcomp = Compton.cmheat/Compton.cmcool;
	}
	else
	{
		/* underflow if Compt cooling rate is zero */
		Compton.lgComUndr = TRUE;
		tcomp = 0.;
	}

	/* check whether energy density temp is greater than compton temp */
	if( tenden.TEnerDen > 1.05*tcomp )
	{
		edgtcm.lgEdnGTcm = TRUE;
	}
	else
	{
		edgtcm.lgEdnGTcm = FALSE;
	}

	/* say some ionization parameters */
	fprintf( ioQQQ, "           U(1.0----):");
	PrintE93( ioQQQ, u.uh);
	fprintf( ioQQQ, "   U(4.0----):");
	PrintE93( ioQQQ,u.uheii );
	fprintf( ioQQQ, "   T(En-Den):");
	PrintE93( ioQQQ,tenden.TEnerDen );
	fprintf( ioQQQ, "   T(Comp):");
	PrintE93( ioQQQ,tcomp );
	fprintf( ioQQQ, "     nuJnu(912A):");
	PrintE93( ioQQQ,prhd.fx1ryd );
	fprintf( ioQQQ, "\n");

	/* some occupation numbers */
	fprintf( ioQQQ, "           Occ(FarIR):");
	PrintE93( ioQQQ, rfield.occnum[0]);
	fprintf( ioQQQ, "   Occ(H n=6):");
	PrintE93( ioQQQ, rfield.occnum[nh.ipHn[0][6]-1]);
	fprintf( ioQQQ, "   Occ(1Ryd):");
	PrintE93( ioQQQ, rfield.occnum[nh.ipHn[0][IP1S]-1]);
	fprintf( ioQQQ, "   Occ(4R):");
	PrintE93( ioQQQ, rfield.occnum[he.nheii-1]);
	fprintf( ioQQQ, "     Occ (Nu-hi):");
	PrintE93( ioQQQ, rfield.occnum[rfield.nflux-1] );
	fprintf( ioQQQ, "\n"); 

	/* now print brightness temps */
	tradio = rfield.occnum[0]*TE1RYD*rfield.anu[0];
	fprintf( ioQQQ, "           Tbr(FarIR):");
	PrintE93( ioQQQ, tradio);
	fprintf( ioQQQ, "   Tbr(H n=6):");
	PrintE93( ioQQQ, rfield.occnum[nh.ipHn[0][6]-1]*TE1RYD*rfield.anu[nh.ipHn[0][6]-1]);
	fprintf( ioQQQ, "   Tbr(1Ryd):");
	PrintE93( ioQQQ, rfield.occnum[nh.ipHn[0][IP1S]-1]*TE1RYD*rfield.anu[nh.ipHn[0][0]-1]);
	fprintf( ioQQQ, "   Tbr(4R):");
	PrintE93( ioQQQ, rfield.occnum[he.nheii-1]*TE1RYD*rfield.anu[he.nheii-1]);
	fprintf( ioQQQ, "     Tbr (Nu-hi):");
	PrintE93( ioQQQ, rfield.occnum[rfield.nflux-1]*TE1RYD*rfield.anu[rfield.nflux-1]);
	fprintf( ioQQQ, "\n");

	if( tradio > 1e9 )
	{
		fprintf( ioQQQ, 
			" >>>The radio brightness temperature is very large,%10.2eK at%10.2ecm.  Is this physical???\n", 
		  tradio, 9.115e-6/rfield.anu[0] );
	}

	/* skip a line */
	fprintf( ioQQQ, "  \n" );

	if( only.lgOnlyHead )
	{
		puts( "[Stop in PrtHeader]" );
		exit(1);
	}


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