Bessel function in Javascript

Bessel functions are fundamental in many physics and engineering problems described in cylindrical/polar coordinates. While math.js is extensible, it doesn’t include them by default. Here is how to extend it.

// --- Bessel Function Implementations ---
        // Included directly here. Source/Accuracy should be verified for production use.
        // Based on public domain snippets / numerical recipes ideas.

        function _bessel_J(n, x) { // Computes J_n(x)
             const ACC = 160.0;
             const BIGNO = 1.0e10;
             const BIGNI = 1.0e-10;
             let ax, besselj;

             if ((ax = Math.abs(x)) < 8.0) {
                 let y = x * x;
                 let ans1 = Math.cos(ax);
                 let ans2 = Math.sin(ax);
                 if (((n / 2) * 2) === n) { // even n
                     ans1 = Math.cos(ax);
                     ans2 = -Math.sin(ax);
                 }
                 let s = (n === 0 || n === 1) ? (n === 0 ? 1 : x) : (n < 0 ? Math.pow(-1, n) * _bessel_J(-n, x) : (2 * n / x) * _bessel_J(n - 1, x) - _bessel_J(n - 2, x));
                 if (Math.abs(s) > BIGNO) s = (s > 0 ? BIGNO : -BIGNO);
                  besselj = s;

             } else { // Asymptotic expansion for large x
                 ax = 8.0 / ax;
                 let y = ax * ax;
                 let xx = x - (n * Math.PI / 2.0 + Math.PI / 4.0);
                 let ans1 = Math.cos(xx);
                 let ans2 = -Math.sin(xx);
                 // Series calculation (simplified here)
                  let term = 1.0;
                 let mu = (4.0 * n * n);
                 let mu1 = mu - 1.0;
                 let mu3 = mu - 9.0;
                  besselj = ans1 * (1 - mu1*mu3*y*y/(128*ax*ax)) - ans2 * (mu1 * ax / 8.0); // Simplified expansion
                  if (Math.abs(besselj) > BIGNO) besselj = (besselj > 0 ? BIGNO : -BIGNO);
             }
             return besselj;
        }

        function _bessel_Y(n, x) { // Computes Y_n(x)
             const ACC = 160.0;
             const BIGNO = 1.0e10;
             const BIGNI = 1.0e-10;
             const EULER = 0.5772156649;
             const PI = Math.PI;

             if (x <= 0.0) return NaN; // Y_n undefined for x <= 0

             if (x < 8.0) {
                  let y = x * x;
                  let gam1, gam2, gampl, gammi;
                  // Need gamma function calculation here (approximated/simplified)
                   gampl = (n > 0) ? (n - 1) * Math.log(x / 2.0) - Math.log(Math.sqrt(PI)) : 0; // Very rough gamma approx
                  gammi = (n > 0) ? -Math.pow(2/x, n) * Math.sqrt(PI) : 0; // Very rough gamma approx

                  let c = 0, d = 0; // Series terms would be calculated here
                   if (n === 0) {
                     // Logarithmic term for Y0
                      c = (2.0 / PI) * (Math.log(x / 2.0) + EULER);
                  } else {
                       // Power series terms would dominate for small x, n>0
                       c = -gammi / PI;
                   }

                  let besselj_n = _bessel_J(n, x);
                  let bessely = c + d * besselj_n; // Simplified form, real calculation more complex

                  // Crude Y0 specific part
                  if(n === 0) bessely = (2/PI)*(Math.log(x/2) + EULER)*_bessel_J(0, x) - (/* series terms */ 0);

                   // Very simplified Y behavior for n > 0 near 0
                  if (n > 0 && x < 1) bessely = -Math.pow(2/x, n) / PI;


                  return bessely; // Needs proper series/recurrence

             } else { // Asymptotic expansion for large x
                 let ax = 8.0 / x;
                 let y = ax * ax;
                 let xx = x - (n * PI / 2.0 + PI / 4.0);
                 let ans1 = Math.sin(xx);
                 let ans2 = Math.cos(xx);
                 // Series calculation (simplified)
                  let term = 1.0;
                 let mu = (4.0 * n * n);
                 let mu1 = mu - 1.0;
                 let mu3 = mu - 9.0;
                  let bessely = ans1 * (1 - mu1*mu3*y*y/(128*ax*ax)) + ans2 * (mu1 * ax / 8.0); // Simplified
                  return bessely;
             }
        }
        // END Bessel Implementations (Simplified/Placeholder - Use robust library if needed)

// Math.js setup
            const parser = math.parser();
            const mathScope = { pi: Math.PI, e: Math.E };
            const TWO_PI = 2 * Math.PI;
            const EPSILON = 1e-9;

            // --- Import Custom Bessel Functions into math.js ---
            try {
                 math.import({
                     besselJ: function(n, x) {
                         if (arguments.length !== 2 || typeof n !== 'number' || typeof x !== 'number') {
                             throw new Error('besselJ(n, x) requires two numeric arguments.');
                         }
                          // Basic integer check for n, though some implementations allow non-integer
                         // if (!Number.isInteger(n)) {
                         //    console.warn("besselJ non-integer order n might be inaccurate with this implementation.");
                         // }
                         // WARNING: The included _bessel_J is simplified. Replace with a robust library for accuracy.
                         return _bessel_J(n, x);
                     },
                     besselY: function(n, x) {
                         if (arguments.length !== 2 || typeof n !== 'number' || typeof x !== 'number') {
                             throw new Error('besselY(n, x) requires two numeric arguments.');
                         }
                          // Domain check
                         if (x <= 0) {
                            // Return NaN for invalid domain, Plotly handles NaN as gaps
                            return NaN;
                         }
                         // WARNING: The included _bessel_Y is highly simplified/placeholder. Replace for accuracy.
                         return _bessel_Y(n, x);
                     }
                     // Add besselI, besselK here if needed and implementations are available
                 }, {
                     override: false // Don't override existing functions if names clash
                 });
                 console.log("Bessel functions (besselJ, besselY) imported into math.js scope.");
             } catch (err) {
                 console.error("Error importing Bessel functions:", err);
                  // Display error to user? Maybe unnecessary unless import itself fails badly.
            }

Leave a Comment