/*HydroRecom photoionization, recombination, radiative rates for model hydrogen atom */
#include <string.h>
#include "cddefines.h"
#include "showme.h"
#include "hydrogenic.h"
#include "opac.h"
#include "nh.h"
#include "trace.h"
#include "otsmin.h"
#include "phycon.h"
#include "ophf.h"
#include "hcaseb.h"
#include "caseb.h"
#include "ioreco.h"
#include "htopoff.h"
#include "elementnames.h"
#include "printe82.h"
#include "hrfit.h"
#include "receff.h"

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

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

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

	/* evaluate recombination escape probability for all levels */
	for( n=IP1S; n <= nhlevel; n++ )
	{
		/* this is escape probability */
		hydro.hrec[ipZ][n][ipRecEsc] = (float)receff(nh.ipHn[ipZ][n]);

		/* otsmin set to zero in zerologic, set to 1 with 
		 * the NO ON THE SPOT command */
		hydro.hrec[ipZ][n][ipRecEsc] = 
			(float)MAX2(hydro.hrec[ipZ][n][ipRecEsc], otsminCom.otsmin);
	}

	/* option for case b conditions, kill ground state recombination */
	if( caseb.lgCaseb )
	{
		hydro.hrec[ipZ][IP1S][ipRecEsc] = 1e-10f;
	}

	/* always reevaluate net escape probability, which depends on background opacities */
	for( n=IP1S; n <= nhlevel; 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*/
		hydro.hrec[ipZ][n][ipRecNetEsc] = 
			(float)MIN2(1.,hydro.hrec[ipZ][n][ipRecEsc]+
		  (1.-hydro.hrec[ipZ][n][ipRecEsc])*ophf.HOpacRatio[ipZ][n]);
	}

	if( (trace.lgTrace && trace.lgHBugFull) && (ipZ == trace.ipZTrace) )
	{
		/* print continuum escape probability */
		fprintf( ioQQQ, "       HydroRecom recomb effic Z%li ",ipZ );
		for( n=IP1S; n <= nhlevel; n++ )
		{
			fprintf( ioQQQ,PrintEfmt("%10.3e", hydro.hrec[ipZ][n][ipRecEsc] ));
		}
		fprintf( ioQQQ, "\n" );

		/* net recombination efficiency factor, including background opacity*/
		fprintf( ioQQQ, "       HydroRecom recomb net effic" );
		for( n=IP1S; n <= nhlevel; n++ )
		{
			fprintf( ioQQQ,PrintEfmt("%10.3e", hydro.hrec[ipZ][n][ipRecNetEsc]) );
		}
		fprintf( ioQQQ, "\n" );

		/* inward continuous optical depths */
		fprintf( ioQQQ, "       HydroRecom in optic dep%10.3e", 
		  opac.tauabs[0][nh.ipHn[ipZ][IP1S]-1] );
		for( n=2; n <= nhlevel; n++ )
		{
			fprintf( ioQQQ,PrintEfmt("%10.3e",  opac.tauabs[0][nh.ipHn[ipZ][n]-1] ));
		}
		fprintf( ioQQQ, "\n" );

		/* outward continuous optical depths*/
		fprintf( ioQQQ, "       HydroRecom out op depth%10.3e", 
		  opac.tauabs[1][nh.ipHn[ipZ][IP1S]-1] );
		for( n=2; n <= nhlevel; n++ )
		{
			fprintf( ioQQQ,PrintEfmt("%10.3e",  opac.tauabs[1][nh.ipHn[ipZ][n]-1] ));
		}
		fprintf( ioQQQ, "\n" );
	}

	if( fabs(1.-TeUsed[ipZ]/phycon.te)> 0.001 )
	{
		TeUsed[ipZ] = phycon.te;
		/* evaluate radiative recom coef */
		/* first the ground state with Dima's routine */
		hrfit(ipZ+1,IP1S,phycon.te,&hydro.hrec[ipZ][IP1S][ipRecRad]);

		/* evaluate excited state rec coef, also do sum over all excited states, 
		 * SumCaseB will be true case B sum */
		SumCaseB = 0.;
		for( n=IP2S; n <= nhlevel; n++ )
		{
			hrfit(ipZ+1,n,phycon.te,&hydro.hrec[ipZ][n][ipRecRad]);
			SumCaseB += hydro.hrec[ipZ][n][ipRecRad];
		}

		/* following returns case A rec coef and saves into HRecRad */
		hrfit(ipZ+1,-1,phycon.te,&hcaseb.HRecRad[ipZ]);

		/* subtract ground to get correct case b radiative sum */
		hcaseb.HRecRad[ipZ] -= hydro.hrec[ipZ][IP1S][ipRecRad];

		/* extra is +/- part needed to be added on to get correct total */
		extra = hcaseb.HRecRad[ipZ] - SumCaseB;

		/* 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. )
		{
			if( strcmp(HTopOff.chHTopType,"scal") == 0 )
			{
				/* this is now scale factor to get correct sum, we will add for n=3 to end
				 * nHTopOff default is 5, change with set nHTopOff */
				SumTopOff = 0;
				for( n=HTopOff.nHTopOff; n <= nhlevel; n++ )
				{
					SumTopOff += hydro.hrec[ipZ][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=HTopOff.nHTopOff; n <= nhlevel; n++ )
				{
					hydro.hrec[ipZ][n][ipRecRad] *= (float)extra;
				}
			}
			else if( strcmp(HTopOff.chHTopType," add") == 0 )
			{
				/* spread it out by adding same to all top levels
				 * this is default */
				extra /= nhlevel - HTopOff.nHTopOff + 1;
				for( n=HTopOff.nHTopOff; n <= nhlevel; n++ )
				{
					hydro.hrec[ipZ][n][ipRecRad] += (float)extra;
				}
			}
			else
			{
				/* this is insanity */
				fprintf( ioQQQ, " HydroRecom has insane chHTopType=%4.4s\n", 
				  HTopOff.chHTopType );
				ShowMe();
				puts( "[Stop in hydrorecom]" );
				exit(1);
			}
		}

		/*begin sanity check */
		lgOK = TRUE;
		for( n=IP1S; n <= nhlevel; n++ )
		{
			if( hydro.hrec[ipZ][n][ipRecRad] <= 0. )
			{
				fprintf( ioQQQ, 
					" HydroRecom non-positive recombination coefficient for Z=%3ld lev n=%3ld rec=%11.2e te=%11.2e\n", 
				  ipZ, n, hydro.hrec[ipZ][n][ipRecRad] , phycon.te);
					lgOK = FALSE;
			}
		}
		/* bail if we found problems */
		if( !lgOK )
		{
			ShowMe();
			puts( "[Stop in HydroRecom]" );
			exit(1);
		}
		/*end sanity check */

		if( (trace.lgTrace && trace.lgHBugFull) && (ipZ == trace.ipZTrace) )
		{
			fprintf( ioQQQ, "       HydroRecom eval rec cof" );
			for( n=IP1S; n <= nhlevel; n++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e", hydro.hrec[ipZ][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( hydro.hrec[ipZ][IP1S][ipRecRad] > 0. );
	assert( hydro.hrec[ipZ][nhlevel][ipRecRad] > 0. );

	/* htotrc(ipZ) is total effective radiative recombination */
	hcaseb.htotrc[ipZ] = 0.;
	for( n=IP1S; n <= nhlevel; n++ )
	{
		hcaseb.htotrc[ipZ] += hydro.hrec[ipZ][n][ipRecRad]*
		  hydro.hrec[ipZ][n][ipRecNetEsc];
	}

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

	if( trace.lgTrace && trace.lgHBug )
	{
		fprintf( ioQQQ, "       HydroRecom Z=%3ld total rec coef", ipZ );
		fprintf( ioQQQ,PrintEfmt("%10.3e", hcaseb.htotrc[ipZ] ));
		fprintf( ioQQQ, " case A=" );
		fprintf( ioQQQ,PrintEfmt("%10.3e", 
			hcaseb.HRecRad[ipZ] + hydro.hrec[ipZ][IP1S][ipRecRad] ) );
		fprintf( ioQQQ, " caseB=");
		fprintf( ioQQQ,PrintEfmt("%10.3e", hcaseb.HRecRad[ipZ] ));
		fprintf( ioQQQ, "\n" );
	}

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