sheetjs/bessel
1
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Initial commit

Browse Source
This commit is contained in:
SheetJS 2013-12-14 02:36:05 -05:00
commit 45c1c7b614
11 changed files with 612 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitignore vendored Normal file
View File

@ -0,0 +1 @@
node_modules

2
.npmignore Normal file
View File

@ -0,0 +1,2 @@
Makefile
test.html

6
.travis.yml Normal file
View File

@ -0,0 +1,6 @@
language: node_js
node_js:
- "0.10"
- "0.8"
before_install:
- "npm install -g mocha"

9
LICENSE Normal file
View File

@ -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.

7
Makefile Normal file
View File

@ -0,0 +1,7 @@
LIBRARY=bessel
$(LIBRARY).js: bessel.md
node_modules/.bin/voc $^ > $@
test mocha:
mocha -R spec

36
README.md Normal file
View File

@ -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>
<script src="bessel.js"></script>
<script>console.log(besselj(1,2));</script>
```
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

209
bessel.js Normal file
View File

@ -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;
}

317
bessel.md Normal file
View File

@ -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 *)
```

19
package.json Normal file
View File

@ -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":""
}
}

4
test.html Normal file
View File

@ -0,0 +1,4 @@
<script src="bessel.js"></script>
<script>
console.log(besselj(1,2));
</script>

2
test.js Normal file
View File

@ -0,0 +1,2 @@
var bessel;
describe('source', function() { it('should load', function() { bessel = require('./'); }); });