/* 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 */
/*HydroRecom photoionization, recombination, radiative rates for model hydrogen atom */
#include "cddefines.h"
#include "hydrogenic.h"
#include "iso.h"
#include "opacity.h"
#include "trace.h"
#include "phycon.h"
#include "ionbal.h"
#include "punch.h"
#include "elementnames.h"
#include "atmdat.h"
#include "rt.h"

void HydroRecom(
	/* nelem is on the C scale, 0 for H */
	long int nelem)
{
	long int n;
	int lgOK;
	double extra, 
	  SumCaseB ,
	  SumTopOff;
	/* will be used to save temperature where rec coef evaluated */
	static double TeUsed[LIMELM]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
		0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};

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

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

	/* evaluate recombination escape probability for all levels */
	for( n=ipH1s; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
	{
		/* this is escape probability */
		iso.RadRecomb[ipH_LIKE][nelem][n][ipRecEsc] = 
			RT_recom_effic(iso.ipIsoLevNIonCon[ipH_LIKE][nelem][n]);

		/* otsmin set to zero in zerologic, set to 1 with 
		 * the NO ON THE SPOT command */
		iso.RadRecomb[ipH_LIKE][nelem][n][ipRecEsc] = 
			MAX2(iso.RadRecomb[ipH_LIKE][nelem][n][ipRecEsc], opac.otsmin);
	}
#	if 0
	n = iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][3];
	if( nelem==0 ) fprintf(ioQQQ,"%i\thydrogrecom optical depths\t%.2e\t%.2e\t%.2e \n",
		nzone,
		opac.TauAbsGeo[0][n-1],
		opac.TauAbsGeo[1][n-1],
		opac.TauAbsGeo[1][n-1] - opac.TauAbsGeo[0][n-1]);
#	endif

	/* option for case b conditions, kill ground state recombination */
	if( opac.lgCaseB )
	{
		iso.RadRecomb[ipH_LIKE][nelem][ipH1s][ipRecEsc] = 1e-10;
	}

	/* always reevaluate net escape probability, which depends on background opacities */
	for( n=ipH1s; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
	{
		/* This is the one that multiplies the radiative recombination rate for level n
		 * following accounts for absorption of diffuse continuum
		 * by .not.hydrogen, includes modified OTS factor*/
		iso.RadRecomb[ipH_LIKE][nelem][n][ipRecNetEsc] = 
			MIN2(1.,iso.RadRecomb[ipH_LIKE][nelem][n][ipRecEsc]+
		  (1.-iso.RadRecomb[ipH_LIKE][nelem][n][ipRecEsc])*iso.ConOpacRatio[ipH_LIKE][nelem][n]);

		ASSERT( iso.RadRecomb[ipH_LIKE][nelem][n][ipRecNetEsc] >= 0. && 
			iso.RadRecomb[ipH_LIKE][nelem][n][ipRecNetEsc] <= 1. );
	}
	/*if( nelem==1 ) fprintf(ioQQQ," he ground %.2e %.2e %.2e \n",
		iso.RadRecomb[ipH_LIKE][nelem][0][ipRecNetEsc], iso.RadRecomb[ipH_LIKE][nelem][0][ipRecEsc],iso.ConOpacRatio[ipH_LIKE][nelem][0]);*/

	if( (trace.lgTrace && trace.lgIsoTraceFull[ipH_LIKE]) && (nelem == trace.ipIsoTrace[ipH_LIKE]) )
	{
		/* print continuum escape probability */
		fprintf( ioQQQ, "       HydroRecom recomb effic Z%li ",nelem );
		for( n=ipH1s; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
		{
			fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRecomb[ipH_LIKE][nelem][n][ipRecEsc] ));
		}
		fprintf( ioQQQ, "\n" );

		/* net recombination efficiency factor, including background opacity*/
		fprintf( ioQQQ, "       HydroRecom recomb net effic" );
		for( n=ipH1s; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
		{
			fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRecomb[ipH_LIKE][nelem][n][ipRecNetEsc]) );
		}
		fprintf( ioQQQ, "\n" );

		/* inward continuous optical depths */
		fprintf( ioQQQ, "       HydroRecom in optic dep%10.3e", 
		  opac.TauAbsGeo[0][iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s]-1] );
		for( n=2; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
		{
			fprintf( ioQQQ,PrintEfmt("%10.3e",  opac.TauAbsGeo[0][iso.ipIsoLevNIonCon[ipH_LIKE][nelem][n]-1] ));
		}
		fprintf( ioQQQ, "\n" );

		/* outward continuous optical depths*/
		fprintf( ioQQQ, "       HydroRecom out op depth%10.3e", 
		  opac.TauAbsGeo[1][iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s]-1] );
		for( n=2; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
		{
			fprintf( ioQQQ,PrintEfmt("%10.3e",  opac.TauAbsGeo[1][iso.ipIsoLevNIonCon[ipH_LIKE][nelem][n]-1] ));
		}
		fprintf( ioQQQ, "\n" );
	}

	if( (fabs(1.-TeUsed[nelem]/phycon.te)> 0.001) || hydro.lgReevalRecom )
	{
		TeUsed[nelem] = phycon.te;
		/* evaluate radiative recom coef */
		/* first the ground state with Dima's routine */
		iso.RadRecomb[ipH_LIKE][nelem][ipH1s][ipRecRad] = atmdat_H_rad_rec(nelem+1,ipH1s,phycon.te);

		/* evaluate excited state rec coef, also do sum over all excited states, 
		 * SumCaseB will be true case B sum */
		SumCaseB = 0.;
		for( n=ipH2s; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
		{
			iso.RadRecomb[ipH_LIKE][nelem][n][ipRecRad] = atmdat_H_rad_rec(nelem+1,n,phycon.te);
			SumCaseB += iso.RadRecomb[ipH_LIKE][nelem][n][ipRecRad];
		}

		/* following returns case A rec coef, then we will subtract ground */
		iso.RadRec_caseB[ipH_LIKE][nelem] = atmdat_H_rad_rec(nelem+1,-1,phycon.te);

		/* summed rad rec rate coef needed for recom cooling */
		ionbal.RR_rate_coef_used[nelem][nelem] = iso.RadRec_caseB[ipH_LIKE][nelem];

		/* save the total recombination rate, only used for printouts in recombmake,
		 * at this point CaseBradrec actually has the total recombination*/
		/*ionbal.RecombinRate[nelem][nelem] = iso.RadRec_caseB[ipH_LIKE][nelem]*(float)dense.eden; */

		/* subtract ground to get correct case b radiative sum */
		iso.RadRec_caseB[ipH_LIKE][nelem] -= iso.RadRecomb[ipH_LIKE][nelem][ipH1s][ipRecRad];

		/* extra is +/- part needed to be added on to get correct total */
		extra = iso.RadRec_caseB[ipH_LIKE][nelem] - SumCaseB;

		/* make sure positive number of topoff levels */
		ASSERT( (iso.numLevels_max[ipH_LIKE][nelem] > iso.nTopOff[ipH_LIKE][nelem]) || 
			iso.lgFullSize[ipH_LIKE][nelem] );

		/* only top off if this extra is positive since for large atom, sum is better
		 * than case b fit/
		 * adjust recombinations to levels to get total sum correct 
		 * default topoff is " add"
		 * with very large atom and high temperatures, extra can be negative,
		 * in which case we want to scale to make sure nothing is negative*/
		if( extra > 0. && !iso.lgFullSize[ipH_LIKE][nelem] )
		{
			if( strcmp(hydro.chHTopType,"scal") == 0 )
			{
				/* this is now scale factor to get correct sum, we will add for n=3 to end
				 * nHTopOff default is 10, change with atom h-like topoff */
				SumTopOff = 0;
				for( n=iso.nTopOff[ipH_LIKE][nelem]; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
				{
					SumTopOff += iso.RadRecomb[ipH_LIKE][nelem][n][ipRecRad];
				}

				/* scale factor to get right sum */
				if( SumTopOff > 1e-30 )
				{
					extra = (1. + extra/SumTopOff);
				}
				else
				{
					extra = 1.;
				}

				/* this is the top-off top off of the h atom */
				for( n=iso.nTopOff[ipH_LIKE][nelem]; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
				{
					iso.RadRecomb[ipH_LIKE][nelem][n][ipRecRad] *= extra;
				}
			}
			else if( strcmp(hydro.chHTopType," add") == 0 )
			{
				/* spread it out by adding same to all top levels
				 * this is default */
				extra /= (iso.numLevels_max[ipH_LIKE][nelem] - iso.nTopOff[ipH_LIKE][nelem]);
				for( n=iso.nTopOff[ipH_LIKE][nelem]; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
				{
					iso.RadRecomb[ipH_LIKE][nelem][n][ipRecRad] += extra;
				}
			}
			else
			{
				/* this is insanity */
				fprintf( ioQQQ, " HydroRecom has insane chHTopType=%4.4s\n", 
				  hydro.chHTopType );
				ShowMe();
				puts( "[Stop in hydrorecom]" );
				cdEXIT(EXIT_FAILURE);
			}
		}

		/*begin sanity check */
		lgOK = TRUE;
		for( n=ipH1s; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
		{
			if( iso.RadRecomb[ipH_LIKE][nelem][n][ipRecRad] <= 0. )
			{
				fprintf( ioQQQ, 
					" HydroRecom non-positive recombination coefficient for Z=%3ld lev n=%3ld rec=%11.2e te=%11.2e\n", 
				  nelem, n, iso.RadRecomb[ipH_LIKE][nelem][n][ipRecRad] , phycon.te);
					lgOK = FALSE;
			}
		}
		/* bail if we found problems */
		if( !lgOK )
		{
			ShowMe();
			puts( "[Stop in HydroRecom]" );
			cdEXIT(EXIT_FAILURE);
		}
		/*end sanity check */

		if( (trace.lgTrace && trace.lgIsoTraceFull[ipH_LIKE]) && (nelem == trace.ipIsoTrace[ipH_LIKE]) )
		{
			fprintf( ioQQQ, "       HydroRecom eval rec cof" );
			for( n=ipH1s; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRecomb[ipH_LIKE][nelem][n][ipRecRad]) );
			}
			fprintf( ioQQQ, "\n" );
		}
	} 
	/* end branch checking on change in temperature */

	/* confirm that we have good rec coef at bottom and top of h atom */
	ASSERT( iso.RadRecomb[ipH_LIKE][nelem][ipH1s][ipRecRad] > 0. );
	ASSERT( iso.RadRecomb[ipH_LIKE][nelem][iso.numLevels_max[ipH_LIKE][nelem]-1][ipRecRad] > 0. );

	/* RadRec_effec is total effective radiative recombination */
	iso.RadRec_effec[ipH_LIKE][nelem] = 0.;
	for( n=ipH1s; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
	{
		iso.RadRec_effec[ipH_LIKE][nelem] += iso.RadRecomb[ipH_LIKE][nelem][n][ipRecRad]*
		  iso.RadRecomb[ipH_LIKE][nelem][n][ipRecNetEsc];
	}

	/* set true when punch recombinatino coefficients command entered */
	if( punch.lgioRecom )
	{
		if( nelem==0 )
		{
			fprintf( punch.ioRecom, "Hydrogenic species\n");
		}
		/* this prints Z on physical, not C, scale */
		fprintf( punch.ioRecom, "H-like %2li %s ", 
			nelem+1 , elementnames.chElementSym[nelem] );
		fprintf( punch.ioRecom,PrintEfmt("%9.2e", iso.RadRec_caseB[ipH_LIKE][nelem] ));
		fprintf( punch.ioRecom, "\n" );
	}

	if( trace.lgTrace && trace.lgHBug )
	{
		fprintf( ioQQQ, "     HydroRecom Z=%3ld total rec coef", nelem );
		fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRec_effec[ipH_LIKE][nelem] ));
		fprintf( ioQQQ, " case A=" );
		fprintf( ioQQQ,PrintEfmt("%10.3e", 
			iso.RadRec_caseB[ipH_LIKE][nelem] + iso.RadRecomb[ipH_LIKE][nelem][ipH1s][ipRecRad] ) );
		fprintf( ioQQQ, " caseB=");
		fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRec_caseB[ipH_LIKE][nelem] ));
		fprintf( ioQQQ, "\n" );
	}

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



