/* 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 */
/*ftocm convert floats to dp complex */
/*cmneg dp complex unary minus operation */
/*cmabs dp complex absolute value */
/*cmadd retn the sum of two dp complex nos. */
/*cmsub retn the difference of two dp complex nos. */
/*cmmul retn the product of two dp complex nos. */
/*cmdiv retn the quotient of two dp complex nos. */
/*cmexp dp complex exponential */
#include "cddefines.h"
#include "math_complex.h"

/*ftocm convert floats to dp complex */
complex ftocm(double r,double i )	
{
	complex dz;
	dz.re = r;
	dz.im = i;
	return( dz );
}

/*cmneg dp complex unary minus operation */
complex cmneg(complex dz )	
{
	dz.re = -dz.re;
	dz.im = -dz.im;
	return( dz );
}

/*cmabs dp complex absolute value */
double cmabs(complex dz )	
{
	double s;

	dz.re = fabs(dz.re);
	dz.im = fabs(dz.im);

	if( dz.re == 0. )
	{
		return dz.im;
	}
	else if ( dz.im == 0. )
	{
		return dz.re;
	}
	else if( fabs( log10(dz.re) - log10(dz.im) ) > 15. )
	{
		return fabs( MAX2(dz.re, dz.im) );
	}
	else if( dz.re<1e-100 || dz.re>1e100 )
	{
		return dz.re * sqrt( 1. + (dz.im/dz.re)*(dz.im/dz.re) );
	}
	else
	{
		s = dz.re*dz.re + dz.im*dz.im;
		return( sqrt(s) );
	}
}

/*cmadd retn the sum of two dp complex nos. */
complex cmadd(complex l,complex r )	
{
	l.re += r.re;
	l.im += r.im;
	return( l );
}

/*cmsub retn the difference of two dp complex nos. */
complex cmsub(complex l,complex r )	
{
	l.re -= r.re;
	l.im -= r.im;
	return( l );
}


/*cmmul retn the product of two dp complex nos. */
complex cmmul(complex l,complex r )	
{
	complex dz;
	dz.re = l.re*r.re - l.im*r.im;
	dz.im = l.re*r.im + l.im*r.re;
	return( dz );
}

/*cmdiv retn the quotient of two dp complex nos. */
complex cmdiv(complex l,complex r )	
{
	complex dz;
	double den;

	if( r.re == 0. && r.im == 0. )
	{
		fprintf( ioQQQ, " Complex division by 0." );
	}

	den = r.re*r.re + r.im*r.im;

	/* chng 07 june 02, by Ryan. I changed the order of operations 
	 * in order to prevent overflows and underflows when r is very small or big	*/
	
	/*dz.re = ( l.re*r.re + l.im*r.im )/den;
	dz.im = ( r.re*l.im - l.re*r.im )/den;	*/

	dz.re = l.re/den*r.re + l.im/den*r.im;
	dz.im = r.re/den*l.im - l.re/den*r.im;

	return( dz );
}

/*cmexp dp complex exponential */
complex cmexp(complex dz )	
{
	double exp_re;

	if( dz.re == 0. )
	{
		exp_re = 1.;
	}
	else
	{
		exp_re = exp(dz.re);
	}

	if( dz.im == 0.0 ) 
	{	
		/* real only complex no. */
		dz.re = exp_re;
		return( dz );
	}
	dz.re = exp_re*cos(dz.im);
	dz.im = exp_re*sin(dz.im);

	return( dz );
}

/*********************************************
 *
 * NB - the functions in this module were originally written
 *      for single precision arithmatic, and no attempt has
 *      been made to upgrade the accuracy of the return values
 *      to full machine precision !
 *
 *      The precision to which the result is calculated is more
 *      or less governed by MACHEPF (at least for expnf). Setting
 *      it to a lesser precision should speed up the routine.
 *
 *      expnf now uses gammafun to calculate factorials of
 *      integer numbers. the routine could be speeded up by
 *      using a smarter algorithm there.
 *
 *********************************************/


static double polevlf( double x, double *coef, int N );
/*@unused@ static double p1evlf( double x, double *coef, int N );*/
static double stirf( double x );

/* chng 28 may 02, by Ryan...moved this to math.h	*/
/*static double gammafun( double x );*/
#if 0
static double lgamf( double x );
#endif


#define EUL 0.57721566490153286060
const double MACHEPF = FLT_EPSILON;
const double PIF = 3.141592653589793238;

/*#include "mconf.h"*/
/*							polevlf.c
 *							p1evlf.c
 *
 *	Evaluate polynomial
 *
 *
 *
 * SYNOPSIS:
 *
 * int N;
 * float x, y, coef[N+1], polevlf[];
 *
 * y = polevlf( x, coef, N );
 *
 *
 *
 * DESCRIPTION:
 *
 * Evaluates polynomial of degree N:
 *
 *                     2          N
 * y  =  C  + C x + C x  +...+ C x
 *        0    1     2          N
 *
 * Coefficients are stored in reverse order:
 *
 * coef[0] = C  , ..., coef[N] = C  .
 *            N                   0
 *
 *  The function p1evl() assumes that coef[N] = 1.0 and is
 * omitted from the array.  Its calling arguments are
 * otherwise the same as polevl().
 *
 *
 * SPEED:
 *
 * In the interest of speed, there are no checks for out
 * of bounds arithmetic.  This routine is used by most of
 * the functions in the library.  Depending on available
 * equipment features, the user may wish to rewrite the
 * program in microcode or assembly language.
 *
 */

/*
Cephes Math Library Release 2.1:  December, 1988
Copyright 1984, 1987, 1988 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/

/*#include "mconf.h"*/

static double polevlf( double x, double *coef, int N )
{
	double ans;
	double *p;
	int i;

	p = coef;
	ans = *p++;

	/*
	  for( i=0; i<N; i++ )
	  ans = ans * x  +  *p++;
	*/

	i = N;
	do
		ans = ans * x  +  *p++;
	while( --i );

	return( ans );
}

/*							p1evl()	*/
/*                                          N
 * Evaluate polynomial when coefficient of x  is 1.0.
 * Otherwise same as polevl.
 */

#if 0
double p1evlf( double x, double *coef, int N )
{
	double ans;
	double *p;
	int i;

	p = coef;
	ans = x + *p++;

	i = N-1;
	do
		ans = ans * x  + *p++;
	while( --i );

	return( ans );
}
#endif
/*							gammafun.c
 *
 *	Gamma function
 *
 *
 *
 * SYNOPSIS:
 *
 * double x, y, gammafun();
 * extern int sgngamf;
 *
 * y = gammafun( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns gamma function of the argument.  The result is
 * correctly signed, and the sign (+1 or -1) is also
 * returned in a global (extern) variable named sgngamf.
 * This same variable is also filled in by the logarithmic
 * gamma function lgam().
 *
 * Arguments between 0 and 10 are reduced by recurrence and the
 * function is approximated by a polynomial function covering
 * the interval (2,3).  Large arguments are handled by Stirling's
 * formula. Negative arguments are made positive using
 * a reflection formula.  
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE       0,-33      100,000     5.7e-7      1.0e-7
 *    IEEE       -33,0      100,000     6.1e-7      1.2e-7
 *
 *
 */

/*							lgamf()
 *
 *	Natural logarithm of gamma function
 *
 *
 *
 * SYNOPSIS:
 *
 * double x, y, lgamf();
 * extern int sgngamf;
 *
 * y = lgamf( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns the base e (2.718...) logarithm of the absolute
 * value of the gamma function of the argument.
 * The sign (+1 or -1) of the gamma function is returned in a
 * global (extern) variable named sgngamf.
 *
 * For arguments greater than 6.5, the logarithm of the gamma
 * function is approximated by the logarithmic version of
 * Stirling's formula.  Arguments between 0 and +6.5 are reduced by
 * by recurrence to the interval [.75,1.25] or [1.5,2.5] of a rational
 * approximation.  The cosecant reflection formula is employed for
 * arguments less than zero.
 *
 * Arguments greater than MAXLGM = 2.035093e36 return MAXNUM and an
 * error message.
 *
 *
 *
 * ACCURACY:
 *
 *
 *
 * arithmetic      domain        # trials     peak         rms
 *    IEEE        -100,+100       500,000    7.4e-7       6.8e-8
 * The error criterion was relative when the function magnitude
 * was greater than one but absolute when it was less than one.
 * The routine has low relative error for positive arguments.
 *
 * The following test used the relative error criterion.
 *    IEEE    -2, +3              100000     4.0e-7      5.6e-8
 *
 */

/*							gamma.c	*/
/*	gamma function	*/

/*
Cephes Math Library Release 2.7:  July, 1998
Copyright 1984, 1987, 1989, 1992, 1998 by Stephen L. Moshier
*/


/* define MAXGAM 34.84425627277176174 */

/* Stirling's formula for the gamma function
 * gamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) ( 1 + 1/x P(1/x) )
 * .028 < 1/x < .1
 * relative error < 1.9e-11
 */
static double STIR[] = {
-2.705194986674176E-003,
 3.473255786154910E-003,
 8.333331788340907E-002,
};

/* choose this number such that MAXSTIR^(MAXSTIR-0.5) = MAX_FLT */
static double MAXSTIR = 26.77;
static double SQTPIF = 2.50662827463100050242; /* sqrt( 2 pi ) */

int sgngamf = 0;

/*float polevlf( float, float *, int );
float p1evlf( float, float *, int );*/

/* Gamma function computed by Stirling's formula,
 * sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x))
 * The polynomial STIR is valid for 33 <= x <= 172.
 */
static double stirf( double x )
{
	double y, w, v;

	w = 1.0/x;
	w = 1.0 + w * polevlf( w, STIR, 2 );
	y = exp( -x );
	if( x > MAXSTIR )
	{ /* Avoid overflow in pow() */
		v = pow( x, 0.5 * x - 0.25 );
		y *= v;
		y *= v;
	}
	else
	{
		y = pow( x, x - 0.5 ) * y;
	}
	y = SQTPIF * y * w;
	return( y );
}


/* gamma(x+2), 0 < x < 1 */
static double P[] = {
 1.536830450601906E-003,
 5.397581592950993E-003,
 4.130370201859976E-003,
 7.232307985516519E-002,
 8.203960091619193E-002,
 4.117857447645796E-001,
 4.227867745131584E-001,
 9.999999822945073E-001,
};

double gammafun( double x )
{
	double p, q, z, nz;
	int lgInvert, lgNegative;

	sgngamf = 1;
	lgNegative = FALSE;
	nz = 0.0;
	if( x < 0.0 )
	{
		lgNegative = TRUE;
		q = -x;
		p = floor(q);
		/*lint -e777 float equal */
		if( p == q ) 
		/*lint +e777 float equal */
		{
			printf(" gammafun finds domain error\n");
			cdEXIT(EXIT_FAILURE);
		}
		if( (int)p%2 == 0 )
			sgngamf = -1;
		nz = q - p;

		if( nz > 0.5 )
		{
			p += 1.0;
			nz = q - p;
		}

		nz = q * sin( PIF * nz );
		if( nz == 0.0 )
		{
			printf(" gammafun finds domain error\n");
			cdEXIT(EXIT_FAILURE);
		}
		nz = fabs(nz);
		x = q;
	}
	if( x >= 10.0 )
	{
		z = stirf(x);
	}

	lgInvert = ( x < 2.0 );

	z = 1.0;
	while( x >= 3.0 )
	{
		x -= 1.0;
		z *= x;
	}
	/*
	while( x < 0.0 )
		{
		if( x > -1.E-4 )
		{
			if( x == 0.0 )
			{
				printf(" gammafun finds domain error\n");
				cdEXIT(EXIT_FAILURE);
			}
			else
			{
				p = z / ((1.0 + EUL * x) * x);
				if( lgNegative )
				{
					p = sgngamf * PIF/(nz * p );
				}
				return(p);
			}
		}
		z *=x;
		x += 1.0;
		}
	*/
	while( x < 2.0 )
	{
		if( x < 1.e-4 )
		{
			if( x == 0.0 )
			{
				printf(" gammafun finds domain error\n");
				cdEXIT(EXIT_FAILURE);
			}
			else
			{
				p = z / ((1.0 + EUL * x) * x);
				if( lgNegative )
				{
					p = sgngamf * PIF/(nz * p );
				}
				return(p);
			}
		}
		z *=x;
		x += 1.0;
	}

	if( lgInvert )
		z = 1.0/z;

	if( x == 2.0 )
		return(z);

	x -= 2.0;
	p = z * polevlf( x, P, 7 );

	if( lgNegative )
	{
		p = sgngamf * PIF/(nz * p );
	}
	return(p);
}




/* log gamma(x+2), -.5 < x < .5 */
/*static double B[] = {
 6.055172732649237E-004,
-1.311620815545743E-003,
 2.863437556468661E-003,
-7.366775108654962E-003,
 2.058355474821512E-002,
-6.735323259371034E-002,
 3.224669577325661E-001,
 4.227843421859038E-001
};*/

/* log gamma(x+1), -.25 < x < .25 */
/*static double C[] = {
 1.369488127325832E-001,
-1.590086327657347E-001,
 1.692415923504637E-001,
-2.067882815621965E-001,
 2.705806208275915E-001,
-4.006931650563372E-001,
 8.224670749082976E-001,
-5.772156501719101E-001
};*/

/* log( sqrt( 2*pi ) ) */
/*static double LS2PI  =  0.91893853320467274178;*/
/*#define MAXLGM 2.035093e36*/
/*static double PIINV =  0.318309886183790671538;*/

/* Logarithm of gamma function */

#if 0
static double lgamf( double x )
{
	double p, q, w, z;
	double nx, tx;
	int i, direction;

	sgngamf = 1;

	if( x < 0.0 )
	{
		q = -x;
		w = lgamf(q); /* note this modifies sgngam! */
		p = floor(q);
		/*lint -e777 float equal */
		if( p == q )
		/*lint +e777 float equal */
			goto loverf;
		i = (int)p;
		if( (i & 1) == 0 )
			sgngamf = -1;
		else
			sgngamf = 1;
		z = q - p;
		if( z > 0.5 )
		{
			p += 1.0;
			z = p - q;
		}
		z = q * sin( PIF * z );
		if( z == 0.0 )
			goto loverf;
		z = -log( PIINV*z ) - w;
		return( z );
	}

	if( x < 6.5 )
	{
		direction = 0;
		z = 1.0;
		tx = x;
		nx = 0.0;
		if( x >= 1.5 )
		{
			while( tx > 2.5 )
			{
				nx -= 1.0;
				tx = x + nx;
				z *=tx;
			}
			x += nx - 2.0;
iv1r5:
			p = x * polevlf( x, B, 7 );
			goto cont;
		}
		if( x >= 1.25 )
		{
			z *= x;
			x -= 1.0; /* x + 1 - 2 */
			direction = 1;
			goto iv1r5;
		}
		if( x >= 0.75 )
		{
			x -= 1.0;
			p = x * polevlf( x, C, 7 );
			q = 0.0;
			goto contz;
		}
		while( tx < 1.5 )
		{
			if( tx == 0.0 )
				goto loverf;
			z *=tx;
			nx += 1.0;
			tx = x + nx;
		}
		direction = 1;
		x += nx - 2.0;
		p = x * polevlf( x, B, 7 );

cont:
		if( z < 0.0 )
		{
			sgngamf = -1;
			z = -z;
		}
		else
		{
			sgngamf = 1;
		}
		q = log(z);
		if( direction )
			q = -q;
contz:
		return( p + q );
	}

	if( x > MAXLGM )
	{
loverf:
		printf(" lgamf finds domain error\n");
		cdEXIT(EXIT_FAILURE);
	}

	/* Note, though an asymptotic formula could be used for x >= 3,
	 * there is cancellation error in the following if x < 6.5.  */
	q = LS2PI - x;
	q += ( x - 0.5 ) * log(x);

	if( x <= 1.0e4 )
	{
		z = 1.0/x;
		p = z * z;
		q += ((    6.789774945028216E-004 * p
			 - 2.769887652139868E-003 ) * p
			+  8.333316229807355E-002 ) * z;
	}
	return( q );
}
#endif


/*							expnf.c
 *
 *		Exponential integral En
 *
 *
 *
 * SYNOPSIS:
 *
 * int n;
 * double x, y, expnf();
 *
 * y = expnf( n, x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Evaluates the exponential integral
 *
 *                 inf.
 *                   -
 *                  | |   -xt
 *                  |    e
 *      E (x)  =    |    ----  dt.
 *       n          |      n
 *                | |     t
 *                 -
 *                  1
 *
 *
 * Both n and x must be nonnegative.
 *
 * The routine employs either a power series, a continued
 * fraction, or an asymptotic formula depending on the
 * relative values of n and x.
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      0, 30       10000       5.6e-7      1.2e-7
 *
 */

/*							expn.c	*/

/* Cephes Math Library Release 2.2:  July, 1992
 * Copyright 1985, 1992 by Stephen L. Moshier
 * Direct inquiries to 30 Frost Street, Cambridge, MA 02140 */


/* this number is 1./MACHEPF */
#define BIG   (1./FLT_EPSILON)


double expnf( int n, double x )
{
	double ans, r, t, yk, xk;
	double pk, pkm1, pkm2, qk, qkm1, qkm2;
	double psi, z;
	int i, k;


	if( n < 0 || x < 0. )
	{
		printf(" expnf finds domain error\n");
		cdEXIT(EXIT_FAILURE);
	}

	if( x > log(FLT_MAX) )
	{
		return( 0.0 );
	}

	if( x == 0.0 )
	{
		if( n < 2 )
		{
			printf(" expnf finds domain error\n");
			cdEXIT(EXIT_FAILURE);
		}
		else
		{
			return( 1.0/((double)n-1.0) );
		}
	}

	if( n == 0 )
	{
		return( exp(-x)/x );
	}

	/*							expn.c	*/
	/*		Expansion for large n		*/

	if( n > 5000 )
	{
		xk = x + n;
		yk = 1.0 / (xk * xk);
		t = (double)n;
		ans = yk * t * (6.0 * x * x  -  8.0 * t * x  +  t * t);
		ans = yk * (ans + t * (t  -  2.0 * x));
		ans = yk * (ans + t);
		ans = (ans + 1.0) * exp( -x ) / xk;
		return(ans);
	}

	if( x <= 1.0 )
	{
		/*							expn.c	*/
		/*		Power series expansion		*/

		psi = -EUL - log(x);
		for( i=1; i<n; i++ )
			psi = psi + 1.0/i;

		z = -x;
		xk = 0.0;
		yk = 1.0;
		pk = 1.0 - n;

		if( n == 1 )
			ans = 0.0;
		else
			ans = 1.0/pk;

		do
		{
			xk += 1.0;
			yk *= z/xk;
			pk += 1.0;
			if( pk != 0.0 )
			{
				ans += yk/pk;
			}
			if( ans != 0.0 )
				t = fabs(yk/ans);
			else
				t = 1.0;
		}
		while( t > MACHEPF );

		/* the following statement is used nowhere ? */
/* 		k = (int)xk; */
		t = (double)n;
		/* using gammafun is wasteful, it only calculates
		 * the factorial of an integer argument ! */
		ans = (powi(z,n-1) * psi / gammafun(t)) - ans;
		return(ans);
	}
	else
	{
		/*							expn.c	*/
		/*		continued fraction		*/
		k = 1;
		pkm2 = 1.0;
		qkm2 = x;
		pkm1 = 1.0;
		qkm1 = x + n;
		ans = pkm1/qkm1;

		do
		{
			k += 1;
			if( k%2 == 1 )
			{
				/* The previous if guarantees that the following two integer
				 * divides never result in roundoff error */
				/*lint -e653 */
				yk = 1.0;
				xk = (double)(n + (k-1)/2);
			}
			else
			{
				yk = x;
				xk = (double)(k/2);
				/*lint +e653 */
			}
			pk = pkm1 * yk  +  pkm2 * xk;
			qk = qkm1 * yk  +  qkm2 * xk;
			if( qk != 0 )
			{
				r = pk/qk;
				t = fabs( (ans - r)/r );
				ans = r;
			}
			else
			{
				t = 1.0;
			}

			pkm2 = pkm1;
			pkm1 = pk;
			qkm2 = qkm1;
			qkm1 = qk;

			if( fabs(pk) > BIG )
			{
				pkm2 *= MACHEPF;
				pkm1 *= MACHEPF;
				qkm2 *= MACHEPF;
				qkm1 *= MACHEPF;
			}
		}
		while( t > MACHEPF );
		ans *= exp( -x );
		return( ans );
	}
}

#undef BIG
