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

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

YouTube Tutorial:

Get Article as PDF

Manas Sharma

I’m a physicist specializing in computational material science with a PhD in Physics from Friedrich-Schiller University Jena, Germany. 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.

manas.bragitoff.com/