/*********************************************************************
   Returns the Bessel functions All of them. Pick up our choice.
   C.A. Bertulani        May/16/2000
*********************************************************************/
typedef double Number;

#include <math.h>
#include<iostream.h>

int main(){
	Number bessj0(Number );
	Number bessj1( Number );
	Number bessj(int, Number );
	Number bessi0( Number );
	Number bessi1( Number );
	Number bessi(int, Number );
	Number bessk0( Number );
	Number bessk1( Number );
	Number bessk(int, Number );
	Number bessy0( Number );
	Number bessy1( Number );
	Number bessy(int, Number );
	void sphbes(int, Number , Number *, Number *, Number *, Number *);
	Number x, sbj, sby, dsbj, dsby;
	int j, n;
	j=1;
	cout << "\n\n Enter x.\n";
	cin >> x;
    for(;;){
		if(j>10){ cout << "\n My patience is over. Stop, please!\n";
			break;
		}
		if(j!=1){
			cout << "\n\n Enter x.\n";
			cin >> x;
		}
		cout << "\n J_0(x): " << bessj0(x);
		cout << "\n J_1(x): " << bessj1(x);
		cout << "\n J_5(x): " << bessj(5,x);
		cout << "\n I_0(x): " << bessi0(x);
		cout << "\n I_1(x): " << bessi1(x);
		cout << "\n I_5(x): " << bessi(5,x);
		cout << "\n K_0(x): " << bessk0(x);
		cout << "\n K_1(x): " << bessk1(x);
		cout << "\n K_5(x): " << bessk(5,x);
		cout << "\n Y_0(x): " << bessy0(x);
		cout << "\n Y_1(x): " << bessy1(x);
		cout << "\n Y_5(x): " << bessy(5,x);
		cout << "\n\n Enter n to calculate spherical Bessel function of order n\n";
		cin >> n;
		sphbes(n, x, &sbj, &sby, &dsbj, &dsby);
		cout << "\n sphbj_n(x): " << sbj;
		cout << "\n sphby_n(x): " << sby;
		cout << "\n d(sphbj_n)/dx: " << dsbj;
		cout << "\n d(sphby_n)/dx: " << dsby;
		j=j+1;
	}
	return 0;
}
/*********************************************************************
   Returns the Bessel function J_0(x) for any real x.
   C.A. Bertulani        May/16/2000
*********************************************************************/
Number bessj0(Number x)
{
	Number ax,z;
	Number xx,y,ans,ans1,ans2;		/* Acumulate polynomials in Number */
									/* precision. */

	if ((ax=fabs(x)) < 8.0) {		/* Direct rational function fit. */
		y=x*x;
		ans1=57568490574.0+y*(-13362590354.0+y*(651619640.7
			+y*(-11214424.18+y*(77392.33017+y*(-184.9052456)))));
		ans2=57568490411.0+y*(1029532985.0+y*(9494680.718
			+y*(59272.64853+y*(267.8532712+y*1.0))));
		ans=ans1/ans2;
	} else {		/* Fitting function. */
		z=8.0/ax;
		y=z*z;
		xx=ax-0.785398164;
		ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
			+y*(-0.2073370639e-5+y*0.2093887211e-6)));
		ans2 = -0.1562499995e-1+y*(0.1430488765e-3
			+y*(-0.6911147651e-5+y*(0.7621095161e-6
			-y*0.934935152e-7)));
		ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
	}
	return ans;
}
/*********************************************************************
   Returns the Bessel function J_1(x) for any real x.
   C.A. Bertulani        May/16/2000
*********************************************************************/
Number bessj1(Number x)
{
	Number ax,z;
	Number xx,y,ans,ans1,ans2;		/* Accumulate polynomials in Number */
									/* precision. */


	if ((ax=fabs(x)) < 8.0) {		/* Rational function approximation. */
		y=x*x;
		ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1
			+y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
		ans2=144725228442.0+y*(2300535178.0+y*(18583304.74
			+y*(99447.43394+y*(376.9991397+y*1.0))));
		ans=ans1/ans2;
	} else {
		z=8.0/ax;
		y=z*z;
		xx=ax-2.356194491;
		ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
			+y*(0.2457520174e-5+y*(-0.240337019e-6))));
		ans2=0.04687499995+y*(-0.2002690873e-3
			+y*(0.8449199096e-5+y*(-0.88228987e-6
			+y*0.105787412e-6)));
		ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
		if (x < 0.0) ans = -ans;
	}
	return ans;
}
/*********************************************************************
   Returns the Bessel function J_n(x) for any real x and n>=2.
   C.A. Bertulani        May/16/2000
*********************************************************************/
#define ACC 40.0		/* Make larger to increase accuracy. */
#define BIGNO 1.0e10
#define BIGNI 1.0e-10

Number bessj(int n, Number x)
{
	Number bessj0(Number x);
	Number bessj1(Number x);
	int j,jsum,m;
	Number ax,bj,bjm,bjp,sum,tox,ans;

	if (n < 2) cerr << "Index n less than 2 in bessj\n";
	ax=fabs(x);
	if (ax == 0.0)
		return 0.0;
	else if (ax > (Number) n) {		/* Upwards recurrence from J_0 and J_1 */
		tox=2.0/ax;
		bjm=bessj0(ax);
		bj=bessj1(ax);
		for (j=1;j<n;j++) {
			bjp=j*tox*bj-bjm;
			bjm=bj;
			bj=bjp;
		}
		ans=bj;
	} else {		/* Downwards recurrence from an even m here computed. */
		tox=2.0/ax;
		m=2*((n+(int) sqrt(ACC*n))/2);
		jsum=0;		/* jsum will alternate between 0 and 1; when it is 1, we */
		bjp=ans=sum=0.0;		/* accumulate in sum the even terms. */
		bj=1.0;
		for (j=m;j>0;j--) {		/* Downward recurrence. */
			bjm=j*tox*bj-bjp;
			bjp=bj;
			bj=bjm;
			if (fabs(bj) > BIGNO) {		/* Renormalize to prevent overflows. */
				bj *= BIGNI;
				bjp *= BIGNI;
				ans *= BIGNI;
				sum *= BIGNI;
			}
			if (jsum) sum += bj;		/* Accumulate the sum.  */
			jsum=!jsum;					/* Change 0 to 1 or vice versa. */
			if (j == n) ans=bjp;		/* Save the unnormalized answer. */
		}
		sum=2.0*sum-bj;
		ans /= sum;						/* Normalize the answer. */
	}
	return x < 0.0 && (n & 1) ? -ans : ans;
}
/*********************************************************************
   Returns the Bessel function I_0(x) for any real x.
   C.A. Bertulani        May/16/2000
*********************************************************************/
#undef ACC
#undef BIGNO
#undef BIGNI
Number bessi0(Number x)
{
	Number ax,ans;
	Number y;				/* Accumulate polynomials in Number precision.  */

	if ((ax=fabs(x)) < 3.75) {	/* Polynomial fit. */	
		y=x/3.75;
		y*=y;
		ans=1.0+y*(3.5156229+y*(3.0899424+y*(1.2067492
			+y*(0.2659732+y*(0.360768e-1+y*0.45813e-2)))));
	} else {
		y=3.75/ax;
		ans=(exp(ax)/sqrt(ax))*(0.39894228+y*(0.1328592e-1
			+y*(0.225319e-2+y*(-0.157565e-2+y*(0.916281e-2
			+y*(-0.2057706e-1+y*(0.2635537e-1+y*(-0.1647633e-1
			+y*0.392377e-2))))))));
	}
	return ans;
}
/*********************************************************************
   Returns the Bessel function K_0(x) for any positive real x.
   C.A. Bertulani        May/16/2000
*********************************************************************/
Number bessk0(Number x)
{
	Number bessi0(Number x);
	Number y,ans;			/* Accumulate polynomials in Number precision.  */

	if (x <= 2.0) {	/* Polynomial fit. */	
		y=x*x/4.0;
		ans=(-log(x/2.0)*bessi0(x))+(-0.57721566+y*(0.42278420
			+y*(0.23069756+y*(0.3488590e-1+y*(0.262698e-2
			+y*(0.10750e-3+y*0.74e-5))))));
	} else {
		y=2.0/x;
		ans=(exp(-x)/sqrt(x))*(1.25331414+y*(-0.7832358e-1
			+y*(0.2189568e-1+y*(-0.1062446e-1+y*(0.587872e-2
			+y*(-0.251540e-2+y*0.53208e-3))))));
	}
	return ans;
}
/*********************************************************************
   Returns the Bessel function I_1(x) for any real x.
   C.A. Bertulani        May/16/2000
*********************************************************************/
Number bessi1(Number x)
{
	Number ax,ans;
	Number y;			/* Accumulate polynomials in Number precision.  */

	if ((ax=fabs(x)) < 3.75) {	/* Polynomial fit. */	
		y=x/3.75;
		y*=y;
		ans=ax*(0.5+y*(0.87890594+y*(0.51498869+y*(0.15084934
			+y*(0.2658733e-1+y*(0.301532e-2+y*0.32411e-3))))));
	} else {
		y=3.75/ax;
		ans=0.2282967e-1+y*(-0.2895312e-1+y*(0.1787654e-1
			-y*0.420059e-2));
		ans=0.39894228+y*(-0.3988024e-1+y*(-0.362018e-2
			+y*(0.163801e-2+y*(-0.1031555e-1+y*ans))));
		ans *= (exp(ax)/sqrt(ax));
	}
	return x < 0.0 ? -ans : ans;
}
/*********************************************************************
   Returns the Bessel function K_1(x) for any positive real x.
   C.A. Bertulani        May/16/2000
*********************************************************************/
Number bessk1(Number x)
{
	Number bessi1(Number x);
	Number y,ans;			/* Accumulate polynomials in Number precision.  */

	if (x <= 2.0) {	/* Polynomial fit. */	
		y=x*x/4.0;
		ans=(log(x/2.0)*bessi1(x))+(1.0/x)*(1.0+y*(0.15443144
			+y*(-0.67278579+y*(-0.18156897+y*(-0.1919402e-1
			+y*(-0.110404e-2+y*(-0.4686e-4)))))));
	} else {
		y=2.0/x;
		ans=(exp(-x)/sqrt(x))*(1.25331414+y*(0.23498619
			+y*(-0.3655620e-1+y*(0.1504268e-1+y*(-0.780353e-2
			+y*(0.325614e-2+y*(-0.68245e-3)))))));
	}
	return ans;
}
/*********************************************************************
   Returns the Bessel function K_n(x) for any positive real x and n>=2.
   C.A. Bertulani        May/16/2000
*********************************************************************/
Number bessk(int n, Number x)
{
	Number bessk0(Number x);
	Number bessk1(Number x);
	int j;
	Number bk,bkm,bkp,tox;

	if (n < 2) cerr << "Index n less than 2 in bessk";
	tox=2.0/x;
	bkm=bessk0(x);		/* Upward recurrence for all x ...   */
	bk=bessk1(x);
	for (j=1;j<n;j++) {		/* ... and there it is.   */
		bkp=bkm+j*tox*bk;
		bkm=bk;
		bk=bkp;
	}
	return bk;
}
/*********************************************************************
   Returns the Bessel function I_n(x) for any real x and n>=2.
   C.A. Bertulani        May/16/2000
*********************************************************************/
#define ACC 40.0		/* Make larger to increase accuracy.  */
#define BIGNO 1.0e10
#define BIGNI 1.0e-10

Number bessi(int n, Number x)
{
	Number bessi0(Number x);
	int j;
	Number bi,bim,bip,tox,ans;

	if (n < 2) cerr << "Index n less than 2 in bessi\n";
	if (x == 0.0)
		return 0.0;
	else {
		tox=2.0/fabs(x);
		bip=ans=0.0;
		bi=1.0;
		for (j=2*(n+(int) sqrt(ACC*n));j>0;j--) {  /* Download recurrence from even m.  */
			bim=bip+j*tox*bi;
			bip=bi;
			bi=bim;
			if (fabs(bi) > BIGNO) {		/* Renormalize to prevent overflows. */
				ans *= BIGNI;
				bi *= BIGNI;
				bip *= BIGNI;
			}
			if (j == n) ans=bip;
		}
		ans *= bessi0(x)/bi;		/* Normalize with bessi0.  */
		return x < 0.0 && (n & 1) ? -ans : ans;
	}
}
#undef ACC
#undef BIGNO
#undef BIGNI
/*********************************************************************
   Returns the Bessel function Y_0(x) for any positive real x.
   C.A. Bertulani        May/16/2000
*********************************************************************/
Number bessy0(Number x)
{
	Number bessj0(Number x);
	Number z;
	Number xx,y,ans,ans1,ans2;	/* Accumulate polynomials in Number precision.  */

	if (x < 8.0) {		/* Direct rational function fit. */
		y=x*x;
		ans1 = -2957821389.0+y*(7062834065.0+y*(-512359803.6
			+y*(10879881.29+y*(-86327.92757+y*228.4622733))));
		ans2=40076544269.0+y*(745249964.8+y*(7189466.438
			+y*(47447.26470+y*(226.1030244+y*1.0))));
		ans=(ans1/ans2)+0.636619772*bessj0(x)*log(x);
	} else {		/* Fitting function. */
		z=8.0/x;
		y=z*z;
		xx=x-0.785398164;
		ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
			+y*(-0.2073370639e-5+y*0.2093887211e-6)));
		ans2 = -0.1562499995e-1+y*(0.1430488765e-3
			+y*(-0.6911147651e-5+y*(0.7621095161e-6
			+y*(-0.934945152e-7))));
		ans=sqrt(0.636619772/x)*(sin(xx)*ans1+z*cos(xx)*ans2);
	}
	return ans;
}
/*********************************************************************
   Returns the Bessel function Y_1(x) for any positive real x.
   C.A. Bertulani        May/16/2000
*********************************************************************/
Number bessy1(Number x)
{
	Number bessj1(Number x);
	Number z;
	Number xx,y,ans,ans1,ans2;	/* Accumulate polynomials in Number precision.  */

	if (x < 8.0) {		/* Rational function approximation.  */
		y=x*x;
		ans1=x*(-0.4900604943e13+y*(0.1275274390e13
			+y*(-0.5153438139e11+y*(0.7349264551e9
			+y*(-0.4237922726e7+y*0.8511937935e4)))));
		ans2=0.2499580570e14+y*(0.4244419664e12
			+y*(0.3733650367e10+y*(0.2245904002e8
			+y*(0.1020426050e6+y*(0.3549632885e3+y)))));
		ans=(ans1/ans2)+0.636619772*(bessj1(x)*log(x)-1.0/x);
	} else {			/* Fitting function. */
		z=8.0/x;
		y=z*z;
		xx=x-2.356194491;
		ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
			+y*(0.2457520174e-5+y*(-0.240337019e-6))));
		ans2=0.04687499995+y*(-0.2002690873e-3
			+y*(0.8449199096e-5+y*(-0.88228987e-6
			+y*0.105787412e-6)));
		ans=sqrt(0.636619772/x)*(sin(xx)*ans1+z*cos(xx)*ans2);
	}
	return ans;
}
/*********************************************************************
   Returns the Bessel function Y_n(x) for any positive real x and n>=2.
   C.A. Bertulani        May/16/2000
*********************************************************************/
Number bessy(int n, Number x)
{
	Number bessy0(Number x);
	Number bessy1(Number x);
	int j;
	Number by,bym,byp,tox;

	if (n < 2) cerr << "Index n less than 2 in bessy";
	tox=2.0/x;
	by=bessy1(x);			/* Starting values for recurrence. */
	bym=bessy0(x);
	for (j=1;j<n;j++) {		/* Recurrence. */
		byp=j*tox*by-bym;
		bym=by;
		by=byp;
	}
	return by;
}
/*********************************************************************
   Returns the spherical Bessel function j_n(x), yn(x), and their
   respective derivatives, for any positive real x and n>=0.
   For details, see Numerical Recipes, Chap. 6.	
  C.A. Bertulani        May/16/2000
*********************************************************************/

#define RTPIO2 1.2533141

void sphbes(int n, Number x, Number *sj, Number *sy, Number *sjp, Number *syp)
{
	void bessjy(Number x, Number xnu, Number *rj, Number *ry, Number *rjp,
		Number *ryp);
	Number factor,order,rj,rjp,ry,ryp;

	if (n < 0 || x <= 0.0) cerr << "bad arguments in sphbes\n";
	order=n+0.5;
	bessjy(x,order,&rj,&ry,&rjp,&ryp);
	factor=RTPIO2/sqrt(x);
	*sj=factor*rj;
	*sy=factor*ry;
	*sjp=factor*rjp-(*sj)/(2.0*x);
	*syp=factor*ryp-(*sy)/(2.0*x);
}
#undef RTPIO2
/*********************************************************************
   Needed to calculate spherical Bessel functions. Returns the Bessel
   fnctions rj=J_nu, ry = Y_nu and their derivatives rjp and ryp, for
   positive x and for xnu (nu) >=0. 
   For details, see Numerical Recipes in C, Chap. 6.		
   C.A. Bertulani        May/16/2000
*********************************************************************/
#include "butil.h"
#define EPS 1.0e-16
#define FPMIN 1.0e-30
#define MAXIT 10000
#define XMIN 2.0
#define PI 3.141592653589793

void bessjy(Number x, Number xnu, Number *rj, Number *ry, Number *rjp, Number *ryp)
{
	void beschb(Number x, Number *gam1, Number *gam2, Number *gampl,
		Number *gammi);
	int i,isign,l,nl;
	Number a,b,br,bi,c,cr,ci,d,del,del1,den,di,dlr,dli,dr,e,f,fact,fact2,
		fact3,ff,gam,gam1,gam2,gammi,gampl,h,p,pimu,pimu2,q,r,rjl,
		rjl1,rjmu,rjp1,rjpl,rjtemp,ry1,rymu,rymup,rytemp,sum,sum1,
		temp,w,x2,xi,xi2,xmu,xmu2;

	if (x <= 0.0 || xnu < 0.0) cerr << "bad arguments in bessjy\n";
	nl=(x < XMIN ? (int)(xnu+0.5) : IMAX(0,(int)(xnu-x+1.5)));
	/* n is the number of downward recurrences of the J's and upward recurrences */
	/* of Y's. xmu lies between -1/2 and 1/2 for x < XMIN, while it is chosen so */
	/* that x is greater than the turning point for x>= XMIN. */
	xmu=xnu-nl;
	xmu2=xmu*xmu;
	xi=1.0/x;
	xi2=2.0*xi;
	w=xi2/PI;			/* The Wronskian.  */
	isign=1;			/* Lentz's method */
	h=xnu*xi;			/* isign keeps track of sign changes in the denominator. */
	if (h < FPMIN) h=FPMIN;
	b=xi2*xnu;
	d=0.0;
	c=h;
	for (i=1;i<=MAXIT;i++) {
		b += xi2;
		d=b-d;
		if (fabs(d) < FPMIN) d=FPMIN;
		c=b-1.0/c;
		if (fabs(c) < FPMIN) c=FPMIN;
		d=1.0/d;
		del=c*d;
		h=del*h;
		if (d < 0.0) isign = -isign;
		if (fabs(del-1.0) < EPS) break;
	}
	if (i > MAXIT) cerr << "x too large in bessjy; try asymptotic expansion\n";
	rjl=isign*FPMIN;	/* Initialize J_nu and J'_nu for downward recurrence. */
	rjpl=h*rjl;
	rjl1=rjl;			/* Store values for later r`escaling. */
	rjp1=rjpl;
	fact=xnu*xi;
	for (l=nl;l>=1;l--) {
		rjtemp=fact*rjl+rjpl;
		fact -= xi;
		rjpl=fact*rjtemp-rjl;
		rjl=rjtemp;
	}
	if (rjl == 0.0) rjl=EPS;
	f=rjpl/rjl;
	if (x < XMIN) {	/* Now have unnormalized J_nu and J'_nu. */
		x2=0.5*x;	/* Use series. */
		pimu=PI*xmu;
		fact = (fabs(pimu) < EPS ? 1.0 : pimu/sin(pimu));
		d = -log(x2);
		e=xmu*d;
		fact2 = (fabs(e) < EPS ? 1.0 : sinh(e)/e);
		beschb(xmu,&gam1,&gam2,&gampl,&gammi);
		ff=2.0/PI*fact*(gam1*cosh(e)+gam2*fact2*d);
		e=exp(e);
		p=e/(gampl*PI);
		q=1.0/(e*PI*gammi);
		pimu2=0.5*pimu;
		fact3 = (fabs(pimu2) < EPS ? 1.0 : sin(pimu2)/pimu2);
		r=PI*pimu2*fact3*fact3;
		c=1.0;
		d = -x2*x2;
		sum=ff+r*q;
		sum1=p;
		for (i=1;i<=MAXIT;i++) {
			ff=(i*ff+p+q)/(i*i-xmu2);
			c *= (d/i);
			p /= (i-xmu);
			q /= (i+xmu);
			del=c*(ff+r*q);
			sum += del;
			del1=c*p-i*del;
			sum1 += del1;
			if (fabs(del) < (1.0+fabs(sum))*EPS) break;
		}
		if (i > MAXIT) cerr << "bessy series failed to converge\n";
		rymu = -sum;
		ry1 = -sum1*xi2;
		rymup=xmu*xi*rymu-ry1;
		rjmu=w/(rymup-f*rymu);
	} else {
		a=0.25-xmu2;
		p = -0.5*xi;
		q=1.0;
		br=2.0*x;
		bi=2.0;
		fact=a*xi/(p*p+q*q);
		cr=br+q*fact;
		ci=bi+p*fact;
		den=br*br+bi*bi;
		dr=br/den;
		di = -bi/den;
		dlr=cr*dr-ci*di;
		dli=cr*di+ci*dr;
		temp=p*dlr-q*dli;
		q=p*dli+q*dlr;
		p=temp;
		for (i=2;i<=MAXIT;i++) {
			a += 2*(i-1);
			bi += 2.0;
			dr=a*dr+br;
			di=a*di+bi;
			if (fabs(dr)+fabs(di) < FPMIN) dr=FPMIN;
			fact=a/(cr*cr+ci*ci);
			cr=br+cr*fact;
			ci=bi-ci*fact;
			if (fabs(cr)+fabs(ci) < FPMIN) cr=FPMIN;
			den=dr*dr+di*di;
			dr /= den;
			di /= -den;
			dlr=cr*dr-ci*di;
			dli=cr*di+ci*dr;
			temp=p*dlr-q*dli;
			q=p*dli+q*dlr;
			p=temp;
			if (fabs(dlr-1.0)+fabs(dli) < EPS) break;
		}
		if (i > MAXIT) cerr << "cf2 failed in bessjy\n";
		gam=(p-f)/q;
		rjmu=sqrt(w/((p-f)*gam+q));
		rjmu=SIGN(rjmu,rjl);
		rymu=rjmu*gam;
		rymup=rymu*(p+q/gam);
		ry1=xmu*xi*rymu-rymup;
	}
	fact=rjmu/rjl;
	*rj=rjl1*fact;		/* Scale original J_nu and Y_nu.  */
	*rjp=rjp1*fact;
	for (i=1;i<=nl;i++) {	/* Upward recurrence of Y_nu. */
		rytemp=(xmu+i)*xi2*ry1-rymu;
		rymu=ry1;
		ry1=rytemp;
	}
	*ry=rymu;
	*ryp=xnu*xi*rymu-ry1;
}
#undef EPS
#undef FPMIN
#undef MAXIT
#undef XMIN
#undef PI
/*********************************************************************
   Needed to calculate spherical Bessel function j_n(x).
   For details, see Numerical Recipes in C, Chap. 6.	
   C.A. Bertulani        May/16/2000
*********************************************************************/
#define NUSE1 7
#define NUSE2 8

void beschb(Number x, Number *gam1, Number *gam2, Number *gampl, Number *gammi)
{
	Number chebev(Number a, Number b, Number c[], int m, Number x);
	Number xx;
	static Number c1[] = {
		-1.142022680371172e0,6.516511267076e-3,
		3.08709017308e-4,-3.470626964e-6,6.943764e-9,
		3.6780e-11,-1.36e-13};
	static Number c2[] = {
		1.843740587300906e0,-0.076852840844786e0,
		1.271927136655e-3,-4.971736704e-6,-3.3126120e-8,
		2.42310e-10,-1.70e-13,-1.0e-15};

	xx=8.0*x*x-1.0;
	*gam1=chebev(-1.0,1.0,c1,NUSE1,xx);
	*gam2=chebev(-1.0,1.0,c2,NUSE2,xx);
	*gampl= *gam2-x*(*gam1);
	*gammi= *gam2+x*(*gam1);
}
#undef NUSE1
#undef NUSE2
/*********************************************************************
   Needed to calculate spherical Bessel function j_n(x).
   For details, see Numerical Recipes in C, Chap. 6.	
   C.A. Bertulani        May/16/2000
*********************************************************************/
Number chebev(Number a, Number b, Number c[], int m, Number x)
{
	Number d=0.0,dd=0.0,sv,y,y2;
	int j;

	if ((x-a)*(x-b) > 0.0) cerr << "x not in range in routine chebev";
	y2=2.0*(y=(2.0*x-a-b)/(b-a));
	for (j=m-1;j>=1;j--) {
		sv=d;
		d=y2*d-dd+c[j];
		dd=sv;
	}
	return y*d-dd+0.5*c[0];
}
/*********************************************************************/