Spherical Bessel Function Integral Representation – C PROGRAM

Spherical Bessel Functions of the first kind of order n can be written in the integral form as follows:
j_n(z)=\frac{z^n}{2^{n+1}n!}\int_0^\pi \cos (z\cos \theta)\sin^{2n+1}\theta d\theta

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:

PhD researcher at Friedrich-Schiller University Jena, Germany. I'm a physicist specializing in theoretical, computational and experimental condensed matter physics. I like to develop Physics related apps and softwares from time to time. Can code in most of the popular languages. Like to share my knowledge in Physics and applications using this Blog and a YouTube channel.



Leave a Reply

Your email address will not be published. Required fields are marked *