/* 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 */
#include "cddefines.h"
#include "physconst.h"
#include "math_complex.h"
#include "gaunt.h"
#include "thirdparty.h"

static double RealF2_1( double alpha, double beta, double gamma, double chi );
static complex Hypergeometric2F1( complex a, complex b, complex c,
		double chi, long *NumRenorms, long *NumTerms );
static complex F2_1( complex alpha, complex beta, complex gamma,
		double chi, long *NumRenormalizations, long *NumTerms );
static complex HyperGeoInt( double v );
static complex qg32complex(	double xl, double xu, complex (*fct)(double) );
static double GauntIntegrand( double y );
static double FreeFreeGaunt( double x );
static double DoBeckert_etal( double etai, double etaf, double chi );
static double DoSutherland( double etai, double etaf, double chi );

/* used to keep intermediate results from over- or underflowing.	*/
static complex Normalization = {1e100, 1e100};
static complex CMinusBMinus1, BMinus1, MinusA;
static double GlobalCHI;
static double Zglobal, HNUglobal, TEglobal;

double CalcThermAveGaunt( double temp, double z, double photon )
{
	double gaunt, u, gamma2;
	
	Zglobal = z;
	HNUglobal = photon;
	TEglobal = temp;

	u = TE1RYD*photon/temp;
	gamma2=TE1RYD*z*z/temp;

	if( log10(u)<-5. )
	{
		if( log10( gamma2 ) < 0. )
		{
			/* Elwert(1954), given as eqn 3.2 in Hummer 88	*/
			gaunt = 0.551329 * ( 0.80888 - log(u) );
		}
		else
		{
			gaunt = -0.551329 * (0.5*log(gamma2) + log(u) + 0.056745);
		}
	}
	else
	{
		/* Perform integration.	*/
		gaunt =  qg32( 0.01, 1., GauntIntegrand );
		gaunt += qg32( 1., 5., GauntIntegrand );
	}
	ASSERT( gaunt>0. && gaunt<100. );

	return gaunt;
}

static double GauntIntegrand( double y )
{
	double value;
	value =	FreeFreeGaunt( y ) * exp(-y);
	return value;
}

static double FreeFreeGaunt( double x )
{
	double Csum, zeta, etaf, etai, chi, gaunt, z, InitialElectronEnergy, FinalElectronEnergy, photon;
	int lgSutherlandOn = FALSE;
	long i;
	
	z = Zglobal;
	photon = HNUglobal;
	ASSERT( z > 0. );
	ASSERT( photon > 0. );

	/* The final electron energy should be xkT + hv and the initial, just xkT	*/
	InitialElectronEnergy = sqrt(x) * TEglobal/TE1RYD;
	FinalElectronEnergy = photon + InitialElectronEnergy;
	ASSERT( InitialElectronEnergy > 0. );

	/* These are the free electron analog to a bound state principal quantum number.*/
	etai = z/sqrt(InitialElectronEnergy);	
	etaf = z/sqrt(FinalElectronEnergy);
	ASSERT( etai > 0. );
	ASSERT( etaf > 0. );
	chi = -4. * etai * etaf / POW2( etai - etaf );
	zeta = etai-etaf;

	if( etai>=130.) 
	{
		if( etaf < 1.7 )
		{
			/* Brussard and van de Hulst (1962), as given in Hummer 88, eqn 2.23b	*/
			gaunt = 1.1027 * (1.-exp(-2.*PI*etaf));
		}
		else if( etaf < 0.1*etai )
		{
			/* Hummer 88, eqn 2.23a	*/
			gaunt = 1. + 0.17282604*pow(etaf,-0.67) - 0.04959570*pow(etaf,-1.33) 
				- 0.01714286*pow(etaf,-2.) + 0.00204498*pow(etaf,-2.67) 
				- 0.00243945*pow(etaf,-3.33) - 0.00120387*pow(etaf,-4.) 
				+ 0.00071814*pow(etaf,-4.67) + 0.00026971*pow(etaf,-5.33);
		}
		else if( zeta > 0.5 )
		{
			/* Grant 1958, as given in Hummer 88, eqn 2.25a	*/
			gaunt = 1. + 0.21775*pow(zeta,-0.67) - 0.01312*pow(zeta, -1.33);
		}
		else 
		{
			double a[10] = {1.47864486, -1.72329012, 0.14420320, 0.05744888, 0.01668957,
				0.01580779,  0.00464268, 0.00385156, 0.00116196, 0.00101967};

			Csum = 0.;
			for( i = 0; i <=9; i++ )
			{
				/* The Chebyshev of the first kind is just a special kind of hypergeometric.	*/
				Csum += a[i]*RealF2_1( (double)(-i), (double)i, 0.5, 0.5*(1.-zeta) );
			}
			gaunt = fabs(0.551329*(0.57721 + log(zeta/2.))*exp(PI*zeta)*Csum);
			ASSERT( gaunt < 10. );
		}
	}
	else if( lgSutherlandOn )
		gaunt = DoSutherland( etai, etaf, chi ); 		
	else
		gaunt = DoBeckert_etal( etai, etaf, chi );

	/*if( gaunt*exp(-x) > 2. && TEglobal < 1000. )
		fprintf( ioQQQ,"ni %.3e nf %.3e chi %.3e u %.3e gam2 %.3e gaunt %.3e x %.3e\n", 
				etai, etaf, chi, TE1RYD*HNUglobal/TEglobal,
				TE1RYD*z*z/TEglobal, gaunt, x);*/
	
	/* TODO	2	- These are liberal bounds, in final product, this
	 * ASSERT should be much more demanding.	*/
	ASSERT( gaunt > 0. && gaunt<BIGFLOAT );
	
	if( gaunt == 0. )
	{
		fprintf( ioQQQ, "Uh-Oh! Gaunt is zero!  Is this okay?\n");
		/* assign some small value */
		gaunt = 1e-5;
	}

	return gaunt;
}


/************************************
 * This part calculates the gaunt factor as per Beckert et al 2000	*/
/* TODO	2	- insert reference	*/
static double DoBeckert_etal( double etai, double etaf, double chi )
{
	double Delta, BeckertGaunt, MaxFReal, LnBeckertGaunt;
	long NumRenorms[2]={0,0}, NumTerms[2]={0,0};
	int IndexMinNumRenorms, IndexMaxNumRenorms;
	complex a,b,c,F[2];
	
	a = ftocm( 1., -etai );
	b = ftocm( 0., -etaf );
	c = ftocm( 1., 0. );
	
	/* evaluate first hypergeometric function.	*/	
	F[0] = Hypergeometric2F1( a, b, c, chi, &NumRenorms[0], &NumTerms[0] );
	
	a.im = -etaf;
	b.im = -etai;
	
	/* evaluate second hypergeometric function.	*/	
	F[1] = Hypergeometric2F1( a, b, c, chi, &NumRenorms[1], &NumTerms[1] );

	/* If there is a significant difference in the number of terms used, 
	 * they should be recalculated with the max	number of terms in initial calculations */
	/* If NumTerms[i]=-1, the hypergeometric was calculated by the use of an integral instead
	 * of series summation...hence NumTerms has no meaning, and no need to recalculate.	*/
	if( ( MAX2(NumTerms[1],NumTerms[0]) - MIN2(NumTerms[1],NumTerms[0]) >= 2  )
		&& NumTerms[1]!=-1 && NumTerms[0]!=-1)
	{
		a = ftocm( 1., -etai );
		b = ftocm( 0., -etaf );
		c = ftocm( 1., 0. );

		NumTerms[0] = MAX2(NumTerms[1],NumTerms[0])+1;
		NumTerms[1] = NumTerms[0];
		NumRenorms[0] = 0;
		NumRenorms[1] = 0;
		
		/* evaluate first hypergeometric function.	*/	
		F[0] = Hypergeometric2F1( a, b, c, chi, &NumRenorms[0], &NumTerms[0] );
		
		a.im = -etaf;
		b.im = -etai;
		
		/* evaluate second hypergeometric function.	*/	
		F[1] = Hypergeometric2F1( a, b, c, chi, &NumRenorms[1], &NumTerms[1] );

		ASSERT( NumTerms[0] == NumTerms[1] );
	}

	/* if magnitude of unNormalized F's are vastly different, zero out the lesser	*/
	if( log10(cmabs(F[0])/cmabs(F[1])) + (NumRenorms[0]-NumRenorms[1])*log10(cmabs(Normalization)) > 10. )
	{
		F[1].im = 0.;
		F[1].re = 0.;
		/*  no longer need to keep track of differences in NumRenorms	*/
		NumRenorms[1] = NumRenorms[0];
	}
	else if( log10(cmabs(F[1])/cmabs(F[0])) + (NumRenorms[1]-NumRenorms[0])*log10(cmabs(Normalization)) > 10. )
	{
		F[0].im = 0.;
		F[0].re = 0.;
		/*  no longer need to keep track of differences in NumRenorms	*/
		NumRenorms[0] = NumRenorms[1];
	}
	
	/* now must fix if NumRenorms[0] != NumRenorms[1], because the next calculation is the
	 * difference of squares...information is lost and cannot be recovered if this calculation
	 * is done with NumRenorms[0] != NumRenorms[1]	*/
	MaxFReal = (fabs(F[1].re)>fabs(F[0].re)) ? fabs(F[1].re):fabs(F[0].re);
	while( NumRenorms[0] != NumRenorms[1] )
	{
		/* but must be very careful to prevent both overflow and underflow.	*/
		if( MaxFReal > 1e50 )
		{
			IndexMinNumRenorms = ( NumRenorms[0] > NumRenorms[1] ) ? 1:0;
			F[IndexMinNumRenorms] = cmdiv( F[IndexMinNumRenorms], Normalization );
			++NumRenorms[IndexMinNumRenorms];
		}
		else
		{
			IndexMaxNumRenorms = ( NumRenorms[0] > NumRenorms[1] ) ? 0:1;
			F[IndexMaxNumRenorms] = cmmul( F[IndexMaxNumRenorms], Normalization );
			--NumRenorms[IndexMaxNumRenorms];
		}
	}

	ASSERT( NumRenorms[0] == NumRenorms[1] );

	/* Okay, now we are guaranteed (?) a small dynamic range, but may still have to renormalize	*/

	/* Are we gonna have an overflow or underflow problem?	*/
	ASSERT( (fabs(F[0].re)<1e+150) && (fabs(F[1].re)<1e+150) && (fabs(F[0].im)<1e+150) && (fabs(F[1].re)<1e+150) );
	ASSERT( (fabs(F[0].re)>1e-150) && ((fabs(F[0].im)>1e-150) || (cmabs(F[0])==0.)) );
	ASSERT( (fabs(F[1].re)>1e-150) && ((fabs(F[1].re)>1e-150) || (cmabs(F[1])==0.)) );

	Delta = cmabs( cmsub( cmmul( F[0], F[0] ), cmmul(  F[1], F[1] ) ) );

	ASSERT( Delta > 0. );

	/* Now multiply by the coefficient in Beckert 2000, eqn 7	*/
	if( etaf > 100. )
	{
		/* must compute logarithmically	if etaf too big for linear computation.	*/
		LnBeckertGaunt = 1.6940360 + log(Delta) + log(etaf) + log(etai) - log(fabs(etai-etaf)) - 6.2831853*etaf;
		LnBeckertGaunt += 2. * NumRenorms[0] * log(cmabs(Normalization));
		BeckertGaunt = exp( LnBeckertGaunt );
		NumRenorms[0] = 0;
	}
	else
	{
		BeckertGaunt = Delta*5.4413981*etaf*etai/fabs(etai - etaf)
			/(1.-exp(-6.2831853*etai) )/( exp(6.2831853*etaf) - 1.);
				
		while( NumRenorms[0] > 0 )
		{
			BeckertGaunt *= cmabs(Normalization);
			BeckertGaunt *= cmabs(Normalization);
			ASSERT( BeckertGaunt < BIGDOUBLE );
			--NumRenorms[0];
		} 
		ASSERT( NumRenorms[0] == 0 );
	}

	ASSERT( NumRenorms[0] == 0 );
	
	/*fprintf( ioQQQ,"etai %.3e etaf %.3e u %.3e B %.3e \n", 
		etai, etaf, TE1RYD * HNUglobal / TEglobal, BeckertGaunt );	*/

	return BeckertGaunt;
}


/************************************
 * This part calculates the gaunt factor as per Sutherland 98	*/
/* TODO	2	- insert reference	*/
static double DoSutherland( double etai, double etaf, double chi )
{
	double Sgaunt, ICoef, weightI1, weightI0;
	long i, NumRenorms[2]={0,0}, NumTerms[2]={0,0};
	complex a,b,c,GCoef,kfac,etasum,G[2],I[2],ComplexFactors,GammaProduct;

	kfac = ftocm( fabs((etaf-etai)/(etaf+etai)), 0. );
	etasum = ftocm( 0., etai + etaf );
	
	GCoef = cmpow(kfac, etasum);
	/* GCoef is a complex vector that should be contained within the unit circle.
	 * and have a non-zero magnitude.  Or is it ON the unit circle?	*/
	ASSERT( fabs(GCoef.re)<1.0 && fabs(GCoef.im)<1.0 && ( GCoef.re!=0. || GCoef.im!=0. ) );

	for( i = 0; i <= 1; i++ )
	{
		a = ftocm( i + 1., -etaf );
		b = ftocm( i + 1., -etai );
		c = ftocm( 2.*i + 2., 0. );
			
		/* First evaluate hypergeometric function.	*/	
		G[i] = Hypergeometric2F1( a, b, c, chi, &NumRenorms[i], &NumTerms[i] );
	}
			
	/* If there is a significant difference in the number of terms used, 
	 * they should be recalculated with the max	number of terms in initial calculations */
	/* If NumTerms[i]=-1, the hypergeometric was calculated by the use of an integral instead
	 * of series summation...hence NumTerms has no meaning, and no need to recalculate.	*/
	if( MAX2(NumTerms[1],NumTerms[0]) - MIN2(NumTerms[1],NumTerms[0]) > 2 
		&& NumTerms[1]!=-1 && NumTerms[0]!=-1  )
	{
		NumTerms[0] = MAX2(NumTerms[1],NumTerms[0]);
		NumTerms[1] = NumTerms[0];
		NumRenorms[0] = 0;
		NumRenorms[1] = 0;
		
		for( i = 0; i <= 1; i++ )
		{
			a = ftocm( i + 1., -etaf );
			b = ftocm( i + 1., -etai );
			c = ftocm( 2.*i + 2., 0. );
				
			G[i] = Hypergeometric2F1( a, b, c, chi, &NumRenorms[i], &NumTerms[i] );
		}

		ASSERT( NumTerms[0] == NumTerms[1] );
	}

	for( i = 0; i <= 1; i++ )
	{
		/* TODO	2	- this check may also too liberal.  */
		ASSERT( fabs(G[i].re)>0. && fabs(G[i].re)<1e100 && fabs(G[i].im)>0. && fabs(G[i].im)<1e100 );

		/* Now multiply by the coefficient in Sutherland 98, eqn 9	*/
		G[i] = cmmul( G[i], GCoef );

		/* This is the coefficient in equation 8 in Sutherland	*/
		/* Karzas and Latter give gammafun(2.*i+2.), Sutherland gives gammafun(2.*i+1.)	*/
		ICoef = 0.25*pow(-chi, (double)i+1.)*exp( 1.5708*fabs(etai-etaf) )/gammafun(2.*i+1.);
		GammaProduct = cmmul( cdgamma(ftocm(i+1.,etai)) , cdgamma(ftocm(i+1.,etaf)) );
		ICoef *= cmabs(GammaProduct);
			
		ASSERT( ICoef > 0. );

		I[i] = cmmul( ftocm(ICoef,0.), G[i] );

		while( NumRenorms[i] > 0 )
		{
			I[i] = cmmul( Normalization, I[i] );
			ASSERT( fabs(I[i].re) < BIGDOUBLE && fabs(I[i].im) < BIGDOUBLE );
			--NumRenorms[i];
		} 
	
		ASSERT( NumRenorms[i] == 0 );
	}

	weightI0 = POW2(etaf+etai);
	weightI1 = 2.*etaf*etai*sqrt(1. + etai*etai)*sqrt(1. + etaf*etaf);

	ComplexFactors = cmmul( I[0], cmsub( cmmul( ftocm(weightI0,0.), I[0] ) , cmmul( ftocm(weightI1,0.), I[1] ) ) );
	
	/* This is Sutherland equation 13	*/
	Sgaunt  = 1.10266 / etai / etaf * cmabs( ComplexFactors );

	return Sgaunt;
}

/* This routine is a wrapper for F2_1	*/
static complex Hypergeometric2F1( complex a, complex b, complex c,
						  double chi, long *NumRenorms, long *NumTerms )
{
	complex a1, b1, c1, a2, b2, c2, Result, Part[2], F[2];
	complex chifac, GammaProduct, Coef, FIntegral;
	/* TODO	2	- pick these interface values and stick with it...best results have been 0.4, 1.5	*/
	double Interface1 = 0.4, Interface2 = 10.;
	long N_Renorms[2], N_Terms[2], IndexMaxNumRenorms, lgDoIntegral = FALSE;
	
	N_Renorms[0] = *NumRenorms;
	N_Renorms[1] = *NumRenorms;
	N_Terms[0] = *NumTerms;
	N_Terms[1] = *NumTerms;
	
	/* positive and zero chi are not possible.	*/
	ASSERT( chi < 0. );

	/* We want to be careful about evaluating the hypergeometric 
	 * in the vicinity of chi=1.  So we employ three different methods...*/

	/* for small chi, we pass the parameters to the hypergeometric function as is.	*/
	if( fabs(chi) < Interface1 )
	{
		Result = F2_1( a, b, c, chi, &*NumRenorms, &*NumTerms );
	}
	/* for large chi, we use a relation given as eqn 5 in Nicholas 89.	*/
	else if( fabs(chi) > Interface2 )
	{
		/* b1 = 1.-c+a; */
		/* c1 = 1.-b+a;	*/
		a1 = a;
		b1 = cmadd( ftocm(1.,0.), cmsub(a,c) );
		c1 = cmadd( ftocm(1.,0.), cmsub(a,b) );
	
		/* b2 = 1.-c+b;	*/
		/* c2 = 1.-a+b;	*/
		a2 = b;
		b2 = cmadd( ftocm(1.,0.), cmsub(b,c) );
		c2 = cmadd( ftocm(1.,0.), cmsub(b,a) );
		
		chifac = ftocm( -chi, 0. );

		F[0] = F2_1(a1,b1,c1,1./chi,&N_Renorms[0], &N_Terms[0]);
		F[1] = F2_1(a2,b2,c2,1./chi,&N_Renorms[1], &N_Terms[1]);

		/* do it again if significant difference in number of terms.	*/
		if( MAX2(N_Terms[1],N_Terms[0]) - MIN2(N_Terms[1],N_Terms[0]) >= 2 )
		{
			N_Terms[0] = MAX2(N_Terms[1],N_Terms[0]);
			N_Terms[1] = N_Terms[0];
			N_Renorms[0] = *NumRenorms;
			N_Renorms[1] = *NumRenorms;
			
			F[0] = F2_1(a1,b1,c1,1./chi,&N_Renorms[0], &N_Terms[0]);
			F[1] = F2_1(a2,b2,c2,1./chi,&N_Renorms[1], &N_Terms[1]);
			ASSERT( N_Terms[0] == N_Terms[1] );
		}

		*NumTerms = MAX2(N_Terms[1],N_Terms[0]);

		/************************************************************************/
		/* Do the first part	*/
		GammaProduct = cdgamma( cmsub(b,a) );
		GammaProduct = cmdiv( GammaProduct, cdgamma( b ) );
		GammaProduct = cmdiv( GammaProduct, cdgamma( cmsub(c,a) ) );
		GammaProduct = cmmul( GammaProduct, cdgamma( c ) );
		
		/* divide the hypergeometric by (-chi)^a and multiply by GammaProduct	*/
		Part[0] = cmmul( cmdiv( F[0], cmpow(chifac,a) ), GammaProduct );

		/************************************************************************/
		/* Do the second part	*/
		GammaProduct = cdgamma( cmsub(a,b) );
		GammaProduct = cmdiv( GammaProduct, cdgamma( a ) );
		GammaProduct = cmdiv( GammaProduct, cdgamma( cmsub(c,b) ) );
		GammaProduct = cmmul( GammaProduct, cdgamma( c ) );
	
		/* divide the hypergeometric by (-chi)^b and multiply by GammaProduct	*/
		Part[1] = cmmul( cmdiv( F[1], cmpow(chifac,b) ), GammaProduct );
		
		/************************************************************************/
		/* Add the two parts to get the result.	*/

		/* First must fix it if N_Renorms[0] != N_Renorms[1]	*/
		if( N_Renorms[0] != N_Renorms[1] )
		{
			IndexMaxNumRenorms = ( N_Renorms[0] > N_Renorms[1] ) ? 0:1;
			Part[IndexMaxNumRenorms] = cmmul( Part[IndexMaxNumRenorms], Normalization );
			--N_Renorms[IndexMaxNumRenorms];
			/* Only allow at most a difference of one in number of renormalizations...
			 * otherwise something is really screwed up	*/
			ASSERT( N_Renorms[0] == N_Renorms[1] );
		}

		*NumRenorms = N_Renorms[0];

		Result = cmadd( Part[0], Part[1] );
	}
	/* And for chi of order 1, we use Nicholas 89, eqn 27.	*/
	else
	{
		/* the hypergeometric integral does not seem to work well.	*/
		if( lgDoIntegral /* && fabs(chi+1.)>0.1 */)
		{
			/* a and b are always interchangeable, assign the lesser to b to 
			 * prevent Coef from blowing up	*/
			if( cmabs(b) > cmabs(a) )
			{
				complex btemp = b;
				b = a;
				a = btemp;
			}
			Coef = cmdiv( cmdiv( cdgamma(c),cdgamma(b) ), cdgamma(cmsub(c,b)) );
			CMinusBMinus1 = cmsub( cmsub(c,b), ftocm(1.,0.) );
			BMinus1 = cmsub( b, ftocm(1.,0.) );
			MinusA = cmmul( a, ftocm(-1.,0.) );
			GlobalCHI = chi;
			FIntegral = qg32complex( 0., 0.5, HyperGeoInt );
			FIntegral = cmadd( FIntegral, qg32complex( 0.5, 1., HyperGeoInt ) );

			Result = cmmul( Coef, FIntegral );
			*NumTerms = -1;
			*NumRenorms = 0;
		}
		else
		{
			/*	Near chi=1 solution	*/
			a1 = a;
			b1 = cmsub( c, b );
			c1 = c;
			chifac = ftocm( 1.-chi, 0. );
			
			Result = cmdiv( F2_1(a1,b1,c1,chi/(chi-1.),&*NumRenorms,&*NumTerms), cmpow(chifac,a) );
		}
	}

	/* Limit the size of the returned value	*/
	while( fabs(Result.re) >= 1e50 )
	{
		Result = cmdiv( Result, Normalization);
		++*NumRenorms;
	}

	return Result;
}

/* This routine calculates hypergeometric functions */
static complex F2_1( complex alpha, complex beta, complex gamma, double chi, long *NumRenormalizations, long *NumTerms )
{
	long  i = 3, MinTerms;
	int lgNotConverged = TRUE;
	complex LastTerm, Term, Sum;

	MinTerms = MAX2( 3, *NumTerms );

	/* This is the first term of the hypergeometric series.	*/
	Sum = cmdiv( ftocm(1.,0.), Normalization );
	++*NumRenormalizations;
	
	/* This is the second term	*/
	LastTerm = cmmul( cmdiv(cmmul(cmmul(Sum,alpha),beta),gamma), ftocm(chi,0.) );

	Sum = cmadd( Sum, LastTerm );
		
	/* Every successive term is easily found by multiplying the last term
	 * by (alpha + i - 2)*(beta + i - 2)*chi/(gamma + i - 2)/(i-1.)	*/
	do{
		alpha.re++;
		beta.re++;
		gamma.re++;

		/* multiply old term by incremented alpha++*beta++/gamma++.  Also multiply by chi/(i-1.)	*/
		Term = cmmul( cmdiv(cmmul(cmmul(LastTerm,alpha),beta),gamma), ftocm(chi/(i-1.),0.) ); 

		Sum = cmadd( Sum, Term );
		
		/* Renormalize if too big	*/
		if( Sum.re > 1e100 )
		{
			Sum = cmdiv( Sum, Normalization );
			LastTerm = cmdiv( Term, Normalization );
			++*NumRenormalizations;
			/* notify of renormalization, and print the number of the term	*/
			fprintf( ioQQQ,"Hypergeometric: Renormalized at term %li.  Sum = %.3e %.3e\n",
				i, Sum.re, Sum.im);
		}
		else
			LastTerm = Term;

		/* Declare converged if this term does not affect Sum by much	*/
		/* Must do this with abs because terms alternate sign.	*/
		if( fabs(LastTerm.re/Sum.re)<0.001 && fabs(LastTerm.im/Sum.im)<0.001 )
			lgNotConverged = FALSE;

		if( *NumRenormalizations >= 5 )
		{
			fprintf( ioQQQ, "We've got too many (%li) renorms!\n",*NumRenormalizations );
		}
		
		++i;

	}while ( lgNotConverged || i<MinTerms );

	*NumTerms = i;

	return Sum;
}

/* This routine calculates hypergeometric functions */
static double RealF2_1( double alpha, double beta, double gamma, double chi )
{
	long  i = 3;
	int lgNotConverged = TRUE;
	double LastTerm, Sum;

	/* This is the first term of the hypergeometric series.	*/
	Sum = 1.;
	
	/* This is the second term	*/
	LastTerm = alpha*beta*chi/gamma;

	Sum += LastTerm;
		
	/* Every successive term is easily found by multiplying the last term
	 * by (alpha + i - 2)*(beta + i - 2)*chi/(gamma + i - 2)/(i-1.)	*/
	do{
		alpha++;
		beta++;
		gamma++;

		/* multiply old term by incremented alpha++*beta++/gamma++.  Also multiply by chi/(i-1.)	*/
		LastTerm *= alpha*beta*chi/gamma/(i-1.); 

		Sum += LastTerm;
		
		/* Declare converged if this term does not affect Sum by much	*/
		/* Must do this with abs because terms alternate sign.	*/
		if( fabs(LastTerm/Sum)<0.001 )
			lgNotConverged = FALSE;

		++i;

	}while ( lgNotConverged );

	return Sum;
}

static complex HyperGeoInt( double v )
{
	complex integrand;
			
	integrand = cmmul( cmmul( cmpow(ftocm(v,0.),BMinus1), cmpow(ftocm(1.-v,0.),CMinusBMinus1) ),
		cmpow(ftocm(1.-v*GlobalCHI,0.),MinusA) );

	return integrand;
}

/*complex 32 point Gaussian quadrature, originally given to Gary F by Jim Lattimer */
/* modified to handle complex numbers by Ryan Porter.	*/
static complex qg32complex(
	double xl, /*lower limit to integration range*/
	double xu, /*upper limit to integration range*/
	/*following is the pointer to the function that will be evaulated*/
	complex (*fct)(double) )
{
	double a, 
	  b, 
	  c;
	complex y;


	/********************************************************************************
	 *                                                                              *
	 *  32-point Gaussian quadrature                                                *
	 *  xl  : the lower limit of integration                                        *
	 *  xu  : the upper limit                                                       *
	 *  fct : the (external) function                                               *
	 *  returns the value of the integral                                           *
	 *                                                                              *
	 * simple call to integrate sine from 0 to pi                                   *
	 * double agn = qg32( 0., 3.141592654 ,  sin );                                 *
	 *                                                                              *
	 *******************************************************************************/

	a = .5*(xu + xl);
	b = xu - xl;
	c = .498631930924740780*b;
	y = cmmul( ftocm(.35093050047350483e-2,0.) , cmadd( (*fct)(a+c), (*fct)(a-c) ) );
	c = b*.49280575577263417;
	y = cmadd( y, cmmul( ftocm(.8137197365452835e-2,0.) , cmadd( (*fct)(a+c), (*fct)(a-c) ) ) );
	c = b*.48238112779375322;
	y = cmadd( y, cmmul( ftocm(.1269603265463103e-1,0.) , cmadd( (*fct)(a+c), (*fct)(a-c) ) ) );
	c = b*.46745303796886984;
	y = cmadd( y, cmmul( ftocm(.17136931456510717e-1,0.), cmadd( (*fct)(a+c), (*fct)(a-c) ) ) );
	c = b*.44816057788302606;
	y = cmadd( y, cmmul( ftocm(.21417949011113340e-1,0.), cmadd( (*fct)(a+c), (*fct)(a-c) ) ) );
	c = b*.42468380686628499;
	y = cmadd( y, cmmul( ftocm(.25499029631188088e-1,0.), cmadd( (*fct)(a+c), (*fct)(a-c) ) ) );
	c = b*.3972418979839712;
	y = cmadd( y, cmmul( ftocm(.29342046739267774e-1,0.), cmadd( (*fct)(a+c), (*fct)(a-c) ) ) );
	c = b*.36609105937014484;
	y = cmadd( y, cmmul( ftocm(.32911111388180923e-1,0.), cmadd( (*fct)(a+c), (*fct)(a-c) ) ) );
	c = b*.3315221334651076;
	y = cmadd( y, cmmul( ftocm(.36172897054424253e-1,0.), cmadd( (*fct)(a+c), (*fct)(a-c) ) ) );
	c = b*.29385787862038116;
	y = cmadd( y, cmmul( ftocm(.39096947893535153e-1,0.), cmadd( (*fct)(a+c), (*fct)(a-c) ) ) );
	c = b*.2534499544661147;
	y = cmadd( y, cmmul( ftocm(.41655962113473378e-1,0.), cmadd( (*fct)(a+c), (*fct)(a-c) ) ) );
	c = b*.21067563806531767;
	y = cmadd( y, cmmul( ftocm(.43826046502201906e-1,0.), cmadd( (*fct)(a+c), (*fct)(a-c) ) ) );
	c = b*.16593430114106382;
	y = cmadd( y, cmmul( ftocm(.45586939347881942e-1,0.), cmadd( (*fct)(a+c), (*fct)(a-c) ) ) );
	c = b*.11964368112606854;
	y = cmadd( y, cmmul( ftocm(.46922199540402283e-1,0.), cmadd( (*fct)(a+c), (*fct)(a-c) ) ) );
	c = b*.7223598079139825e-1;
	y = cmadd( y, cmmul( ftocm(.47819360039637430e-1,0.), cmadd( (*fct)(a+c), (*fct)(a-c) ) ) );
	c = b*.24153832843869158e-1;
	y = cmadd( y, cmmul( ftocm(.4827004425736390e-1,0.) , cmadd( (*fct)(a+c), (*fct)(a-c) ) ) );
	y = cmmul( y, ftocm( b, 0.) );

	/* the answer */

	return( y );
}

