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:
Ph.D. researcher at Friedrich-Schiller University Jena, Germany. I’m a physicist specializing in computational material science. I write efficient codes for simulating light-matter interactions at atomic scales. I like to develop Physics, DFT, and Machine Learning related apps and software from time to time. Can code in most of the popular languages. I like to share my knowledge in Physics and applications using this Blog and a YouTube channel.