Spherical Bessel Functions of the first kind of order n can be written in the integral form as follows:

Using the Simpson’s Method for Numerical Integration, we can write a program that can calculate the spherical Bessel functions.

### PROGRAM:

/******************************************* *SPHERICAL BESSEL FUNCTIONS USING INTEGRALS* *******************************************/ #include<stdio.h> #include<math.h> //function to calculate the factorial double factorial(int n){ int i; double fact=1; for(i=n;i>=1;i--){ fact=fact*i; } return fact; } /*Function definition to perform integration by Simpson's 1/3rd Rule */ double simpsons(double z,int n,double f(double theta,double z,int n),double a,double b,int intervals){ double h,integral,theta,sum=0; int i; h=fabs(b-a)/intervals; for(i=1;i<intervals;i++){ theta=a+i*h; if(i%2==0){ sum=sum+2*f(theta,z,n); } else{ sum=sum+4*f(theta,z,n); } } integral=(h/3)*(f(theta,z,n)+f(theta,z,n)+sum); return integral; } //function that contains the integral part of the spherical bessel function double f(double theta,double z, int n){ return cos(z*cos(theta))*pow(sin(theta),2*n+1); } //function that returns the value of the spherical bessel function for a given x and n double Jn(double z, int n){ return pow(z,n)/(pow(2,n+1)*factorial(n))*simpsons(z,n,f,0,M_PI,1000); } main(){ double z; FILE *fp=NULL; fp=fopen("sphericalBesselIntegral.txt","w"); for(z=0;z<=10;z=z+0.01){ fprintf(fp,"%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",z,Jn(z,0),Jn(z,1),Jn(z,2),Jn(z,3),Jn(z,4),Jn(z,5)); } } //GNUPLOT COMMAND: //plot 'sphericalBesselIntegral.txt' u 1:2 w l t "j0(z)", '' u 1:3 w l t "j1(z)", '' u 1:4 w l t "j2(z)", '' u 1:5 w l t "j3(z)", '' u 1:6 w l t "j4(z)", '' u 1:7 w l t "j5(z)"

### OUTPUT:

The data-points are stored in a file called `sphericalBesselIntegral.txt`

.

They can be plotted on any software.

Copy-pasting the data to Excel and then plotting it gives the following result:

One can plot the file using **Gnuplot **command:

`plot 'sphericalBesselIntegral.txt' u 1:2 w l t "j0(z)", '' u 1:3 w l t "j1(z)", '' u 1:4 w l t "j2(z)", '' u 1:5 w l t "j3(z)", '' u 1:6 w l t "j4(z)", '' u 1:7 w l t "j5(z)"`

**Output**: