From 45c1c7b614b20e916abf31f949cb30198343175b Mon Sep 17 00:00:00 2001 From: SheetJS Date: Sat, 14 Dec 2013 02:36:05 -0500 Subject: [PATCH] Initial commit --- .gitignore | 1 + .npmignore | 2 + .travis.yml | 6 + LICENSE | 9 ++ Makefile | 7 ++ README.md | 36 ++++++ bessel.js | 209 +++++++++++++++++++++++++++++++++ bessel.md | 317 +++++++++++++++++++++++++++++++++++++++++++++++++++ package.json | 19 +++ test.html | 4 + test.js | 2 + 11 files changed, 612 insertions(+) create mode 100644 .gitignore create mode 100644 .npmignore create mode 100644 .travis.yml create mode 100644 LICENSE create mode 100644 Makefile create mode 100644 README.md create mode 100644 bessel.js create mode 100644 bessel.md create mode 100644 package.json create mode 100644 test.html create mode 100644 test.js diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..3c3629e --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +node_modules diff --git a/.npmignore b/.npmignore new file mode 100644 index 0000000..5637c09 --- /dev/null +++ b/.npmignore @@ -0,0 +1,2 @@ +Makefile +test.html diff --git a/.travis.yml b/.travis.yml new file mode 100644 index 0000000..a1f55b5 --- /dev/null +++ b/.travis.yml @@ -0,0 +1,6 @@ +language: node_js +node_js: + - "0.10" + - "0.8" +before_install: + - "npm install -g mocha" diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..93d22d9 --- /dev/null +++ b/LICENSE @@ -0,0 +1,9 @@ +Copyright (C) 2013 SheetJS + +The MIT License (MIT) + +Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..bc8c57c --- /dev/null +++ b/Makefile @@ -0,0 +1,7 @@ +LIBRARY=bessel + +$(LIBRARY).js: bessel.md + node_modules/.bin/voc $^ > $@ + +test mocha: + mocha -R spec diff --git a/README.md b/README.md new file mode 100644 index 0000000..3ce6439 --- /dev/null +++ b/README.md @@ -0,0 +1,36 @@ +# Bessel.JS + +Pure-JS implementation of the Bessel functions (J,Y,I,K), for node and browser + +The standard notation is used here: + + - J is the Bessel function of the first kind + - Y is the Bessel function of the second kind + - I is the modified Bessel function of the first kind + - K is the modified Bessel function of the first kind + +# Usage + +The functions `besselj`, `bessely`, `besseli`, `besselk` are exposed when you include +the script `bessel.js`: + +```html> + + +``` + +See `test.html` for an example + +In node, those four functions are exported: + +```js> +var besselj01 = require('bessel').besselj(0,1); +``` + +Each function follows Excel semantics `(value, function-order)`. For example, + +```js> +bessel.besselj(1.5, 1) +``` + +is the value of the bessel function J1 at the point x=1.5 diff --git a/bessel.js b/bessel.js new file mode 100644 index 0000000..88c0eef --- /dev/null +++ b/bessel.js @@ -0,0 +1,209 @@ +var M = Math; +Array.prototype.horner = function(v) { return this.reduce(function(z,w){return v * z + w;},0); }; +function _bessel_iter(x, n, f0, f1, sign) { + if(!sign) sign = -1; + var tdx = 2 / x, f2; + if(n === 0) return f0; + if(n === 1) return f1; + for(var o = 1; o != n; ++o) { + f2 = f1 * o * tdx + sign * f0; + f0 = f1; f1 = f2; + } + return f1; +} +function _bessel_wrap(bessel0, bessel1, name, nonzero, sign) { + return function bessel(x,n) { + if(n === 0) return bessel0(x); + if(n === 1) return bessel1(x); + if(n < 0) throw name + ': Order (' + n + ') must be nonnegative'; + if(nonzero == 1 && x === 0) throw name + ': Undefined when x == 0'; + if(nonzero == 2 && x <= 0) throw name + ': Undefined when x <= 0'; + var b0 = bessel0(x), b1 = bessel1(x); + return _bessel_iter(x, n, b0, b1, sign); + }; +} +var besselj = (function() { + var b0_a1a = [57568490574.0,-13362590354.0,651619640.7,-11214424.18,77392.33017,-184.9052456].reverse(); + var b0_a2a = [57568490411.0,1029532985.0,9494680.718,59272.64853,267.8532712,1.0].reverse(); + var b0_a1b = [1.0, -0.1098628627e-2, 0.2734510407e-4, -0.2073370639e-5, 0.2093887211e-6].reverse(); + var b0_a2b = [-0.1562499995e-1, 0.1430488765e-3, -0.6911147651e-5, 0.7621095161e-6, -0.934935152e-7].reverse(); + var W = 0.636619772; // 2 / Math.PI + + function bessel0(x) { + var a, a1, a2, y = x * x, xx = M.abs(x) - 0.785398164; + if(M.abs(x) < 8) { + a1 = b0_a1a.horner(y); + a2 = b0_a2a.horner(y); + a = a1/a2; + } + else { + y = 64 / y; + a1 = b0_a1b.horner(y); + a2 = b0_a2b.horner(y); + a = M.sqrt(W/M.abs(x))*(M.cos(xx)*a1-M.sin(xx)*a2*8/M.abs(x)); + } + return a; + } + var b1_a1a = [72362614232.0,-7895059235.0,242396853.1,-2972611.439, 15704.48260, -30.16036606].reverse(); + var b1_a2a = [144725228442.0, 2300535178.0, 18583304.74, 99447.43394, 376.9991397, 1.0].reverse(); + var b1_a1b = [1.0, 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, -0.240337019e-6].reverse(); + var b1_a2b = [0.04687499995, -0.2002690873e-3, 0.8449199096e-5, -0.88228987e-6, 0.105787412e-6].reverse(); + function bessel1(x) { + var a, a1, a2, y = x*x, xx = M.abs(x) - 2.356194491; + if(Math.abs(x)< 8) { + a1 = x*b1_a1a.horner(y); + a2 = b1_a2a.horner(y); + a = a1 / a2; + } else { + y = 64 / y; + a1=b1_a1b.horner(y); + a2=b1_a2b.horner(y); + a=M.sqrt(W/M.abs(x))*(M.cos(xx)*a1-M.sin(xx)*a2*8/M.abs(x)); + if(x < 0) a = -a; + } + return a; + } + return function besselj(x, n) { + n = Math.round(n); + if(n === 0) return bessel0(M.abs(x)); + if(n === 1) return bessel1(M.abs(x)); + if(n < 0) throw 'BESSELJ: Order (' + n + ') must be nonnegative'; + if(M.abs(x) === 0) return 0; + + var ret, j, tox = 2 / M.abs(x), m, jsum, sum, bjp, bj, bjm; + if(M.abs(x) > n) { + ret = _bessel_iter(x, n, bessel0(M.abs(x)), bessel1(M.abs(x)),-1); + } else { + m=2*M.floor((n+M.floor(M.sqrt(40*n)))/2); + jsum=0; + bjp=ret=sum=0.0; + bj=1.0; + for (j=m;j>0;j--) { + bjm=j*tox*bj-bjp; + bjp=bj; + bj=bjm; + if (M.abs(bj) > 1E10) { + bj *= 1E-10; + bjp *= 1E-10; + ret *= 1E-10; + sum *= 1E-10; + } + if (jsum) sum += bj; + jsum=!jsum; + if (j == n) ret=bjp; + } + sum=2.0*sum-bj; + ret /= sum; + } + return x < 0 && (n%2) ? -ret : ret; + }; +})(); +var bessely = (function() { + var b0_a1a = [-2957821389.0, 7062834065.0, -512359803.6, 10879881.29, -86327.92757, 228.4622733].reverse(); + var b0_a2a = [40076544269.0, 745249964.8, 7189466.438, 47447.26470, 226.1030244, 1.0].reverse(); + var b0_a1b = [1.0, -0.1098628627e-2, 0.2734510407e-4, -0.2073370639e-5, 0.2093887211e-6].reverse(); + var b0_a2b = [-0.1562499995e-1, 0.1430488765e-3, -0.6911147651e-5, 0.7621095161e-6, -0.934945152e-7].reverse(); + + var W = 0.636619772; + function bessel0(x) { + var a, a1, a2, y = x * x, xx = x - 0.785398164; + if(x < 8) { + a1 = b0_a1a.horner(y); + a2 = b0_a2a.horner(y); + a = a1/a2 + W * besselj(x,0) * M.log(x); + } else { + y = 64 / y; + a1 = b0_a1b.horner(y); + a2 = b0_a2b.horner(y); + a = M.sqrt(W/x)*(M.sin(xx)*a1+M.cos(xx)*a2*8/x); + } + return a; + } + + var b1_a1a = [-0.4900604943e13, 0.1275274390e13, -0.5153438139e11, 0.7349264551e9, -0.4237922726e7, 0.8511937935e4].reverse(); + var b1_a2a = [0.2499580570e14, 0.4244419664e12, 0.3733650367e10, 0.2245904002e8, 0.1020426050e6, 0.3549632885e3, 1].reverse(); + var b1_a1b = [1.0, 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, -0.240337019e-6].reverse(); + var b1_a2b = [0.04687499995, -0.2002690873e-3, 0.8449199096e-5, -0.88228987e-6, 0.105787412e-6].reverse(); + function bessel1(x) { + var a, a1, a2, y = x*x, xx = x - 2.356194491; + if(x < 8) { + a1 = x*b1_a1a.horner(y); + a2 = b1_a2a.horner(y); + a = a1/a2 + W * (besselj(x,1) * M.log(x) - 1 / x); + } else { + y = 64 / y; + a1=b1_a1b.horner(y); + a2=b1_a2b.horner(y); + a=M.sqrt(W/x)*(M.sin(xx)*a1+M.cos(xx)*a2*8/x); + } + return a; + } + + return _bessel_wrap(bessel0, bessel1, 'BESSELY', 1, -1); +})(); +var besseli = (function() { + var b0_a = [1.0, 3.5156229, 3.0899424, 1.2067492, 0.2659732, 0.360768e-1, 0.45813e-2].reverse(); + var b0_b = [0.39894228, 0.1328592e-1, 0.225319e-2, -0.157565e-2, 0.916281e-2, -0.2057706e-1, 0.2635537e-1, -0.1647633e-1, 0.392377e-2].reverse(); + function bessel0(x) { + if(x <= 3.75) return b0_a.horner(x*x/(3.75*3.75)); + return M.exp(M.abs(x))/M.sqrt(M.abs(x))*b0_b.horner(3.75/M.abs(x)); + } + + var b1_a = [0.5, 0.87890594, 0.51498869, 0.15084934, 0.2658733e-1, 0.301532e-2, 0.32411e-3].reverse(); + var b1_b = [0.39894228, -0.3988024e-1, -0.362018e-2, 0.163801e-2, -0.1031555e-1, 0.2282967e-1, -0.2895312e-1, 0.1787654e-1, -0.420059e-2].reverse(); + function bessel1(x) { + if(x < 3.75) return x * b1_a.horner(x*x/(3.75*3.75)); + return (x < 0 ? -1 : 1) * M.exp(M.abs(x))/M.sqrt(M.abs(x))*b1_b.horner(3.75/M.abs(x)); + } + + return function besseli(x, n) { + n = Math.round(n); + if(n === 0) return bessel0(x); + if(n == 1) return bessel1(x); + if(n < 0) throw 'BESSELI Order (' + n + ') must be nonnegative'; + if(M.abs(x) === 0) return 0; + + var ret, j, tox = 2 / M.abs(x), m, bip, bi, bim; + m=2*M.round((n+M.round(M.sqrt(40*n)))/2); + bip=ret=0.0; + bi=1.0; + for (j=m;j>0;j--) { + bim=j*tox*bi + bip; + bip=bi; bi=bim; + if (M.abs(bi) > 1E10) { + bi *= 1E-10; + bip *= 1E-10; + ret *= 1E-10; + } + if(j == n) ret = bip; + } + ret *= besseli(x, 0) / bi; + return x < 0 && (n%2) ? -ret : ret; + }; + +})(); + +var besselk = (function() { + var b0_a = [-0.57721566, 0.42278420, 0.23069756, 0.3488590e-1, 0.262698e-2, 0.10750e-3, 0.74e-5].reverse(); + var b0_b = [1.25331414, -0.7832358e-1, 0.2189568e-1, -0.1062446e-1, 0.587872e-2, -0.251540e-2, 0.53208e-3].reverse(); + function bessel0(x) { + if(x <= 2) return -M.log(x/2)*besseli(x,0) + b0_a.horner(x*x/4); + return M.exp(-x)/M.sqrt(x)*b0_b.horner(2/x); + } + + var b1_a = [1.0, 0.15443144, -0.67278579, -0.18156897, -0.1919402e-1, -0.110404e-2, -0.4686e-4].reverse(); + var b1_b = [1.25331414, 0.23498619, -0.3655620e-1, 0.1504268e-1, -0.780353e-2, 0.325614e-2, -0.68245e-3].reverse(); + function bessel1(x) { + if(x <= 2) return M.log(x/2)*besseli(x,1) + (1/x)*b1_a.horner(x*x/4); + return M.exp(-x)/M.sqrt(x)*b1_b.horner(2/x); + } + + return _bessel_wrap(bessel0, bessel1, 'BESSELK', 2, 1); +})(); +if(typeof exports !== "undefined") { + exports.besselj = besselj; + exports.bessely = bessely; + exports.besseli = besseli; + exports.besselk = besselk; +} + diff --git a/bessel.md b/bessel.md new file mode 100644 index 0000000..ace1082 --- /dev/null +++ b/bessel.md @@ -0,0 +1,317 @@ +# Code + +The code is embedded in this document! [VOC](https://npmjs.org/package/voc) can +run through this file and generate `bessel.js`. + +I really dislike writing `Math` so I use: + +```js +var M = Math; +``` + +## Horner Method + +The methods use an approximating polynomial and evaluate using Horner's method. +In true JS form, let us define an `Array` prototype method: + +``` +Array.prototype.horner = function(v) { return this.reduce(function(z,w){return v * z + w;},0); }; +``` + +## Recurrence + +It can be shown that the four Bessel functions satisfy (on their support): + +```tex> +B_{n} (x) = \frac{2n}{x} B_{n-1}(x) - B_{n-2}(x) +``` + +So rather than go back and try to find solution for each order, we will build +solutions for `n=0` and `n=1` and then apply the recurrence. The helper: + +```js +function _bessel_iter(x, n, f0, f1, sign) { + if(!sign) sign = -1; + var tdx = 2 / x, f2; + if(n === 0) return f0; + if(n === 1) return f1; + for(var o = 1; o != n; ++o) { + f2 = f1 * o * tdx + sign * f0; + f0 = f1; f1 = f2; + } + return f1; +} +``` + +We can directly generate the JS function given the basic solutions `bessel0` and +`bessel1` by leveraging `_bessel_iter` from above. We have to add a few sanity +checks since `Y_n` is undefined at 0 and `K_n` is real only when `x>0` + +``` +function _bessel_wrap(bessel0, bessel1, name, nonzero, sign) { + return function bessel(x,n) { + if(n === 0) return bessel0(x); + if(n === 1) return bessel1(x); + if(n < 0) throw name + ': Order (' + n + ') must be nonnegative'; + if(nonzero == 1 && x === 0) throw name + ': Undefined when x == 0'; + if(nonzero == 2 && x <= 0) throw name + ': Undefined when x <= 0'; + var b0 = bessel0(x), b1 = bessel1(x); + return _bessel_iter(x, n, b0, b1, sign); + }; +} +``` + +## Individual Solutions + +To determine each individual solution, we first calculate a Chebyshev polynomial +based on the regime (lower and higher-value approximations). This module uses +the constants from the [GNU Scientific Library](https://gnu.org/s/gsl), and from +the venerable [Numerical Recipes book](http://www.nr.com/) but I have +independently verified the constants using Mathematica. + +``` +var besselj = (function() { +``` + +The constants are named `b[01]_[ab]([12][ab])?` with the blocks corresponding to +the function order (e.g. `b0_` refers to order 0), variable name in the function +and conditional. + +``` + var b0_a1a = [57568490574.0,-13362590354.0,651619640.7,-11214424.18,77392.33017,-184.9052456].reverse(); + var b0_a2a = [57568490411.0,1029532985.0,9494680.718,59272.64853,267.8532712,1.0].reverse(); + var b0_a1b = [1.0, -0.1098628627e-2, 0.2734510407e-4, -0.2073370639e-5, 0.2093887211e-6].reverse(); + var b0_a2b = [-0.1562499995e-1, 0.1430488765e-3, -0.6911147651e-5, 0.7621095161e-6, -0.934935152e-7].reverse(); +``` + +I noticed some strange oddities when leaning on `Math.PI`, so it is cached: + +``` + var W = 0.636619772; // 2 / Math.PI + + function bessel0(x) { + var a, a1, a2, y = x * x, xx = M.abs(x) - 0.785398164; +``` + +For small `x`, the direct Laurent approximation gives better results. + +``` + if(M.abs(x) < 8) { + a1 = b0_a1a.horner(y); + a2 = b0_a2a.horner(y); + a = a1/a2; + } +``` + +For larger `x`, the Chebyshev approach is taken: + +``` + else { + y = 64 / y; + a1 = b0_a1b.horner(y); + a2 = b0_a2b.horner(y); + a = M.sqrt(W/M.abs(x))*(M.cos(xx)*a1-M.sin(xx)*a2*8/M.abs(x)); + } + return a; + } +``` + +A similar approach is taken for the first-order bessel function + +``` + var b1_a1a = [72362614232.0,-7895059235.0,242396853.1,-2972611.439, 15704.48260, -30.16036606].reverse(); + var b1_a2a = [144725228442.0, 2300535178.0, 18583304.74, 99447.43394, 376.9991397, 1.0].reverse(); + var b1_a1b = [1.0, 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, -0.240337019e-6].reverse(); + var b1_a2b = [0.04687499995, -0.2002690873e-3, 0.8449199096e-5, -0.88228987e-6, 0.105787412e-6].reverse(); + function bessel1(x) { + var a, a1, a2, y = x*x, xx = M.abs(x) - 2.356194491; + if(Math.abs(x)< 8) { + a1 = x*b1_a1a.horner(y); + a2 = b1_a2a.horner(y); + a = a1 / a2; + } else { + y = 64 / y; + a1=b1_a1b.horner(y); + a2=b1_a2b.horner(y); + a=M.sqrt(W/M.abs(x))*(M.cos(xx)*a1-M.sin(xx)*a2*8/M.abs(x)); + if(x < 0) a = -a; + } + return a; + } +``` + +For large values of x, the aforementioned iteration is fine, but for small +values the expressions quickly blow up. Hence a more careful iteration is used: + +``` + return function besselj(x, n) { + n = Math.round(n); + if(n === 0) return bessel0(M.abs(x)); + if(n === 1) return bessel1(M.abs(x)); + if(n < 0) throw 'BESSELJ: Order (' + n + ') must be nonnegative'; + if(M.abs(x) === 0) return 0; + + var ret, j, tox = 2 / M.abs(x), m, jsum, sum, bjp, bj, bjm; + if(M.abs(x) > n) { + ret = _bessel_iter(x, n, bessel0(M.abs(x)), bessel1(M.abs(x)),-1); + } else { + m=2*M.floor((n+M.floor(M.sqrt(40*n)))/2); + jsum=0; + bjp=ret=sum=0.0; + bj=1.0; + for (j=m;j>0;j--) { + bjm=j*tox*bj-bjp; + bjp=bj; + bj=bjm; + if (M.abs(bj) > 1E10) { + bj *= 1E-10; + bjp *= 1E-10; + ret *= 1E-10; + sum *= 1E-10; + } + if (jsum) sum += bj; + jsum=!jsum; + if (j == n) ret=bjp; + } + sum=2.0*sum-bj; + ret /= sum; + } + return x < 0 && (n%2) ? -ret : ret; + }; +})(); +``` + +The second kind function `Y` is a bit less finicky: + +``` +var bessely = (function() { + var b0_a1a = [-2957821389.0, 7062834065.0, -512359803.6, 10879881.29, -86327.92757, 228.4622733].reverse(); + var b0_a2a = [40076544269.0, 745249964.8, 7189466.438, 47447.26470, 226.1030244, 1.0].reverse(); + var b0_a1b = [1.0, -0.1098628627e-2, 0.2734510407e-4, -0.2073370639e-5, 0.2093887211e-6].reverse(); + var b0_a2b = [-0.1562499995e-1, 0.1430488765e-3, -0.6911147651e-5, 0.7621095161e-6, -0.934945152e-7].reverse(); + + var W = 0.636619772; + function bessel0(x) { + var a, a1, a2, y = x * x, xx = x - 0.785398164; + if(x < 8) { + a1 = b0_a1a.horner(y); + a2 = b0_a2a.horner(y); + a = a1/a2 + W * besselj(x,0) * M.log(x); + } else { + y = 64 / y; + a1 = b0_a1b.horner(y); + a2 = b0_a2b.horner(y); + a = M.sqrt(W/x)*(M.sin(xx)*a1+M.cos(xx)*a2*8/x); + } + return a; + } + + var b1_a1a = [-0.4900604943e13, 0.1275274390e13, -0.5153438139e11, 0.7349264551e9, -0.4237922726e7, 0.8511937935e4].reverse(); + var b1_a2a = [0.2499580570e14, 0.4244419664e12, 0.3733650367e10, 0.2245904002e8, 0.1020426050e6, 0.3549632885e3, 1].reverse(); + var b1_a1b = [1.0, 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, -0.240337019e-6].reverse(); + var b1_a2b = [0.04687499995, -0.2002690873e-3, 0.8449199096e-5, -0.88228987e-6, 0.105787412e-6].reverse(); + function bessel1(x) { + var a, a1, a2, y = x*x, xx = x - 2.356194491; + if(x < 8) { + a1 = x*b1_a1a.horner(y); + a2 = b1_a2a.horner(y); + a = a1/a2 + W * (besselj(x,1) * M.log(x) - 1 / x); + } else { + y = 64 / y; + a1=b1_a1b.horner(y); + a2=b1_a2b.horner(y); + a=M.sqrt(W/x)*(M.sin(xx)*a1+M.cos(xx)*a2*8/x); + } + return a; + } + + return _bessel_wrap(bessel0, bessel1, 'BESSELY', 1, -1); +})(); +``` + +And the modified Bessel functions are even easier: + +``` +var besseli = (function() { + var b0_a = [1.0, 3.5156229, 3.0899424, 1.2067492, 0.2659732, 0.360768e-1, 0.45813e-2].reverse(); + var b0_b = [0.39894228, 0.1328592e-1, 0.225319e-2, -0.157565e-2, 0.916281e-2, -0.2057706e-1, 0.2635537e-1, -0.1647633e-1, 0.392377e-2].reverse(); + function bessel0(x) { + if(x <= 3.75) return b0_a.horner(x*x/(3.75*3.75)); + return M.exp(M.abs(x))/M.sqrt(M.abs(x))*b0_b.horner(3.75/M.abs(x)); + } + + var b1_a = [0.5, 0.87890594, 0.51498869, 0.15084934, 0.2658733e-1, 0.301532e-2, 0.32411e-3].reverse(); + var b1_b = [0.39894228, -0.3988024e-1, -0.362018e-2, 0.163801e-2, -0.1031555e-1, 0.2282967e-1, -0.2895312e-1, 0.1787654e-1, -0.420059e-2].reverse(); + function bessel1(x) { + if(x < 3.75) return x * b1_a.horner(x*x/(3.75*3.75)); + return (x < 0 ? -1 : 1) * M.exp(M.abs(x))/M.sqrt(M.abs(x))*b1_b.horner(3.75/M.abs(x)); + } + + return function besseli(x, n) { + n = Math.round(n); + if(n === 0) return bessel0(x); + if(n == 1) return bessel1(x); + if(n < 0) throw 'BESSELI Order (' + n + ') must be nonnegative'; + if(M.abs(x) === 0) return 0; + + var ret, j, tox = 2 / M.abs(x), m, bip, bi, bim; + m=2*M.round((n+M.round(M.sqrt(40*n)))/2); + bip=ret=0.0; + bi=1.0; + for (j=m;j>0;j--) { + bim=j*tox*bi + bip; + bip=bi; bi=bim; + if (M.abs(bi) > 1E10) { + bi *= 1E-10; + bip *= 1E-10; + ret *= 1E-10; + } + if(j == n) ret = bip; + } + ret *= besseli(x, 0) / bi; + return x < 0 && (n%2) ? -ret : ret; + }; + +})(); + +var besselk = (function() { + var b0_a = [-0.57721566, 0.42278420, 0.23069756, 0.3488590e-1, 0.262698e-2, 0.10750e-3, 0.74e-5].reverse(); + var b0_b = [1.25331414, -0.7832358e-1, 0.2189568e-1, -0.1062446e-1, 0.587872e-2, -0.251540e-2, 0.53208e-3].reverse(); + function bessel0(x) { + if(x <= 2) return -M.log(x/2)*besseli(x,0) + b0_a.horner(x*x/4); + return M.exp(-x)/M.sqrt(x)*b0_b.horner(2/x); + } + + var b1_a = [1.0, 0.15443144, -0.67278579, -0.18156897, -0.1919402e-1, -0.110404e-2, -0.4686e-4].reverse(); + var b1_b = [1.25331414, 0.23498619, -0.3655620e-1, 0.1504268e-1, -0.780353e-2, 0.325614e-2, -0.68245e-3].reverse(); + function bessel1(x) { + if(x <= 2) return M.log(x/2)*besseli(x,1) + (1/x)*b1_a.horner(x*x/4); + return M.exp(-x)/M.sqrt(x)*b1_b.horner(2/x); + } + + return _bessel_wrap(bessel0, bessel1, 'BESSELK', 2, 1); +})(); +``` + +## Export Magic + +Now we need to export: + +``` +if(typeof exports !== "undefined") { + exports.besselj = besselj; + exports.bessely = bessely; + exports.besseli = besseli; + exports.besselk = besselk; +} +``` + +# Testing + +Fortunately, these functions are available in Excel, so you can test there (or +use the same functions in Mathematica). Note that this uses Excel semantics, +which are reverse from Mathematica: + +```mathematica> +BesselJ[n, x] (* Mathematica *) +``` diff --git a/package.json b/package.json new file mode 100644 index 0000000..a887a12 --- /dev/null +++ b/package.json @@ -0,0 +1,19 @@ +{ + "name": "bessel", + "version": "0.1.1", + "author": "SheetJS", + "description": "Bessel Functions in pure JS", + "keywords": ["bessel", "math", "specfun"], + "main": "bessel.js", + "repository": { "type": "git", "url": "https://github.com/SheetJS/bessel" }, + "license": "MIT", + "scripts": { + "test": "make test" + }, + "dependencies": { + "voc":"" + }, + "devDependencies": { + "mocha":"" + } +} diff --git a/test.html b/test.html new file mode 100644 index 0000000..53c67ef --- /dev/null +++ b/test.html @@ -0,0 +1,4 @@ + + diff --git a/test.js b/test.js new file mode 100644 index 0000000..836aec3 --- /dev/null +++ b/test.js @@ -0,0 +1,2 @@ +var bessel; +describe('source', function() { it('should load', function() { bessel = require('./'); }); });