using System;

namespace Science.Mathematics.SpecialFunction
{
	/// <summary>
	/// Bessel
	/// </summary>
	public class Bessel
	{
		public Bessel()
		{
		}
		public static double J0(double x)
		{
			double ax,z;
			double xx,y,ans,ans1,ans2;

			if ((ax=Math.Abs(x)) < 8.0) 
			{
				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 
			{
				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=Math.Sqrt(0.636619772/ax)*(Math.Cos(xx)*ans1-z*Math.Sin(xx)*ans2);
			}
			return ans;
		}
		public static double Y0(double x)
		{
			double z;
			double xx,y,ans,ans1,ans2;

			if (x < 8.0) 
			{
				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*J0(x)*Math.Log(x);
			} 
			else 
			{
				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=Math.Sqrt(0.636619772/x)*(Math.Sin(xx)*ans1+z*Math.Cos(xx)*ans2);
			}
			return ans;
		}
		public static double J1(double x)
		{
			double ax,z;
			double xx,y,ans,ans1,ans2;

			if ((ax=Math.Abs(x)) < 8.0) 
			{
				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=Math.Sqrt(0.636619772/ax)*(Math.Cos(xx)*ans1-z*Math.Sin(xx)*ans2);
				if (x < 0.0) ans = -ans;
			}
			return ans;
		}
		public static double Y1(double x)
		{
			double z;
			double xx,y,ans,ans1,ans2;

			if (x < 8.0) 
			{
				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*(J1(x)*Math.Log(x)-1.0/x);
			} 
			else 
			{
				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=Math.Sqrt(0.636619772/x)*(Math.Sin(xx)*ans1+z*Math.Cos(xx)*ans2);
			}
			return ans;
		}
		public static double J(int n, double x)
		{
			double ACC = 40.0;
			double BIGNO = 1.0e10;
			double BIGNI = 1.0e-10;

			int j,m;
			bool jsum;
			double ax,bj,bjm,bjp,sum,tox,ans;

            if (n < 2) try { throw new Science.Error(); }
                catch (Science.Error e)
                {
                    e.Write("Index n less than 2 in bessj");
                }

			ax=Math.Abs(x);
			if (ax == 0.0)
				return 0.0;
			else if (ax > (double) n) 
			{
				tox=2.0/ax;
				bjm=J0(ax);
				bj=J1(ax);
				for (j=1;j<n;j++) 
				{
					bjp=j*tox*bj-bjm;
					bjm=bj;
					bj=bjp;
				}
				ans=bj;
			} 
			else 
			{
				tox=2.0/ax;
				m=2*((n+(int) Math.Sqrt(ACC*n))/2);
				jsum=false;
				bjp=ans=sum=0.0;
				bj=1.0;
				for (j=m;j>0;j--) 
				{
					bjm=j*tox*bj-bjp;
					bjp=bj;
					bj=bjm;
					if (Math.Abs(bj) > BIGNO) 
					{
						bj *= BIGNI;
						bjp *= BIGNI;
						ans *= BIGNI;
						sum *= BIGNI;
					}
					if (jsum) sum += bj;
					jsum=!jsum;
					if (j == n) ans=bjp;
				}
				sum=2.0*sum-bj;
				ans /= sum;
			}
			return x < 0.0 && (n%2)== 1 ? -ans : ans;
		}
		public static double Y(int n, double x)
		{
			int j;
			double by,bym,byp,tox;

            if (n < 2) try { throw new Science.Error(); }
                catch (Science.Error e)
                {
                    e.Write("Index n less than 2 in bessy");
                }

			tox=2.0/x;
			by=Y1(x);
			bym=Y0(x);
			for (j=1;j<n;j++) 
			{
				byp=j*tox*by-bym;
				bym=by;
				by=byp;
			}
			return by;
		}
		public static double I0(double x)
		{
			double ax,ans;
			double y;

			if ((ax=Math.Abs(x)) < 3.75) 
			{
				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=(Math.Exp(ax)/Math.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;
		}
		public static double K0(double x)
		{
			double y,ans;

			if (x <= 2.0) 
			{
				y=x*x/4.0;
				ans=(-Math.Log(x/2.0)*I0(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=(Math.Exp(-x)/Math.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;
		}
		public static double I1(double x)
		{
			double ax,ans;
			double y;

			if ((ax=Math.Abs(x)) < 3.75) 
			{
				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 *= (Math.Exp(ax)/Math.Sqrt(ax));
			}
			return x < 0.0 ? -ans : ans;
		}
		public static double K1(double x)
		{
			double y,ans;

			if (x <= 2.0) 
			{
				y=x*x/4.0;
				ans=(Math.Log(x/2.0)*I1(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=(Math.Exp(-x)/Math.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;
		}
		public static double I(int n, double x)
		{
			double ACC = 40.0;
			double BIGNO = 1.0e10;
			double BIGNI = 1.0e-10;

			int j;
			double bi,bim,bip,tox,ans;

            if (n < 2) try { throw new Science.Error(); }
                catch (Science.Error e)
                {
                    e.Write("Index n less than 2 in bessi");
                }

			if (x == 0.0)
				return 0.0;
			else 
			{
				tox=2.0/Math.Abs(x);
				bip=ans=0.0;
				bi=1.0;
				for (j=2*(n+(int) Math.Sqrt(ACC*n));j>0;j--) 
				{
					bim=bip+j*tox*bi;
					bip=bi;
					bi=bim;
					if (Math.Abs(bi) > BIGNO) 
					{
						ans *= BIGNI;
						bi *= BIGNI;
						bip *= BIGNI;
					}
					if (j == n) ans=bip;
				}
				ans *= I0(x)/bi;
				return x < 0.0 && (n%2) == 1 ? -ans : ans;
			}
		}
		public static double K(int n, double x)
		{
			int j;
			double bk,bkm,bkp,tox;

            if (n < 2) try { throw new Science.Error(); }
                catch (Science.Error e)
                {
                    e.Write("Index n less than 2 in bessk");
                }

			tox=2.0/x;
			bkm=K0(x);
			bk=K1(x);
			for (j=1;j<n;j++) 
			{
				bkp=bkm+j*tox*bk;
				bkm=bk;
				bk=bkp;
			}
			return bk;
		}
		
		public static double J(double nu, double x)
		{
			return bessjy(x,nu,1);
		}
		public static double Y(double nu, double x)
		{
			return bessjy(x,nu,2);
		}
		public static double DerivativeOfJ(double nu, double x)
		{
			return bessjy(x,nu,3);
		}
		public static double DerivativeOfY(double nu, double x)
		{
			return bessjy(x,nu,4);
		}

		private static double bessjy(double x, double xnu, int returnwhat)
		{
			double EPS = 1.0e-16;
			double FPMIN = 1.0e-30;
			int MAXIT = 10000;
			double XMIN = 2.0;
			double PI = 3.141592653589793;

			double rj, ry, rjp, ryp;

			int i,isign,l,nl;
			double a,b,br,bi,c,cr,ci,d,del,del1,den,di,dlr,dli,dr,e,f,fact,fact2,
				fact3,ff,gam,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) try { throw new Science.Error(); }
                catch (Science.Error e1)
                {
                    e1.Write("bad arguments in bessjy");
                }

			nl=(x < XMIN ? (int)(xnu+0.5) : Math.Max(0,(int)(xnu-x+1.5)));
			xmu=xnu-nl;
			xmu2=xmu*xmu;
			xi=1.0/x;
			xi2=2.0*xi;
			w=xi2/PI;
			isign=1;
			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=b-d;
				if (Math.Abs(d) < FPMIN) d=FPMIN;
				c=b-1.0/c;
				if (Math.Abs(c) < FPMIN) c=FPMIN;
				d=1.0/d;
				del=c*d;
				h=del*h;
				if (d < 0.0) isign = -isign;
				if (Math.Abs(del-1.0) < EPS) break;
			}

            if (i > MAXIT) try { throw new Science.Error(); }
                catch (Science.Error e2)
                {
                    e2.Write("x too large in bessjy; try asymptotic expansion");
                }

			rjl=isign*FPMIN;
			rjpl=h*rjl;
			rjl1=rjl;
			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) 
			{
				x2=0.5*x;
				pimu=PI*xmu;
				fact = (Math.Abs(pimu) < EPS ? 1.0 : pimu/Math.Sin(pimu));
				d = -Math.Log(x2);
				e=xmu*d;
				fact2 = (Math.Abs(e) < EPS ? 1.0 : Math.Sinh(e)/e);


				int NUSE1 = 7;
				int NUSE2 = 8;
				double gam1, gam2, gampl, gammi;

				double xx;
				double[] c1 = {	-1.142022680371172e0,6.516511267076e-3,
								  3.08709017308e-4,-3.470626964e-6,6.943764e-9,
								  3.6780e-11,-1.36e-13};
				double[] 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;


				ff=2.0/PI*fact*(gam1*Math.Cosh(e)+gam2*fact2*d);
				e=Math.Exp(e);
				p=e/(gampl*PI);
				q=1.0/(e*PI*gammi);
				pimu2=0.5*pimu;
				fact3 = (Math.Abs(pimu2) < EPS ? 1.0 : Math.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 (Math.Abs(del) < (1.0+Math.Abs(sum))*EPS) break;
				}

                if (i > MAXIT) try { throw new Science.Error(); }
                    catch (Science.Error e3)
                    {
                        e3.Write("bessy series failed to converge");
                    }
            	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 (Math.Abs(dr)+Math.Abs(di) < FPMIN) dr=FPMIN;
					fact=a/(cr*cr+ci*ci);
					cr=br+cr*fact;
					ci=bi-ci*fact;
					if (Math.Abs(cr)+Math.Abs(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 (Math.Abs(dlr-1.0)+Math.Abs(dli) < EPS) break;
				}
                if (i > MAXIT) try { throw new Science.Error(); }
                    catch (Science.Error e4)
                    {
                        e4.Write("cf2 failed in bessjy");
                    }

                gam=(p-f)/q;
				rjmu=Math.Sqrt(w/((p-f)*gam+q));
				rjmu=rjmu*Math.Sign(rjl);
				rymu=rjmu*gam;
				rymup=rymu*(p+q/gam);
				ry1=xmu*xi*rymu-rymup;
			}
			fact=rjmu/rjl;
			rj=rjl1*fact;
			rjp=rjp1*fact;
			for (i=1;i<=nl;i++) 
			{
				rytemp=(xmu+i)*xi2*ry1-rymu;
				rymu=ry1;
				ry1=rytemp;
			}
			ry=rymu;
			ryp=xnu*xi*rymu-ry1;

			if(returnwhat == 1) return rj;
			else if(returnwhat == 2) return ry;
			else if(returnwhat == 3) return rjp;
			else return ryp;
		} 
		private static double Chebev(double a, double b, double[] c, int m, double x)
		{
			double d=0.0,dd=0.0,sv,y,y2;
			int j;

            if ((x - a) * (x - b) > 0.0) try { throw new Science.Error(); }
                catch (Science.Error e5)
                {
                    e5.Write("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];
		}

		public static double I(double nu, double x)
		{
			return bessik(x,nu,1);
		}
		public static double K(double nu, double x)
		{
			return bessik(x,nu,2);
		}
		public static double DerivativeOfI(double nu, double x)
		{
			return bessik(x,nu,3);
		}
		public static double DerivativeOfK(double nu, double x)
		{
			return bessik(x,nu,4);
		}
		
    	private static double bessik(double x, double xnu, int returnwhat)
		{
			double EPS = 1.0e-16;
			double FPMIN = 1.0e-30;
			int MAXIT = 10000;
			double XMIN = 2.0;
			double PI = 3.141592653589793;

			double ri, rk, rip, rkp;

			int i,l,nl;
			double a,a1,b,c,d,del,del1,delh,dels,e,f,fact,fact2,ff,
				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) try { throw new Science.Error(); }
                catch (Science.Error e6)
                {
                    e6.Write("bad arguments in bessik");
                }

            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 (Math.Abs(del-1.0) < EPS) break;
			}

            if (x <= 0.0 || xnu < 0.0) try { throw new Science.Error(); }
                catch (Science.Error e7)
                {
                    e7.Write("x too large in bessik; try asymptotic expansion");
                }

			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 = (Math.Abs(pimu) < EPS ? 1.0 : pimu/Math.Sin(pimu));
				d = -Math.Log(x2);
				e=xmu*d;
				fact2 = (Math.Abs(e) < EPS ? 1.0 : Math.Sinh(e)/e);

				int NUSE1 = 7;
				int NUSE2 = 8;
				double gam1, gam2, gampl, gammi;

				double xx;
				double[] c1 = {	-1.142022680371172e0,6.516511267076e-3,
								  3.08709017308e-4,-3.470626964e-6,6.943764e-9,
								  3.6780e-11,-1.36e-13};
				double[] 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;


				ff=fact*(gam1*Math.Cosh(e)+gam2*fact2*d);
				sum=ff;
				e=Math.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 (Math.Abs(del) < Math.Abs(sum)*EPS) break;
				}

                if (i > MAXIT) try { throw new Science.Error(); }
                    catch (Science.Error e8)
                    {
                        e8.Write("bessk series failed to converge");
                    }

				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 (Math.Abs(dels/s) < EPS) break;
				}

                if (i > MAXIT) try { throw new Science.Error(); }
                    catch (Science.Error e9)
                    {
                        e9.Write("bessik: failure to converge in cf2");
                    }

				h=a1*h;
				rkmu=Math.Sqrt(PI/(2.0*x))*Math.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;

			if(returnwhat == 1) return ri;
			else if(returnwhat == 2) return rk;
			else if(returnwhat == 3) return rip;
			else return rkp;
		}

		public static double Ai(double x)
		{
			return airy(x,1);
		}
		public static double Bi(double x)
		{
			return airy(x,2);
		}
		public static double DerivativeOfAi(double x)
		{
			return airy(x,3);
		}
		public static double DerivativeOfBi(double x)
		{
			return airy(x,4);
		}

     	private static double airy(double x, int what)
		{
			double THIRD = 1.0/3.0;
			double TWOTHR = 2.0/3.0;
			double ONOVRT = 0.57735027;

			double ai, bi, aip, bip;
			double absx,rootx,z;

			absx=Math.Abs(x);
			rootx=Math.Sqrt(absx);
			z=TWOTHR*absx*rootx;
			if (x > 0.0) 
			{
				ai=rootx*ONOVRT*bessik(z,THIRD,2)/Math.PI;
				bi=rootx*(bessik(z,THIRD,2)/Math.PI+2.0*ONOVRT*bessik(z,THIRD,1));
				aip = -x*ONOVRT*bessik(z,TWOTHR,2)/Math.PI;
				bip=x*(bessik(z,TWOTHR,2)/Math.PI+2.0*ONOVRT*bessik(z,TWOTHR,1));
			} 
			else if (x < 0.0) 
			{
				ai=0.5*rootx*(bessjy(z,THIRD,1)-ONOVRT*bessjy(z,THIRD,2));
				bi = -0.5*rootx*(bessjy(z,THIRD,2)+ONOVRT*bessjy(z,THIRD,1));
				aip=0.5*absx*(ONOVRT*bessjy(z,TWOTHR,2)+bessjy(z,TWOTHR,1));
				bip=0.5*absx*(ONOVRT*bessjy(z,TWOTHR,1)-bessjy(z,TWOTHR,2));
			} 
			else 
			{
				ai=0.35502805;
				bi=ai/ONOVRT;
				aip = -0.25881940;
				bip = -aip/ONOVRT;
			}
			if(what == 1) return ai;
			else if(what == 2) return bi;
			else if(what == 3) return aip;
			else return bip;
		}
	
		public static double j(int n, double x)
		{
			return sphbes(n,x,1);
		}
		public static double y(int n, double x)
		{
			return sphbes(n,x,2);
		}
		public static double DerivativeOfj(int n, double x)
		{
			return sphbes(n,x,3);
		}
		public static double DerivativeOfy(int n, double x)
		{
			return sphbes(n,x,4);
		}

    	private static double sphbes(int n, double x, int what)
		{
			double RTPIO2 = 1.2533141;
			double sj, sy, sjp, syp;
			double factor,order;


            if (n < 0 || x <= 0.0) try { throw new Science.Error(); }
                catch (Science.Error e10)
                {
                    e10.Write("bad arguments in sphbes");
                }

			order=n+0.5;
			factor=RTPIO2/Math.Sqrt(x);
			sj=factor*bessjy(x,order,1);
			sy=factor*bessjy(x,order,2);
			sjp=factor*bessjy(x,order,3)-sj/(2.0*x);
			syp=factor*bessjy(x,order,4)-sy/(2.0*x);

			if(what == 1) return sj;
			else if(what == 2) return sy;
			else if(what == 3) return sjp;
			else return syp;
		}
	}
}