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:

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.


**********SPHERICAL BESSEL FUNCTIONS************

/*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 */
		out = (2*n-1)*bn(z,n-1)/z-bn(z,n-2);
	return out;
	double z;
	int n;

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):

Spherical Bessel Functions
Jn(z) from n=0 to 5

YouTube Tutorial:

[wpedon id="7041" align="center"]

Leave a Reply

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

This site uses Akismet to reduce spam. Learn how your comment data is processed.