document bessel function

BlackHolePerturbationToolkit · Dec 4, 2024 · 0952973 · 0952973
1 parent f1f0250
commit 0952973
Showing 1 changed file with 12 additions and 3 deletions.
diff --git a/src/gpuAAK.cu b/src/gpuAAK.cu
@@ -19,16 +19,18 @@ static double bessel_j0( double x ) {
    double xx,y,ans,ans1,ans2;
 
    if ((ax=fabs(x)) < 8.0) {
+      // Padé approximation (composed of numerator ans1, denominator ans2)
       y=x*x;
       ans1=57568490574.0+y*(-13362590354.0+y*(651619640.7
          +y*(-11214424.18+y*(77392.33017+y*(-184.9052456)))));
       ans2=57568490411.0+y*(1029532985.0+y*(9494680.718
          +y*(59272.64853+y*(267.8532712+y*1.0))));
       ans=ans1/ans2;
    } else {
+      // Asymptotic expansion that generally improves as x grows
       z=8.0/ax;
       y=z*z;
-      xx=ax-0.785398164;
+      xx=ax-0.785398164;  // minus pi/4
       ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
          +y*(-0.2073370639e-5+y*0.2093887211e-6)));
       ans2 = -0.1562499995e-1+y*(0.1430488765e-3
@@ -40,22 +42,23 @@ static double bessel_j0( double x ) {
 }
 
 static double bessel_j1( double x ) {
-   // Numerical Recipes, pp. 274-280.
 
    double ax,z;
    double xx,y,ans,ans1,ans2;
 
    if ((ax=fabs(x)) < 8.0) {
+      // Padé approximation (composed of numerator ans1, denominator ans2)
       y=x*x;
       ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1
          +y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
       ans2=144725228442.0+y*(2300535178.0+y*(18583304.74
          +y*(99447.43394+y*(376.9991397+y*1.0))));
       ans=ans1/ans2;
    } else {
+      // Asymptotic expansion that generally improves as x grows
       z=8.0/ax;
       y=z*z;
-      xx=ax-2.356194491;
+      xx=ax-2.356194491;  // minus 3pi/4
       ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
          +y*(0.2457520174e-5+y*(-0.240337019e-6))));
       ans2=0.04687499995+y*(-0.2002690873e-3
@@ -69,11 +72,13 @@ static double bessel_j1( double x ) {
 
 double bessel_jn( int n, double x ) {
    // Numerical Recipes, pp. 274-280.
+   // Constructs J_n for arbitrary n with a recurrence relation
 
    int    j, jsum, m;
    double ax, bj, bjm, bjp, sum, tox, ans;
 
    ax=fabs(x);
+   // handle special cases
    if (n == 0)
       return( bessel_j0(ax) );
    if (n == 1)
@@ -82,6 +87,7 @@ double bessel_jn( int n, double x ) {
    if (ax == 0.0)
       return 0.0;
    else if (ax > (double) n) {
+      // forward recurrence (assumed stable above n)
       tox=2.0/ax;
       bjm=bessel_j0(ax);
       bj=bessel_j1(ax);
@@ -92,6 +98,7 @@ double bessel_jn( int n, double x ) {
       }
       ans=bj;
    } else {
+      // backwards recurrence for small x
       tox=2.0/ax;
       m=2*((n+(int) sqrt(ACC*n))/2);
       jsum=0;
@@ -114,6 +121,8 @@ double bessel_jn( int n, double x ) {
       sum=2.0*sum-bj;
       ans /= sum;
    }
+
+   // adjust sign for negative x if n odd
    return  x < 0.0 && n%2 == 1 ? -ans : ans;
 }