/* 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 */
/*HydroLevel calls iso_level to solve for ionization balance 
 * and level populations of model hydrogen atom */
/*PrtHydroTrace1 print trace info for hydrogen-like species */
#include "cddefines.h"
#include "taulines.h"
#include "iso.h"
#include "dense.h"
#include "secondaries.h"
#include "trace.h"
#include "phycon.h"
#include "ionbal.h"
#include "hydrogenic.h"

/*PrtHydroTrace1a print trace info for hydrogen-like species */
STATIC void PrtHydroTrace1a(long nelem )
{
	double colfrc, 
	  phtfrc, 
	  secfrc;

	DEBUG_ENTRY( "PrtHydroTrace1a()" );

	/* identify how atom is ionized for full trace */
	if( iso_sp[ipH_LIKE][nelem].xIonSimple > 0. )
	{
		/* fraction of ionization due to photoionization */
		phtfrc = iso_sp[ipH_LIKE][nelem].fb[ipH1s].gamnc/((dense.eden*(iso_sp[ipH_LIKE][nelem].RadRec_effec + 
			ionbal.CotaRate[nelem]) )*
			iso_sp[ipH_LIKE][nelem].xIonSimple);

		/* fraction of ionization due to collisional ionization */
		colfrc = (iso_sp[ipH_LIKE][nelem].fb[ipH1s].ColIoniz*dense.EdenHCorr )/
			((dense.eden*(iso_sp[ipH_LIKE][nelem].RadRec_effec + 
			ionbal.CotaRate[0]) )*
			iso_sp[ipH_LIKE][nelem].xIonSimple);

		/* fraction of ionization due to secondary ionization */
		secfrc = secondaries.csupra[nelem][nelem]/((dense.eden*(iso_sp[ipH_LIKE][nelem].RadRec_effec + 
			ionbal.CotaRate[0]) )*
			iso_sp[ipH_LIKE][nelem].xIonSimple);
	}
	else
	{
		phtfrc = 0.;
		colfrc = 0.;
		secfrc = 0.;
	}

	fprintf( ioQQQ, "     HydroLevel Z=%2ld called, simple II/I=",nelem);
	PrintE93( ioQQQ, iso_sp[ipH_LIKE][nelem].xIonSimple);
	fprintf( ioQQQ," PhotFrc:");
	PrintE93( ioQQQ,phtfrc);
	fprintf(ioQQQ," ColFrc:");
	PrintE93( ioQQQ,colfrc);
	fprintf( ioQQQ," SecFrc");
	PrintE93( ioQQQ, secfrc);
	fprintf( ioQQQ,"  Te:");
	PrintE93( ioQQQ,phycon.te);
	fprintf( ioQQQ," eden:");
	PrintE93( ioQQQ,dense.eden);
	fprintf( ioQQQ,"\n"); 
	return;
}

/*PrtHydroTrace1 print trace info for hydrogen-like species */
STATIC void PrtHydroTrace1(long nelem )
{
	long int ipHi , ipLo , i;

	DEBUG_ENTRY( "PrtHydroTrace1()" );

	fprintf( ioQQQ, 
		"       HydroLevel%3ld finds arrays, with optical depths defined? %li induced 2ph=%12.3e\n", 
		nelem, iteration, iso_sp[ipH_LIKE][nelem].TwoNu[0].induc_up );
	/* 06 aug 28, from numLevels_max to _local. */
	for( ipHi=ipH2s; ipHi < iso_sp[ipH_LIKE][nelem].numLevels_local; ipHi++ )
	{
		fprintf( ioQQQ, "up:%2ld", ipHi );
		fprintf( ioQQQ, "lo" );
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			fprintf( ioQQQ, "%9ld", ipLo );
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "%3ld", ipHi );
		fprintf( ioQQQ, " A*esc" );
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			fprintf( ioQQQ,PrintEfmt("%9.2e",  iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Aul()*
				iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Pesc() ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "%3ld", ipHi );
		fprintf( ioQQQ, " A*ees" );
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			fprintf( ioQQQ,PrintEfmt("%9.2e",  iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Aul()*
				iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Pelec_esc() ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "%3ld", ipHi );
		fprintf( ioQQQ, " tauin" );
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			fprintf( ioQQQ,PrintEfmt("%9.2e",  iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().TauIn() ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "%3ld", ipHi );
		fprintf( ioQQQ, " t tot" );
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().TauTot() ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "%3ld", ipHi );
		fprintf( ioQQQ, " Esc  " );
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Pesc() ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "%3ld", ipHi );
		fprintf( ioQQQ, " Eesc " );
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Pelec_esc() ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "%3ld", ipHi );
		fprintf( ioQQQ, " Dest " );
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			fprintf( ioQQQ,PrintEfmt("%9.2e",  iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Pdest()) );
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "%3ld", ipHi );
		fprintf( ioQQQ, " A*dst" );
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			fprintf( ioQQQ,PrintEfmt("%9.2e",  iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Aul()*
				iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Pdest() ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "%3ld", ipHi );
		fprintf( ioQQQ, " StrkE" );
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			fprintf( ioQQQ,PrintEfmt("%9.2e",  iso_sp[ipH_LIKE][nelem].ex[ipHi][ipLo].pestrk_up ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "%3ld", ipHi );
		fprintf( ioQQQ, " B(ul)" );
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().pump()*
				iso_sp[ipH_LIKE][nelem].st[ipLo].g()/iso_sp[ipH_LIKE][nelem].st[ipHi].g() ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "%3ld", ipHi );
		fprintf( ioQQQ, " tcont" );
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			fprintf( ioQQQ,PrintEfmt("%9.2e",  iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().TauCon() ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "%3ld", ipHi );
		fprintf( ioQQQ, " C(ul)" );
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Coll().ColUL( colliders ) ));
		}
		fprintf( ioQQQ, "\n" );

		if( ipHi == 2 )
		{
			fprintf( ioQQQ, "    FeIIo");
			fprintf( ioQQQ,PrintEfmt("%9.2e", 
				hydro.dstfe2lya* iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Aul() ));
			fprintf( ioQQQ, "\n");
		}
	}

	fprintf( ioQQQ, "         " );
	/* 06 aug 28, from numLevels_max to _local. */
	for( i=1; i < iso_sp[ipH_LIKE][nelem].numLevels_local; i++ )
	{
		fprintf( ioQQQ, "%9ld", i );
	}
	fprintf( ioQQQ, "\n" );
	return;
}

/*HydroLevel calls iso_level to solve for ionization balance 
 * and level populations of model hydrogen atom */
void HydroLevel(long int nelem)
{
	long int i; 
	int ipISO = ipH_LIKE;

	DEBUG_ENTRY( "HydroLevel()" );

	/* check that we were called with valid charge */
	ASSERT( nelem >= 0);
	ASSERT( nelem < LIMELM );

	/* option to print some rates */
	if( (trace.lgTrace && trace.lgIsoTraceFull[ipISO]) && (nelem == trace.ipIsoTrace[ipISO]) )
		PrtHydroTrace1(nelem);

	if( trace.lgHBug && trace.lgTrace )
		PrtHydroTrace1a(nelem);

	/* this is main trace h-like printout */
	if( (trace.lgIsoTraceFull[ipISO] && trace.lgTrace) && (nelem == trace.ipIsoTrace[ipISO]) )
	{
		fprintf( ioQQQ, "       HLEV HGAMNC" );
		PrintE93( ioQQQ, iso_sp[ipISO][nelem].fb[ipH1s].gamnc );
		/* 06 aug 28, from numLevels_max to _local. */
		for( i=ipH2s; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].fb[i].gamnc ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "       HLEV TOTCAP" );
		/* 06 aug 28, from numLevels_max to _local. */
		for( i=1; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].fb[i].RateCont2Level/dense.eden ));
		}
		fprintf( ioQQQ," tot");
		fprintf( ioQQQ,PrintEfmt("%10.2e", ionbal.RateRecomTot[nelem][nelem-ipISO]/dense.eden ) );
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "       HLEV IND Rc" );
		/* 06 aug 28, from numLevels_max to _local. */
		for( i=ipH1s; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].fb[i].RecomInducRate*iso_sp[ipISO][nelem].fb[i].PopLTE ));
		}
		fprintf( ioQQQ, "\n" );

		/* incuded recombination rate coefficients */
		fprintf( ioQQQ, "       IND Rc LTE " );
		/* 06 aug 28, from numLevels_max to _local. */
		for( i=ipH1s; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e",
				iso_sp[ipISO][nelem].fb[i].gamnc*iso_sp[ipISO][nelem].fb[i].PopLTE ));
		}
		fprintf( ioQQQ, "\n" );

		/* LTE level populations */
		fprintf( ioQQQ, "       HLEV   HLTE" );
		/* 06 aug 28, from numLevels_max to _local. */
		for( i=ipH1s; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].fb[i].PopLTE ));
		}
		fprintf( ioQQQ, "\n" );

		/* fraction of total ionization due to collisions from given level */
		fprintf( ioQQQ, "       HLEVfr cion" );
		/* 06 aug 28, from numLevels_max to _local. */
		for( i=ipH1s; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", 
				iso_sp[ipISO][nelem].fb[i].ColIoniz*dense.EdenHCorr/MAX2(1e-30,iso_sp[ipISO][nelem].fb[i].RateLevel2Cont) ) );
		}
		fprintf( ioQQQ, "\n" );

		/* fraction of total ionization due to photoionization from given level */
		if( ionbal.RateRecomTot[nelem][nelem]> 0. )
		{
			fprintf( ioQQQ, "       HLEVfrPhIon" );
			/* 06 aug 28, from numLevels_max to _local. */
			for( i=ipH1s; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
			{
				fprintf(ioQQQ,PrintEfmt("%9.2e", 
					iso_sp[ipISO][nelem].fb[i].gamnc/MAX2(1e-30,iso_sp[ipISO][nelem].fb[i].RateLevel2Cont) ) );
			}
			fprintf( ioQQQ, "\n" );
		}

		fprintf( ioQQQ, "       HLEV     HN" );
		/* 06 aug 28, from numLevels_max to _local. */
		for( i=ipH1s; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].st[i].Pop() ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "       HLEV   b(n)" );
		/* 06 aug 28, from numLevels_max to _local. */
		for( i=ipH1s; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].fb[i].DepartCoef ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "       HLEV   X12tot");
		fprintf(ioQQQ,PrintEfmt("%9.2e", secondaries.x12tot ));
		fprintf( ioQQQ," Grn dest:");
		fprintf(ioQQQ,PrintEfmt("%9.2e",
		  ionbal.RateIoniz[nelem][nelem][nelem+1] ));
		fprintf(ioQQQ, "\n"); 
	}

	if( trace.lgTrace )
	{
		/* iso.RecomTotal[nelem] is gross rec coef, computed here while filling in matrix
		 * elements, all physical processes included. 
		 * RadRec_effec is total effective radiative only */
		fprintf( ioQQQ, "       HydroLevel Z:%2ld return %s te=",
			nelem,
			iso_sp[ipISO][nelem].chTypeAtomUsed );
		PrintE93( ioQQQ,phycon.te);
		fprintf( ioQQQ," density=%.4e", dense.xIonDense[nelem][nelem-ipISO] );

		fprintf( ioQQQ," simple=%.4e",iso_sp[ipISO][nelem].xIonSimple);

		fprintf( ioQQQ," b1=%.2e",iso_sp[ipISO][nelem].fb[ipH1s].DepartCoef);

		fprintf( ioQQQ," ion rate=%.4e",ionbal.RateIonizTot(nelem,nelem-ipISO) );

		fprintf( ioQQQ," TotRec=%.4e",ionbal.RateRecomTot[nelem][nelem-ipISO]);

		fprintf( ioQQQ," RadRec=%.4e",iso_sp[ipISO][nelem].RadRec_effec);
		fprintf( ioQQQ, "\n");
	}
	return;
}