diff --git a/src/gpuAAK.cu b/src/gpuAAK.cu
index ca55a396..29ff3fbf 100644
--- a/src/gpuAAK.cu
+++ b/src/gpuAAK.cu
@@ -19,6 +19,7 @@ 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)))));
@@ -26,9 +27,10 @@ static double bessel_j0( double x ) {
          +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,12 +42,12 @@ 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))))));
@@ -53,9 +55,10 @@ static double bessel_j1( double x ) {
          +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;
 }