/* This file contains the functions to calculate Bessel functions. */

#include <bessel.h>

/* bessJ0 calculates the bessel function J_0 */

double bessJ0(double x)
{
	double ax,z;
	double xx,y,ans,ans1,ans2;

	if ((ax=fabs(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=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
	}
	return ans;
}


/* neumannN0 calculates the Neumann function N_0 */

double neumannN0(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*bessJ0(x)*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=sqrt(0.636619772/x)*(sin(xx)*ans1+z*cos(xx)*ans2);
	}
	return ans;
}


/* Hankel1_0 calculates the Hankel function H(1)_0 */

dcomplex Hankel1_0(double x)
{
     return Complex(bessJ0(x), neumannN0(x));
}

/*************************************************************************/

/* bessel_spher calculates spherical bessel function */

double bessel_spher(int l, double x)
{
     double s;
  
     switch(l)
     {
     case 0 : 
            { 
            s = sin(x) / x;
            break;
            }
     case 1 : 
            { 
            s = ( sin(x) / pow(x,2.) - cos(x) / x );
            break;
            }
     case 2 : 
            { 
            s = ( sin(x) * ( 3.0 / pow(x, 3.) - 1 / x ) - 
                  cos(x) * ( 3.0 / pow(x, 2.)) );
            break;
            }
     }

     return s;      
}