/*********************************************************************
   Computes Airy fnctions Ai(x), Bi(x), and their derivatives,
   for any real x.
   C.A. Bertulani        May/16/2000
*********************************************************************/
typedef double Number;

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

int main(){
	void airy(Number x, Number *, Number *, Number *aip, Number *);
	Number x, ai, bi, aip, bip;
	int j;
	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;
		}
		airy(x, &ai, &bi, &aip, &bip);
		cout << "\n Ai(x): " << ai;
		cout << "\n Bi(x): " << bi;
		cout << "\n d(Ai)/dx: " << aip;
		cout << "\n d(Bi)/dx: " << bip;
		j=j+1;
	}
	return 0;
}
/*********************************************************************
   Returns Airy fnctions Ai(x), Bi(x), and their derivatives,
   for any real x.
   C.A. Bertulani        May/16/2000
*********************************************************************/
#define PI 3.1415927
#define THIRD (1.0/3.0)
#define TWOTHR (2.0*THIRD)
#define ONOVRT 0.57735027

void airy(Number x, Number *ai, Number *bi, Number *aip, Number *bip)
{
	void bessik(Number x, Number xnu, Number *ri, Number *rk, Number *rip,
		Number *rkp);
	void bessjy(Number x, Number xnu, Number *rj, Number *ry, Number *rjp,
		Number *ryp);
	Number absx,ri,rip,rj,rjp,rk,rkp,rootx,ry,ryp,z;

	absx=fabs(x);
	rootx=sqrt(absx);
	z=TWOTHR*absx*rootx;
	if (x > 0.0) {
		bessik(z,THIRD,&ri,&rk,&rip,&rkp);
		*ai=rootx*ONOVRT*rk/PI;
		*bi=rootx*(rk/PI+2.0*ONOVRT*ri);
		bessik(z,TWOTHR,&ri,&rk,&rip,&rkp);
		*aip = -x*ONOVRT*rk/PI;
		*bip=x*(rk/PI+2.0*ONOVRT*ri);
	} else if (x < 0.0) {
		bessjy(z,THIRD,&rj,&ry,&rjp,&ryp);
		*ai=0.5*rootx*(rj-ONOVRT*ry);
		*bi = -0.5*rootx*(ry+ONOVRT*rj);
		bessjy(z,TWOTHR,&rj,&ry,&rjp,&ryp);
		*aip=0.5*absx*(ONOVRT*ry+rj);
		*bip=0.5*absx*(ONOVRT*rj-ry);
	} else {
		*ai=0.35502805;
		*bi=(*ai)/ONOVRT;
		*aip = -0.25881940;
		*bip = -(*aip)/ONOVRT;
	}
}
#undef PI
#undef THIRD
#undef TWOTHR
#undef ONOVRT
/*********************************************************************
   Needed to calculate Airy fnctions. Returns the Bessel
   functions ri=I_nu, rk = K_nu and their derivatives rip and rkp, for
   positive x and for xnu (nu) >=0. For details, see Numerical
   Recipes, chapter 6.
   C.A. Bertulani        May/16/2000
*********************************************************************/
#define EPS 1.0e-16
#define FPMIN 1.0e-30
#define MAXIT 10000
#define XMIN 2.0
#define PI 3.141592653589793

void bessik(Number x, Number xnu, Number *ri, Number *rk, Number *rip, Number *rkp)
{
	void beschb(Number x, Number *gam1, Number *gam2, Number *gampl,
		Number *gammi);
	int i,l,nl;
	Number a,a1,b,c,d,del,del1,delh,dels,e,f,fact,fact2,ff,gam1,gam2,
		gammi,gampl,h,p,pimu,q,q1,q2,qnew,ril,ril1,rimu,rip1,ripl,
		ritemp,rk1,rkmu,rkmup,rktemp,s,sum,sum1,x2,xi,xi2,xmu,xmu2;

	if (x <= 0.0 || xnu < 0.0) cerr << "bad arguments in bessik\n";
	nl=(int)(xnu+0.5);
	xmu=xnu-nl;
	xmu2=xmu*xmu;
	xi=1.0/x;
	xi2=2.0*xi;
	h=xnu*xi;
	if (h < FPMIN) h=FPMIN;
	b=xi2*xnu;
	d=0.0;
	c=h;
	for (i=1;i<=MAXIT;i++) {
		b += xi2;
		d=1.0/(b+d);
		c=b+1.0/c;
		del=c*d;
		h=del*h;
		if (fabs(del-1.0) < EPS) break;
	}
	if (i > MAXIT) cerr << "x too large in bessik; try asymptotic expansion\n";
	ril=FPMIN;
	ripl=h*ril;
	ril1=ril;
	rip1=ripl;
	fact=xnu*xi;
	for (l=nl;l>=1;l--) {
		ritemp=fact*ril+ripl;
		fact -= xi;
		ripl=fact*ritemp+ril;
		ril=ritemp;
	}
	f=ripl/ril;
	if (x < XMIN) {
		x2=0.5*x;
		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=fact*(gam1*cosh(e)+gam2*fact2*d);
		sum=ff;
		e=exp(e);
		p=0.5*e/gampl;
		q=0.5/(e*gammi);
		c=1.0;
		d=x2*x2;
		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;
			sum += del;
			del1=c*(p-i*ff);
			sum1 += del1;
			if (fabs(del) < fabs(sum)*EPS) break;
		}
		if (i > MAXIT) cerr << "bessk series failed to converge\n";
		rkmu=sum;
		rk1=sum1*xi2;
	} else {
		b=2.0*(1.0+x);
		d=1.0/b;
		h=delh=d;
		q1=0.0;
		q2=1.0;
		a1=0.25-xmu2;
		q=c=a1;
		a = -a1;
		s=1.0+q*delh;
		for (i=2;i<=MAXIT;i++) {
			a -= 2*(i-1);
			c = -a*c/i;
			qnew=(q1-b*q2)/a;
			q1=q2;
			q2=qnew;
			q += c*qnew;
			b += 2.0;
			d=1.0/(b+a*d);
			delh=(b*d-1.0)*delh;
			h += delh;
			dels=q*delh;
			s += dels;
			if (fabs(dels/s) < EPS) break;
		}
		if (i > MAXIT) cerr << "bessik: failure to converge in cf2\n";
		h=a1*h;
		rkmu=sqrt(PI/(2.0*x))*exp(-x)/s;
		rk1=rkmu*(xmu+x+0.5-h)*xi;
	}
	rkmup=xmu*xi*rkmu-rk1;
	rimu=xi/(f*rkmu-rkmup);
	*ri=(rimu*ril1)/ril;
	*rip=(rimu*rip1)/ril;
	for (i=1;i<=nl;i++) {
		rktemp=(xmu+i)*xi2*rk1+rkmu;
		rkmu=rk1;
		rk1=rktemp;
	}
	*rk=rkmu;
	*rkp=xnu*xi*rkmu-rk1;
}
#undef EPS
#undef FPMIN
#undef MAXIT
#undef XMIN
#undef PI
/*********************************************************************
   Needed to calculate Airy functions. Returns the Bessel
   functions 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 Airy functions.
   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 Airy functions.
   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];
}
/*********************************************************************/