Spherical Bessel Functions – C Program

In this post I will show you how to calculate and plot the Spehrical Bessel Functions ( $j_n(z)$) of the first kind using C and Gnuplot.

We will use the following information: $j_0(z)=\frac{\sin z}{z}$ $j_1(z)=\frac{\sin z}{z^2}-\frac{\cos z}{z}$
and the recurrence relation: $j_{n-1}(z)+j_{n+1}(z)=(2n+1)\frac{j_n(z)}{z}$

We will create a program that calculates the values of the Bessel Function at various z values and for different n and store these values in a txt file. Then just plot it using Gnuplot.

We will create two functions called ‘b0’ and ‘b1’, that contain the definition of $j_0(z) , j_1(z)$ respectively.
Then we will create a function ‘bn’ that will use the first two functions and recursion to find the value of Bessel function for different z,n.
NOTE: I am using a slightly modified form of the recurrence relation. To get the form I am using, just replace n by n-1.

C PROGRAM:

/***********************************************
**********SPHERICAL BESSEL FUNCTIONS************
***********************************************/
#include<stdio.h>
#include<math.h>

/*Define j0(z) */
double b0(double z){
return sin(z)/z;
}

/*Define j1(z) */
double b1(double z){
return sin(z)/(z*z)-cos(z)/z;
}

/*Define jn(z) */
double bn(double z,int n){
double out;
if (n==0){
out = b0(z);
}
else if(n==1){
out = b1(z);
}
/*using recurrence relation */
else{
out = (2*n-1)*bn(z,n-1)/z-bn(z,n-2);
}
return out;
}
main(){
double z;
int n;
FILE *fp=NULL;
fp=fopen("bessel.txt","w");
for(z=0.01;z<=20;z=z+0.01){
//fprintf(fp,"%lf\t%lf\n",z,bn(z,3));
fprintf(fp,"%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",z,bn(z,0),bn(z,1),bn(z,2),bn(z,3),bn(z,4),bn(z,5));
}

}

When you run the above C, it will generate a file called ‘bessel.txt’ which would contain 7 columns of data-points.
The first column contains the ‘z’ values and the rest of them are for $j_0(z), j_1(z),....$

These can be easily plotted using Gnuplot by using the following commands:

Gnuplot Command:

->set xlabel "z"
->plot 'bessel.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 (Gnuplot): 