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:

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 *