/*ee1 first exponential integral */
#include "cddefines.h"
#include "powi.h"
#include "ee1.h"
#include "sexp.h"

double ee1(double x)
{
	double ans, 
	  bot, 
	  ee1_v, 
	  top;
	static double a[6]={-.57721566,.99999193,-.24991055,.05519968,-.00976004,
	  .00107857};
	static double b[4]={8.57332874,18.059016973,8.634760893,.267773734};
	static double c[4]={9.573322345,25.632956149,21.099653083,3.95849692};

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

	/* computes the exponential integral E1(x),
	 * from abramowitz and stegun
	 * returns error condition for negative argument,
	 * returns zero when limit on sexp exceeded
	 * */

	/* error - this does not do complex numbers */
	if( x <= 0 )
	{
		fprintf( ioQQQ, " negative argument in function ee1, x<0\n" );
		puts( "[Stop in ee1]" );
		exit(1);
	}

	/* branch for x less than 1 */
	else if( x < 1. )
	{
		ans = a[0] + a[1]*x + a[2]*x*x + a[3]*powi(x,3) + a[4]*powi(x,4) + 
		  a[5]*powi(x,5) - log(x);
	}

	/* branch for x greater than, or equal to, one */
	else
	{
		top = powi(x,4) + b[0]*powi(x,3) + b[1]*x*x + b[2]*x + b[3];
		bot = powi(x,4) + c[0]*powi(x,3) + c[1]*x*x + c[2]*x + c[3];
		/* 
		 *>>>chng 98 dec 17
		 * Jason's original implementation did not require the exp since
		 * the routines that used this knew if was not present.  
		 * put in for safely.
		 
		ans = top/bot/x;*/
		ans = top/bot/x*sexp(x);
	}

	ee1_v = ans;

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