diff --git a/.flowconfig b/.flowconfig
new file mode 100644
index 0000000..39d82f6
--- /dev/null
+++ b/.flowconfig
@@ -0,0 +1,19 @@
+[ignore]
+.*/node_modules/.*
+.*/dist/.*
+.*/test.js
+
+.*/bits/.*
+.*/ctest/.*
+.*/misc/.*
+.*/perf/.*
+
+.*/bessel.js
+
+[include]
+bessel.flow.js
+
+[libs]
+bits/10_types.js
+
+[options]
diff --git a/.gitignore b/.gitignore
index 3c3629e..865fd19 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1 +1,2 @@
node_modules
+misc/coverage.html
diff --git a/.jscs.json b/.jscs.json
new file mode 100644
index 0000000..5965004
--- /dev/null
+++ b/.jscs.json
@@ -0,0 +1,6 @@
+{
+ "requireCommaBeforeLineBreak": true,
+ "disallowTrailingWhitespace": true,
+ "disallowTrailingComma": true
+}
+
diff --git a/.jshintrc b/.jshintrc
new file mode 100644
index 0000000..5cd205b
--- /dev/null
+++ b/.jshintrc
@@ -0,0 +1,4 @@
+{
+ "bitwise": false,
+ "curly": false
+}
diff --git a/.travis.yml b/.travis.yml
index dba96e7..8beb852 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -1,7 +1,15 @@
language: node_js
node_js:
- - "0.11"
+ - "6"
+ - "5"
+ - "4.2"
+ - "0.12"
- "0.10"
- "0.8"
before_install:
- - "npm install -g mocha"
+ - "npm install -g npm@next"
+ - "npm install -g mocha voc"
+ - "npm install blanket"
+ - "npm install coveralls mocha-lcov-reporter"
+after_success:
+ - "make coveralls-spin"
diff --git a/Makefile b/Makefile
index bc8c57c..756e598 100644
--- a/Makefile
+++ b/Makefile
@@ -1,7 +1,66 @@
-LIBRARY=bessel
+LIB=bessel
+REQS=
+ADDONS=
+AUXTARGETS=
+HTMLLINT=index.html
-$(LIBRARY).js: bessel.md
- node_modules/.bin/voc $^ > $@
+ULIB=$(shell echo $(LIB) | tr a-z A-Z)
+DEPS=$(sort $(wildcard bits/*.js))
+TARGET=$(LIB).js
+FLOWTARGET=$(LIB).flow.js
-test mocha:
- mocha -R spec
+## Main Targets
+
+.PHONY: all
+all: $(TARGET) $(AUXTARGETS) ## Build library and auxiliary scripts
+
+$(TARGET) $(AUXTARGETS): %.js : %.flow.js
+ node -e 'process.stdout.write(require("fs").readFileSync("$<","utf8").replace(/^[ \t]*\/\*[:#][^*]*\*\/\s*(\n)?/gm,"").replace(/\/\*[:#][^*]*\*\//gm,""))' > $@
+
+$(FLOWTARGET): $(DEPS)
+ cat $^ | tr -d '\15\32' > $@
+
+bits/01_version.js: package.json
+ echo "BESSEL.version = '"`grep version package.json | awk '{gsub(/[^0-9a-z\.-]/,"",$$2); print $$2}'`"';" > $@
+
+.PHONY: clean
+clean: ## Remove targets and build artifacts
+ rm -f $(TARGET) $(FLOWTARGET)
+
+## Testing
+
+.PHONY: test mocha
+test mocha: test.js $(TARGET) ## Run test suite
+ mocha -R spec -t 20000
+
+## Code Checking
+
+.PHONY: lint
+lint: $(TARGET) $(AUXTARGETS) ## Run jshint and jscs checks
+ @jshint --show-non-errors $(TARGET) $(AUXTARGETS)
+ @jshint --show-non-errors package.json
+ @jshint --show-non-errors --extract=always $(HTMLLINT)
+ @jscs $(TARGET) $(AUXTARGETS)
+
+.PHONY: flow
+flow: lint ## Run flow checker
+ @flow check --all --show-all-errors
+
+.PHONY: cov
+cov: misc/coverage.html ## Run coverage test
+
+misc/coverage.html: $(TARGET) test.js
+ mocha --require blanket -R html-cov -t 20000 > $@
+
+.PHONY: coveralls
+coveralls: ## Coverage Test + Send to coveralls.io
+ mocha --require blanket --reporter mocha-lcov-reporter -t 20000 | node ./node_modules/coveralls/bin/coveralls.js
+
+
+.PHONY: help
+help:
+ @grep -hE '(^[a-zA-Z_-][ a-zA-Z_-]*:.*?|^#[#*])' $(MAKEFILE_LIST) | bash misc/help.sh
+
+#* To show a spinner, append "-spin" to any target e.g. cov-spin
+%-spin:
+ @make $* & bash misc/spin.sh $$!
diff --git a/README.md b/README.md
index 3ce6439..3e8133d 100644
--- a/README.md
+++ b/README.md
@@ -1,36 +1,92 @@
-# Bessel.JS
+# bessel
-Pure-JS implementation of the Bessel functions (J,Y,I,K), for node and browser
+Pure-JS implementation of Bessel functions J,Y,I,K (for the browser and nodejs).
+Emphasis on correctness and performance for integer order.
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
+ - `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 second kind
-# Usage
+## Installation
-The functions `besselj`, `bessely`, `besseli`, `besselk` are exposed when you include
-the script `bessel.js`:
+With [npm](https://www.npmjs.org/package/bessel):
-```html>
+```bash
+$ npm install bessel
+```
+
+In the browser:
+
+```html
-
```
-See `test.html` for an example
+The browser exposes a variable `BESSEL`
-In node, those four functions are exported:
+The script will manipulate `module.exports` if available (e.g. in a CommonJS
+`require` context). This is not always desirable. To prevent the behavior,
+define `DO_NOT_EXPORT_BESSEL`
-```js>
-var besselj01 = require('bessel').besselj(0,1);
+## Usage
+
+In all cases, the relevant function takes two arguments (`value`, `order`).
+
+The return value is a JS number. `NaN` signals an error in calculation.
+
+- `BESSEL.besselj(x, n)` computes `J_{n}(x)`
+
+- `BESSEL.bessely(x, n)` computes `Y_{n}(x)`
+
+- `BESSEL.besseli(x, n)` computes `I_{n}(x)`
+
+- `BESSEL.besselk(x, n)` computes `K_{n}(x)`
+
+For example:
+
+```js
+> // var BESSEL = require('bessel'); // uncomment this line if in node
+> BESSEL.besselj(1.5,0) // 0.5118276712499389
+> BESSEL.bessely(1.5,0) // 0.38244892476502895
+> BESSEL.besseli(1.5,0) // 1.6467232021476754
+> BESSEL.besselk(1.5,0) // 0.2138055693236539
+
+> var Y = BESSEL.bessely
+> Y(Math.PI, 5) + Y(Math.PI, 3) - (2 * 4 / Math.PI) * Y(Math.PI, 4) // 0
```
-Each function follows Excel semantics `(value, function-order)`. For example,
+## Testing
-```js>
-bessel.besselj(1.5, 1)
+`make test` will run the nodejs-based test.
+
+To generate the `excel.tsv` test cases, make a 6-column Excel Sheet:
+
+| x | n | `BESSELI` | `BESSELJ` | `BESSELK` | `BESSELY` |
+|---|---|:--------------:|:--------------:|:--------------:|:--------------:|
+| x | n |`BESSELI(A1,B1)`|`BESSELJ(A1,B1)`|`BESSELK(A1,B1)`|`BESSELY(A1,B1)`|
+
+To generate the `mma.tsv` test cases, use the Mathematica Bessel Functions:
+
+```mathematica
+(* Bessel_[value, order] *)
+F[x_,n_]:={x/2,n,BesselI[n,x/2], BesselJ[n,x/2], BesselK[n,x/2], BesselY[n,x/2]}
+Do[Print[ExportString[N[F[x,n],10],"csv"]], {n,1,3}, {x,1,20} ]
```
-is the value of the bessel function J1 at the point x=1.5
+Note: Each function follows Excel semantics `(value, order)`. Other platforms
+like Mathematica reverse the argument order.
+
+## License
+
+Please consult the attached LICENSE file for details. All rights not explicitly
+granted by the MIT License are reserved by the Original Author.
+
+## Badges
+
+[![Build Status](https://travis-ci.org/SheetJS/bessel.svg?branch=master)](https://travis-ci.org/SheetJS/bessel)
+
+[![Coverage Status](http://img.shields.io/coveralls/SheetJS/bessel/master.svg)](https://coveralls.io/r/SheetJS/bessel?branch=master)
+
+[![Analytics](https://ga-beacon.appspot.com/UA-36810333-1/SheetJS/bessel?pixel)](https://github.com/SheetJS/bessel)
diff --git a/bessel.flow.js b/bessel.flow.js
new file mode 100644
index 0000000..2336c91
--- /dev/null
+++ b/bessel.flow.js
@@ -0,0 +1,252 @@
+/* bessel.js (C) 2013-present SheetJS -- http://sheetjs.com */
+/* vim: set ts=2: */
+/*exported BESSEL */
+var BESSEL;
+/*:: declare var DO_NOT_EXPORT_BESSEL: any; */
+/*:: declare var define: any; */
+(function (factory) {
+ /*jshint ignore:start */
+ if(typeof DO_NOT_EXPORT_BESSEL === 'undefined') {
+ if('object' === typeof exports) {
+ factory(exports);
+ } else if ('function' === typeof define && define.amd) {
+ define(function () {
+ var module = {};
+ factory(module);
+ return module;
+ });
+ } else {
+ factory(BESSEL = {});
+ }
+ } else {
+ factory(BESSEL = {});
+ }
+ /*jshint ignore:end */
+}(function(BESSEL) {
+BESSEL.version = '0.3.0';
+var M = Math;
+
+/*::
+type BesselN = (x:number) => number;
+type BesselF = (x:number, n:number) => number;
+*/
+function _horner(arr/*:Array*/, v/*:number*/)/*:number*/ { for(var i = 0, z = 0; i < arr.length; ++i) z = v * z + arr[i]; return z; }
+function _bessel_iter(x/*:number*/, n/*:number*/, f0/*:number*/, f1/*:number*/, sign/*:number*/)/*:number*/ {
+ if(n === 0) return f0;
+ if(n === 1) return f1;
+ var tdx = 2 / x, f2 = f1;
+ for(var o = 1; o < n; ++o) {
+ f2 = f1 * o * tdx + sign * f0;
+ f0 = f1; f1 = f2;
+ }
+ return f2;
+}
+function _bessel_wrap(bessel0/*:BesselN*/, bessel1/*:BesselN*/, name/*:string*/, nonzero/*:number*/, sign/*:number*/)/*:BesselF*/ {
+ return function bessel(x/*:number*/,n/*:number*/) {
+ if(nonzero) {
+ if(x === 0) return (nonzero == 1 ? -Infinity : Infinity);
+ else if(x < 0) return NaN;
+ }
+ if(n === 0) return bessel0(x);
+ if(n === 1) return bessel1(x);
+ if(n < 0) return NaN;
+ n|=0;
+ var b0/*:number*/ = bessel0(x), b1/*:number*/ = bessel1(x);
+ return _bessel_iter(x, n, b0, b1, sign);
+ };
+}
+var besselj/*:BesselF*/ = (function() {
+ var W = 0.636619772; // 2 / Math.PI
+
+ 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();
+
+ function bessel0(x/*:number*/)/*:number*/ {
+ var a=0, a1=0, a2=0, y = x * x;
+ if(x < 8) {
+ a1 = _horner(b0_a1a, y);
+ a2 = _horner(b0_a2a, y);
+ a = a1 / a2;
+ } else {
+ var xx = x - 0.785398164;
+ y = 64 / y;
+ a1 = _horner(b0_a1b, y);
+ a2 = _horner(b0_a2b, y);
+ a = M.sqrt(W/x)*(M.cos(xx)*a1-M.sin(xx)*a2*8/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/*:number*/)/*:number*/ {
+ var a=0, a1=0, a2=0, y = x*x, xx = M.abs(x) - 2.356194491;
+ if(Math.abs(x)< 8) {
+ a1 = x*_horner(b1_a1a, y);
+ a2 = _horner(b1_a2a, y);
+ a = a1 / a2;
+ } else {
+ y = 64 / y;
+ a1=_horner(b1_a1b, y);
+ a2=_horner(b1_a2b, 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/*:number*/, n/*:number*/)/*:number*/ {
+ n = Math.round(n);
+ if(!isFinite(x)) return isNaN(x) ? x : 0;
+ if(n < 0) return ((n%2)?-1:1)*besselj(x, -n);
+ if(x < 0) return ((n%2)?-1:1)*besselj(-x, n);
+ if(n === 0) return bessel0(x);
+ if(n === 1) return bessel1(x);
+ if(x === 0) return 0;
+
+ var ret=0.0;
+ if(x > n) {
+ ret = _bessel_iter(x, n, bessel0(x), bessel1(x),-1);
+ } else {
+ var m=2*M.floor((n+M.floor(M.sqrt(40*n)))/2);
+ var jsum=false;
+ var bjp=0.0, sum=0.0;
+ var bj=1.0, bjm = 0.0;
+ var tox = 2 / x;
+ for (var 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 ret;
+ };
+})();
+var bessely/*:BesselF*/ = (function() {
+ var W = 0.636619772;
+
+ 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();
+
+ function bessel0(x/*:number*/)/*:number*/ {
+ var a=0, a1=0, a2=0, y = x * x, xx = x - 0.785398164;
+ if(x < 8) {
+ a1 = _horner(b0_a1a, y);
+ a2 = _horner(b0_a2a, y);
+ a = a1/a2 + W * besselj(x,0) * M.log(x);
+ } else {
+ y = 64 / y;
+ a1 = _horner(b0_a1b, y);
+ a2 = _horner(b0_a2b, 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/*:number*/)/*:number*/ {
+ var a=0, a1=0, a2=0, y = x*x, xx = x - 2.356194491;
+ if(x < 8) {
+ a1 = x*_horner(b1_a1a, y);
+ a2 = _horner(b1_a2a, y);
+ a = a1/a2 + W * (besselj(x,1) * M.log(x) - 1 / x);
+ } else {
+ y = 64 / y;
+ a1=_horner(b1_a1b, y);
+ a2=_horner(b1_a2b, 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/*:BesselF*/ = (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/*:number*/)/*:number*/ {
+ if(x <= 3.75) return _horner(b0_a, x*x/(3.75*3.75));
+ return M.exp(M.abs(x))/M.sqrt(M.abs(x))*_horner(b0_b, 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/*:number*/)/*:number*/ {
+ if(x < 3.75) return x * _horner(b1_a, x*x/(3.75*3.75));
+ return (x < 0 ? -1 : 1) * M.exp(M.abs(x))/M.sqrt(M.abs(x))*_horner(b1_b, 3.75/M.abs(x));
+ }
+
+ return function besseli(x/*:number*/, n/*:number*/)/*:number*/ {
+ n = Math.round(n);
+ if(n === 0) return bessel0(x);
+ if(n === 1) return bessel1(x);
+ if(n < 0) return NaN;
+ if(M.abs(x) === 0) return 0;
+ if(x == Infinity) return Infinity;
+
+ var ret = 0.0, j, tox = 2 / M.abs(x), bip = 0.0, bi=1.0, bim=0.0;
+ var m=2*M.round((n+M.round(M.sqrt(40*n)))/2);
+ 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/*:BesselF*/ = (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/*:number*/)/*:number*/ {
+ if(x <= 2) return -M.log(x/2) * besseli(x,0) + _horner(b0_a, x*x/4);
+ return M.exp(-x) / M.sqrt(x) * _horner(b0_b, 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/*:number*/)/*:number*/ {
+ if(x <= 2) return M.log(x/2) * besseli(x,1) + (1/x) * _horner(b1_a, x*x/4);
+ return M.exp(-x)/M.sqrt(x)*_horner(b1_b, 2/x);
+ }
+
+ return _bessel_wrap(bessel0, bessel1, 'BESSELK', 2, 1);
+})();
+BESSEL.besselj = besselj;
+BESSEL.bessely = bessely;
+BESSEL.besseli = besseli;
+BESSEL.besselk = besselk;
+}));
diff --git a/bessel.js b/bessel.js
index c260b52..c69f73e 100644
--- a/bessel.js
+++ b/bessel.js
@@ -1,55 +1,85 @@
+/* bessel.js (C) 2013-present SheetJS -- http://sheetjs.com */
+/* vim: set ts=2: */
+/*exported BESSEL */
+var BESSEL;
+(function (factory) {
+ /*jshint ignore:start */
+ if(typeof DO_NOT_EXPORT_BESSEL === 'undefined') {
+ if('object' === typeof exports) {
+ factory(exports);
+ } else if ('function' === typeof define && define.amd) {
+ define(function () {
+ var module = {};
+ factory(module);
+ return module;
+ });
+ } else {
+ factory(BESSEL = {});
+ }
+ } else {
+ factory(BESSEL = {});
+ }
+ /*jshint ignore:end */
+}(function(BESSEL) {
+BESSEL.version = '0.3.0';
var M = Math;
-function _horner(arr, v) { return arr.reduce(function(z,w){return v * z + w;},0); };
+
+function _horner(arr, v) { for(var i = 0, z = 0; i < arr.length; ++i) z = v * z + arr[i]; return z; }
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) {
+ var tdx = 2 / x, f2 = f1;
+ for(var o = 1; o < n; ++o) {
f2 = f1 * o * tdx + sign * f0;
f0 = f1; f1 = f2;
}
- return f1;
+ return f2;
}
function _bessel_wrap(bessel0, bessel1, name, nonzero, sign) {
return function bessel(x,n) {
+ if(nonzero) {
+ if(x === 0) return (nonzero == 1 ? -Infinity : Infinity);
+ else if(x < 0) return NaN;
+ }
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';
+ if(n < 0) return NaN;
+ n|=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
+ 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();
+
function bessel0(x) {
- var a, a1, a2, y = x * x, xx = M.abs(x) - 0.785398164;
- if(M.abs(x) < 8) {
+ var a=0, a1=0, a2=0, y = x * x;
+ if(x < 8) {
a1 = _horner(b0_a1a, y);
a2 = _horner(b0_a2a, y);
- a = a1/a2;
- }
- else {
+ a = a1 / a2;
+ } else {
+ var xx = x - 0.785398164;
y = 64 / y;
a1 = _horner(b0_a1b, y);
a2 = _horner(b0_a2b, y);
- a = M.sqrt(W/M.abs(x))*(M.cos(xx)*a1-M.sin(xx)*a2*8/M.abs(x));
+ a = M.sqrt(W/x)*(M.cos(xx)*a1-M.sin(xx)*a2*8/x);
}
return a;
}
- var b1_a1a = [72362614232.0,-7895059235.0,242396853.1,-2972611.439, 15704.48260, -30.16036606].reverse();
+
+ 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;
+ var a=0, a1=0, a2=0, y = x*x, xx = M.abs(x) - 2.356194491;
if(Math.abs(x)< 8) {
a1 = x*_horner(b1_a1a, y);
a2 = _horner(b1_a2a, y);
@@ -63,22 +93,26 @@ var besselj = (function() {
}
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;
+ if(!isFinite(x)) return isNaN(x) ? x : 0;
+ if(n < 0) return ((n%2)?-1:1)*besselj(x, -n);
+ if(x < 0) return ((n%2)?-1:1)*besselj(-x, n);
+ if(n === 0) return bessel0(x);
+ if(n === 1) return bessel1(x);
+ if(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);
+ var ret=0.0;
+ if(x > n) {
+ ret = _bessel_iter(x, n, bessel0(x), bessel1(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--) {
+ var m=2*M.floor((n+M.floor(M.sqrt(40*n)))/2);
+ var jsum=false;
+ var bjp=0.0, sum=0.0;
+ var bj=1.0, bjm = 0.0;
+ var tox = 2 / x;
+ for (var j=m;j>0;j--) {
bjm=j*tox*bj-bjp;
bjp=bj;
bj=bjm;
@@ -95,18 +129,19 @@ var besselj = (function() {
sum=2.0*sum-bj;
ret /= sum;
}
- return x < 0 && (n%2) ? -ret : ret;
+ return ret;
};
})();
var bessely = (function() {
+ var W = 0.636619772;
+
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;
+ var a=0, a1=0, a2=0, y = x * x, xx = x - 0.785398164;
if(x < 8) {
a1 = _horner(b0_a1a, y);
a2 = _horner(b0_a2a, y);
@@ -124,8 +159,9 @@ var bessely = (function() {
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;
+ var a=0, a1=0, a2=0, y = x*x, xx = x - 2.356194491;
if(x < 8) {
a1 = x*_horner(b1_a1a, y);
a2 = _horner(b1_a2a, y);
@@ -144,6 +180,7 @@ var bessely = (function() {
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 _horner(b0_a, x*x/(3.75*3.75));
return M.exp(M.abs(x))/M.sqrt(M.abs(x))*_horner(b0_b, 3.75/M.abs(x));
@@ -151,6 +188,7 @@ var besseli = (function() {
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 * _horner(b1_a, x*x/(3.75*3.75));
return (x < 0 ? -1 : 1) * M.exp(M.abs(x))/M.sqrt(M.abs(x))*_horner(b1_b, 3.75/M.abs(x));
@@ -159,14 +197,13 @@ var besseli = (function() {
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(n === 1) return bessel1(x);
+ if(n < 0) return NaN;
if(M.abs(x) === 0) return 0;
+ if(x == Infinity) return Infinity;
- 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;
+ var ret = 0.0, j, tox = 2 / M.abs(x), bip = 0.0, bi=1.0, bim=0.0;
+ var m=2*M.round((n+M.round(M.sqrt(40*n)))/2);
for (j=m;j>0;j--) {
bim=j*tox*bi + bip;
bip=bi; bi=bim;
@@ -186,24 +223,24 @@ var besseli = (function() {
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) + _horner(b0_a, x*x/4);
- return M.exp(-x)/M.sqrt(x)*_horner(b0_b, 2/x);
+ if(x <= 2) return -M.log(x/2) * besseli(x,0) + _horner(b0_a, x*x/4);
+ return M.exp(-x) / M.sqrt(x) * _horner(b0_b, 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)*_horner(b1_a, x*x/4);
+ if(x <= 2) return M.log(x/2) * besseli(x,1) + (1/x) * _horner(b1_a, x*x/4);
return M.exp(-x)/M.sqrt(x)*_horner(b1_b, 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;
-}
-
+BESSEL.besselj = besselj;
+BESSEL.bessely = bessely;
+BESSEL.besseli = besseli;
+BESSEL.besselk = besselk;
+}));
diff --git a/bessel.md b/bessel.md
index b532abc..c31f947 100644
--- a/bessel.md
+++ b/bessel.md
@@ -14,7 +14,7 @@ var M = Math;
The methods use an approximating polynomial and evaluate using Horner's method:
```
-function _horner(arr, v) { return arr.reduce(function(z,w){return v * z + w;},0); };
+function _horner(arr/*:Array*/, v/*:number*/) { return arr.reduce(function(z,w){return v * z + w;},0); }
```
## Recurrence
@@ -29,7 +29,7 @@ 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) {
+function _bessel_iter(x/*:number*/, n/*:number*/, f0/*:number*/, f1/*:number*/, sign/*:?number*/) {
if(!sign) sign = -1;
var tdx = 2 / x, f2;
if(n === 0) return f0;
@@ -48,7 +48,7 @@ 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) {
+ return function bessel(x/*:number*/,n/*:number*/) {
if(n === 0) return bessel0(x);
if(n === 1) return bessel1(x);
if(n < 0) throw name + ': Order (' + n + ') must be nonnegative';
@@ -143,7 +143,7 @@ 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) {
+ return function besselj(x/*:number*/, n/*:number*/) {
n = Math.round(n);
if(n === 0) return bessel0(M.abs(x));
if(n === 1) return bessel1(M.abs(x));
@@ -246,7 +246,7 @@ var besseli = (function() {
return (x < 0 ? -1 : 1) * M.exp(M.abs(x))/M.sqrt(M.abs(x))*_horner(b1_b, 3.75/M.abs(x));
}
- return function besseli(x, n) {
+ return function besseli(x/*:number*/, n/*:number*/) {
n = Math.round(n);
if(n === 0) return bessel0(x);
if(n == 1) return bessel1(x);
diff --git a/bits/00_header.js b/bits/00_header.js
new file mode 100644
index 0000000..1b6f408
--- /dev/null
+++ b/bits/00_header.js
@@ -0,0 +1,25 @@
+/* bessel.js (C) 2013-present SheetJS -- http://sheetjs.com */
+/* vim: set ts=2: */
+/*exported BESSEL */
+var BESSEL;
+/*:: declare var DO_NOT_EXPORT_BESSEL: any; */
+/*:: declare var define: any; */
+(function (factory) {
+ /*jshint ignore:start */
+ if(typeof DO_NOT_EXPORT_BESSEL === 'undefined') {
+ if('object' === typeof exports) {
+ factory(exports);
+ } else if ('function' === typeof define && define.amd) {
+ define(function () {
+ var module = {};
+ factory(module);
+ return module;
+ });
+ } else {
+ factory(BESSEL = {});
+ }
+ } else {
+ factory(BESSEL = {});
+ }
+ /*jshint ignore:end */
+}(function(BESSEL) {
diff --git a/bits/01_version.js b/bits/01_version.js
new file mode 100644
index 0000000..025a37a
--- /dev/null
+++ b/bits/01_version.js
@@ -0,0 +1 @@
+BESSEL.version = '0.3.0';
diff --git a/bits/02_helpers.js b/bits/02_helpers.js
new file mode 100644
index 0000000..726d059
--- /dev/null
+++ b/bits/02_helpers.js
@@ -0,0 +1,2 @@
+var M = Math;
+
diff --git a/bits/10_types.js b/bits/10_types.js
new file mode 100644
index 0000000..3b845db
--- /dev/null
+++ b/bits/10_types.js
@@ -0,0 +1,4 @@
+/*::
+type BesselN = (x:number) => number;
+type BesselF = (x:number, n:number) => number;
+*/
diff --git a/bits/20_helpers.js b/bits/20_helpers.js
new file mode 100644
index 0000000..441adec
--- /dev/null
+++ b/bits/20_helpers.js
@@ -0,0 +1 @@
+function _horner(arr/*:Array*/, v/*:number*/)/*:number*/ { for(var i = 0, z = 0; i < arr.length; ++i) z = v * z + arr[i]; return z; }
diff --git a/bits/30_iter.js b/bits/30_iter.js
new file mode 100644
index 0000000..308feba
--- /dev/null
+++ b/bits/30_iter.js
@@ -0,0 +1,10 @@
+function _bessel_iter(x/*:number*/, n/*:number*/, f0/*:number*/, f1/*:number*/, sign/*:number*/)/*:number*/ {
+ if(n === 0) return f0;
+ if(n === 1) return f1;
+ var tdx = 2 / x, f2 = f1;
+ for(var o = 1; o < n; ++o) {
+ f2 = f1 * o * tdx + sign * f0;
+ f0 = f1; f1 = f2;
+ }
+ return f2;
+}
diff --git a/bits/40_wrapper.js b/bits/40_wrapper.js
new file mode 100644
index 0000000..c12ab9d
--- /dev/null
+++ b/bits/40_wrapper.js
@@ -0,0 +1,14 @@
+function _bessel_wrap(bessel0/*:BesselN*/, bessel1/*:BesselN*/, name/*:string*/, nonzero/*:number*/, sign/*:number*/)/*:BesselF*/ {
+ return function bessel(x/*:number*/,n/*:number*/) {
+ if(nonzero) {
+ if(x === 0) return (nonzero == 1 ? -Infinity : Infinity);
+ else if(x < 0) return NaN;
+ }
+ if(n === 0) return bessel0(x);
+ if(n === 1) return bessel1(x);
+ if(n < 0) return NaN;
+ n|=0;
+ var b0/*:number*/ = bessel0(x), b1/*:number*/ = bessel1(x);
+ return _bessel_iter(x, n, b0, b1, sign);
+ };
+}
diff --git a/bits/50_besselj.js b/bits/50_besselj.js
new file mode 100644
index 0000000..0d18dc3
--- /dev/null
+++ b/bits/50_besselj.js
@@ -0,0 +1,83 @@
+var besselj/*:BesselF*/ = (function() {
+ var W = 0.636619772; // 2 / Math.PI
+
+ 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();
+
+ function bessel0(x/*:number*/)/*:number*/ {
+ var a=0, a1=0, a2=0, y = x * x;
+ if(x < 8) {
+ a1 = _horner(b0_a1a, y);
+ a2 = _horner(b0_a2a, y);
+ a = a1 / a2;
+ } else {
+ var xx = x - 0.785398164;
+ y = 64 / y;
+ a1 = _horner(b0_a1b, y);
+ a2 = _horner(b0_a2b, y);
+ a = M.sqrt(W/x)*(M.cos(xx)*a1-M.sin(xx)*a2*8/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/*:number*/)/*:number*/ {
+ var a=0, a1=0, a2=0, y = x*x, xx = M.abs(x) - 2.356194491;
+ if(Math.abs(x)< 8) {
+ a1 = x*_horner(b1_a1a, y);
+ a2 = _horner(b1_a2a, y);
+ a = a1 / a2;
+ } else {
+ y = 64 / y;
+ a1=_horner(b1_a1b, y);
+ a2=_horner(b1_a2b, 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/*:number*/, n/*:number*/)/*:number*/ {
+ n = Math.round(n);
+ if(!isFinite(x)) return isNaN(x) ? x : 0;
+ if(n < 0) return ((n%2)?-1:1)*besselj(x, -n);
+ if(x < 0) return ((n%2)?-1:1)*besselj(-x, n);
+ if(n === 0) return bessel0(x);
+ if(n === 1) return bessel1(x);
+ if(x === 0) return 0;
+
+ var ret=0.0;
+ if(x > n) {
+ ret = _bessel_iter(x, n, bessel0(x), bessel1(x),-1);
+ } else {
+ var m=2*M.floor((n+M.floor(M.sqrt(40*n)))/2);
+ var jsum=false;
+ var bjp=0.0, sum=0.0;
+ var bj=1.0, bjm = 0.0;
+ var tox = 2 / x;
+ for (var 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 ret;
+ };
+})();
diff --git a/bits/60_bessely.js b/bits/60_bessely.js
new file mode 100644
index 0000000..c8a6013
--- /dev/null
+++ b/bits/60_bessely.js
@@ -0,0 +1,45 @@
+var bessely/*:BesselF*/ = (function() {
+ var W = 0.636619772;
+
+ 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();
+
+ function bessel0(x/*:number*/)/*:number*/ {
+ var a=0, a1=0, a2=0, y = x * x, xx = x - 0.785398164;
+ if(x < 8) {
+ a1 = _horner(b0_a1a, y);
+ a2 = _horner(b0_a2a, y);
+ a = a1/a2 + W * besselj(x,0) * M.log(x);
+ } else {
+ y = 64 / y;
+ a1 = _horner(b0_a1b, y);
+ a2 = _horner(b0_a2b, 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/*:number*/)/*:number*/ {
+ var a=0, a1=0, a2=0, y = x*x, xx = x - 2.356194491;
+ if(x < 8) {
+ a1 = x*_horner(b1_a1a, y);
+ a2 = _horner(b1_a2a, y);
+ a = a1/a2 + W * (besselj(x,1) * M.log(x) - 1 / x);
+ } else {
+ y = 64 / y;
+ a1=_horner(b1_a1b, y);
+ a2=_horner(b1_a2b, 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);
+})();
diff --git a/bits/70_besseli.js b/bits/70_besseli.js
new file mode 100644
index 0000000..c67f1a9
--- /dev/null
+++ b/bits/70_besseli.js
@@ -0,0 +1,43 @@
+var besseli/*:BesselF*/ = (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/*:number*/)/*:number*/ {
+ if(x <= 3.75) return _horner(b0_a, x*x/(3.75*3.75));
+ return M.exp(M.abs(x))/M.sqrt(M.abs(x))*_horner(b0_b, 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/*:number*/)/*:number*/ {
+ if(x < 3.75) return x * _horner(b1_a, x*x/(3.75*3.75));
+ return (x < 0 ? -1 : 1) * M.exp(M.abs(x))/M.sqrt(M.abs(x))*_horner(b1_b, 3.75/M.abs(x));
+ }
+
+ return function besseli(x/*:number*/, n/*:number*/)/*:number*/ {
+ n = Math.round(n);
+ if(n === 0) return bessel0(x);
+ if(n === 1) return bessel1(x);
+ if(n < 0) return NaN;
+ if(M.abs(x) === 0) return 0;
+ if(x == Infinity) return Infinity;
+
+ var ret = 0.0, j, tox = 2 / M.abs(x), bip = 0.0, bi=1.0, bim=0.0;
+ var m=2*M.round((n+M.round(M.sqrt(40*n)))/2);
+ 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;
+ };
+
+})();
+
diff --git a/bits/80_besselk.js b/bits/80_besselk.js
new file mode 100644
index 0000000..ec35b61
--- /dev/null
+++ b/bits/80_besselk.js
@@ -0,0 +1,19 @@
+var besselk/*:BesselF*/ = (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/*:number*/)/*:number*/ {
+ if(x <= 2) return -M.log(x/2) * besseli(x,0) + _horner(b0_a, x*x/4);
+ return M.exp(-x) / M.sqrt(x) * _horner(b0_b, 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/*:number*/)/*:number*/ {
+ if(x <= 2) return M.log(x/2) * besseli(x,1) + (1/x) * _horner(b1_a, x*x/4);
+ return M.exp(-x)/M.sqrt(x)*_horner(b1_b, 2/x);
+ }
+
+ return _bessel_wrap(bessel0, bessel1, 'BESSELK', 2, 1);
+})();
diff --git a/bits/90_exports.js b/bits/90_exports.js
new file mode 100644
index 0000000..4940334
--- /dev/null
+++ b/bits/90_exports.js
@@ -0,0 +1,4 @@
+BESSEL.besselj = besselj;
+BESSEL.bessely = bessely;
+BESSEL.besseli = besseli;
+BESSEL.besselk = besselk;
diff --git a/bits/99_footer.js b/bits/99_footer.js
new file mode 100644
index 0000000..be4600a
--- /dev/null
+++ b/bits/99_footer.js
@@ -0,0 +1 @@
+}));
diff --git a/index.html b/index.html
new file mode 100644
index 0000000..ce125bf
--- /dev/null
+++ b/index.html
@@ -0,0 +1,84 @@
+
+
+
+
+
+ Bessel Live Demo
+
+
+
+
+ Bessel Live Demo
+ Source Code
+ Issues? Something look weird? Click here and report an issue
+
+
+
+
+
+
+
diff --git a/misc/help.sh b/misc/help.sh
new file mode 100755
index 0000000..7a1c4c4
--- /dev/null
+++ b/misc/help.sh
@@ -0,0 +1,42 @@
+#!/bin/bash
+# make_help.sh -- process listing of targets and special items in Makefile
+# Copyright (C) 2016-present SheetJS
+#
+# usage in makefile: pipe the output of the following command:
+# @grep -hE '(^[a-zA-Z_-][ a-zA-Z_-]*:.*?|^#[#*])' $(MAKEFILE_LIST)
+#
+# lines starting with "## " are treated as subtitles
+# lines starting with "#* " are treated as plaintext comments
+# multiple targets with "## " after the ":" are rendered as separate targets
+# if the presumed default target is labeled, it will be assigned a unique color
+
+awk '
+BEGIN{recipes=0;}
+ !/#[#*] .*$/ {next;}
+ {multi=0; isrecipe=0;}
+ /^[^#]*:/ {isrecipe=1; ++recipes;}
+ /^[^ :]* .*:/ {multi=1}
+ multi==0 && isrecipe>0 { if(recipes > 1) print; else print $0, "[default]"; next}
+ isrecipe == 0 {print; next}
+ multi>0 {
+ k=split($0, msg, "##"); m=split($0, a, ":"); n=split(a[1], b, " ");
+ for(i=1; i<=n; ++i) print b[i] ":", "##" msg[2], (recipes==1 && i==1 ? "[default]" : "")
+ }
+END {}
+' | if [[ -t 1 ]]; then
+awk '
+BEGIN {FS = ":.*?## "}
+ {color=36}
+ /\[default\]/ {color=35}
+ NF==1 && /^##/ {color=34}
+ NF==1 && /^#\*/ {color=20; $1 = substr($1, 4)}
+ {printf "\033[" color "m%-20s\033[0m %s\n", $1, $2;}
+END{}' -
+else
+awk '
+BEGIN {FS = ":.*?## "}
+ /^#\* / {$1 = substr($1, 4)}
+ {printf "%-20s %s\n", $1, $2;}
+END{}' -
+fi
+
diff --git a/misc/spin.sh b/misc/spin.sh
new file mode 100755
index 0000000..471dfee
--- /dev/null
+++ b/misc/spin.sh
@@ -0,0 +1,14 @@
+#!/bin/bash
+# spin.sh -- show a spinner (for coverage test)
+# Copyright (C) 2014-present SheetJS
+
+wpid=$1
+delay=1
+str="|/-\\"
+while [ $(ps -a|awk '$1=='$wpid' {print $1}') ]; do
+ t=${str#?}
+ printf " [%c]" "$str"
+ str=$t${str%"$t"}
+ sleep $delay
+ printf "\b\b\b\b"
+done
diff --git a/package.json b/package.json
index ec15861..6176285 100644
--- a/package.json
+++ b/package.json
@@ -1,19 +1,26 @@
{
"name": "bessel",
- "version": "0.2.0",
+ "version": "0.3.0",
"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",
+ "description": "Pure-JS Bessel Functions",
+ "keywords": [ "bessel", "math", "specfun" ],
+ "main": "./bessel",
+ "dependencies": {
+ },
+ "devDependencies": {
+ "mocha":"",
+ "voc":""
+ },
+ "repository": { "type":"git", "url":"git://github.com/SheetJS/bessel.git" },
"scripts": {
"test": "make test"
},
- "dependencies": {
- "voc":""
+ "config": {
+ "blanket": {
+ "pattern": "bessel.js"
+ }
},
- "devDependencies": {
- "mocha":""
- }
+ "bugs": { "url": "https://github.com/SheetJS/bessel/issues" },
+ "license": "MIT",
+ "engines": { "node": ">=0.8" }
}
diff --git a/test.html b/test.html
deleted file mode 100644
index 53c67ef..0000000
--- a/test.html
+++ /dev/null
@@ -1,4 +0,0 @@
-
-
diff --git a/test.js b/test.js
index bb45acb..7dd4150 100644
--- a/test.js
+++ b/test.js
@@ -1,22 +1,124 @@
-var assert = require('assert');
-var bessel;
-describe('When using bessel functions', function() {
- before(function() {
- bessel = require('./bessel');
- });
- it('It must load besseli', function() {
- assert.equal(typeof bessel.besseli, 'function');
- });
- it('It must compute the Modified Bessel function at 1.5 with an order of 1 (0.981666)', function() {
- assert(Math.floor((bessel.besseli(1.5, 1) - 0.981666) * 10e6) < 10);
- });
- it('It must compute the Bessel function at 1.9 with an order of 2 (0.329926)', function() {
- assert(Math.floor((bessel.besselj(1.9, 2) - 0.329926) * 10e6) < 10);
- });
- it('It must compute the Modified Bessel function at 1.5 with an order of 1 (0.277388)', function() {
- assert(Math.floor((bessel.besselk(1.5, 1) - 0.277388) * 10e6) < 10);
- });
- it('It must compute Weber\'s Bessel function at 2.5 with an order of 1 (0.145918)', function() {
- assert(Math.floor((bessel.besselk(2.5, 1) - 0.981666) * 10e6) < 10);
+/* vim: set ts=2: */
+var X;
+if(typeof require !== 'undefined') {
+ assert = require('assert');
+ describe('source',function(){
+ it('should load',function(){X=require('./');});
+ it('all bessel functions', function() {
+ assert.equal(typeof X.besseli, 'function');
+ assert.equal(typeof X.besselj, 'function');
+ assert.equal(typeof X.besselk, 'function');
+ assert.equal(typeof X.bessely, 'function');
+ });
});
+} else { X = BESSEL; }
+
+function approx(a, b, t) { return Math.abs((a - b) / (Math.abs(a) + Math.abs(b))) <= t; }
+
+function chk(func, x, n, v, F, t) {
+ var V = func(x,n);
+ if(!approx(V, v, t) && v != V && !(isNaN(v) && isNaN(V))) throw new Error(F + "(" + x + "," + n + ") = " + V + " !~ " + v);
+}
+
+var funcs = ['besseli', 'besselj', 'besselk', 'bessely'];
+
+var tests = {
+'besseli': [
+ [1.5, 1, 0.981666],
+ [1.9, 2, 0.603272435],
+ [2.5, 1, 2.516716242]
+],
+'besselj': [
+ [1.5, 1, 0.557936508],
+ [1.9, 2, 0.329925829],
+ [2.5, 1, 0.497094103]
+],
+'besselk': [
+ [1.5, 1, 0.277388],
+ [1.9, 2, 0.296909301],
+ [2.5, 1, 0.073890816]
+],
+'bessely': [
+ [1.5, 1, -0.412308627],
+ [1.9, 2, -0.669878674],
+ [2.5, 1, 0.14591814]
+] };
+
+describe('correctness', function() {
+ funcs.forEach(function(F) { it(F, function() {
+ tests[F].forEach(function(t) { chk(X[F], t[0], t[1], t[2], F, 1e-6); });
+ });});
+
+ if(typeof require != 'undefined') {
+ var fs = require('fs');
+ it('excel', function() {
+ var xl = fs.readFileSync("test_files/excel.tsv", 'ascii').split("\n").map(function(l) { return l.split("\t");});
+ xl.forEach(function(l) { if(l.length < 6) return;
+ var x=Number(l[0]), n=Number(l[1]);
+ var i=Number(l[2]), j=Number(l[3]), k=Number(l[4]), y=Number(l[5]);
+ chk(X.besseli, x, n, i, 'besseli', 1e-4);
+ chk(X.besselj, x, n, j, 'besselj', 1e-4);
+ chk(X.besselk, x, n, k, 'besselk', 1e-4);
+ chk(X.bessely, x, n, y, 'bessely', 1e-4);
+ });
+ });
+ it('mma', function() {
+ var mma = fs.readFileSync("test_files/mma.tsv", 'ascii').split("\n").map(function(l) { return l.split("\t");});
+ mma.forEach(function(l) { if(l.length < 6) return;
+ var x=Number(l[0]), n=Number(l[1]);
+ var i=Number(l[2]), j=Number(l[3]), k=Number(l[4]), y=Number(l[5]);
+ chk(X.besseli, x, n, i, 'besseli', 1e-5);
+ chk(X.besselj, x, n, j, 'besselj', 1e-5);
+ chk(X.besselk, x, n, k, 'besselk', 1e-5);
+ chk(X.bessely, x, n, y, 'bessely', 1e-5);
+ });
+ });
+ }
+});
+
+var rand = [
+ 0.3946086812775882, 3.946086812775882, 39.46086812775882,
+ 0.03119419917289612, 0.3119419917289612, 3.119419917289612,
+ 0.18311903629136062, 1.8311903629136062, 18.311903629136062,
+ 0.5317681126599452, 5.317681126599452, 53.17681126599452,
+ 0.37724541457032235, 3.7724541457032235, 37.724541457032235,
+ 0.8287802926845655, 8.287802926845655, 82.87802926845655,
+ 0.7197724866379658, 7.197724866379658, 71.97724866379658
+];
+
+describe('properties', function() {
+ var t = 1e-4;
+ it('limiting behavior', function() {
+ chk(X.besselj, 0, 0, 1, 'besselj', 1e-6);
+ chk(X.besseli, 0, 0, 1, 'besseli', 1e-6);
+ chk(X.besselk, 0, 0, Infinity, 'besselk', 1e-6);
+ chk(X.bessely, 0, 0, -Infinity, 'bessely', 1e-6);
+ for(var i = 1; i < 20; ++i) {
+ chk(X.besselj, 0, i, 0, 'besselj', 1e-6);
+ chk(X.besseli, 0, i, 0, 'besseli', 1e-6);
+ }
+ });
+ it('besselj', function() {
+ var F = "BESSELJ", f = X.besselj;
+ rand.forEach(function(r) {
+ for(var n = 0; n < 20; ++n) {
+ var pp = f( r, n), pn = f( r, -n), np = f(-r, n), nn = f(-r, -n);
+ if((n%2)) {
+ /* besselj odd if n is odd */
+ if(!approx(np, -pp, t)) throw new Error(F + " np[" + np + "] != -pp[" + (-pp) + "] (" + r + "," + n + ")");
+ if(!approx(nn, -pn, t)) throw new Error(F + " nn[" + nn + "] != -pn[" + (-pn) + "] (" + r + "," + n + ")");
+ /* asymmetric in n */
+ if(!approx(pn, -pp, t)) throw new Error(F + " pn[" + pn + "] != -pp[" + (-pp) + "] (" + r + "," + n + ")");
+ if(!approx(nn, -np, t)) throw new Error(F + " nn[" + nn + "] != -pn[" + (-pn) + "] (" + r + "," + n + ")");
+ } else {
+ /* besselj even if n is even */
+ if(!approx(np, pp, t)) throw new Error(F + " np[" + np + "] != pp[" + (pp) + "] (" + r + "," + n + ")");
+ if(!approx(nn, pn, t)) throw new Error(F + " nn[" + nn + "] != pn[" + (pn) + "] (" + r + "," + n + ")");
+ /* symmetric in n */
+ if(!approx(pn, pp, t)) throw new Error(F + " pn[" + pn + "] != pp[" + (pp) + "] (" + r + "," + n + ")");
+ if(!approx(nn, np, t)) throw new Error(F + " nn[" + nn + "] != pn[" + (pn) + "] (" + r + "," + n + ")");
+ }
+ }
+ });
+ });
});
diff --git a/test_files/excel.tsv b/test_files/excel.tsv
new file mode 100644
index 0000000..0f03b09
--- /dev/null
+++ b/test_files/excel.tsv
@@ -0,0 +1,60 @@
+0.5 1 0.257894303 0.242268458 1.656441128 -1.471472392
+1 1 0.565159098 0.440050586 0.601907232 -0.781212821
+1.5 1 0.981666428 0.557936508 0.277387804 -0.412308627
+2 1 1.590636857 0.576724808 0.13986588 -0.107032432
+2.5 1 2.516716242 0.497094103 0.073890816 0.145918138
+3 1 3.953370217 0.339058958 0.040156431 0.324674424
+3.5 1 6.205834932 0.137377527 0.022239393 0.410188417
+4 1 9.759465157 -0.066043328 0.012483499 0.397925709
+4.5 1 15.38922293 -0.231060432 0.007078095 0.300997321
+5 1 24.33564185 -0.327579139 0.004044613 0.14786314
+5.5 1 38.58816378 -0.341438216 0.002325569 -0.023758243
+6 1 61.34193694 -0.276683859 0.00134392 -0.175010349
+6.5 1 97.73501304 -0.153841304 0.000779894 -0.274091281
+7 1 156.0390965 -0.004682826 0.000454182 -0.302667244
+7.5 1 249.5843679 0.135248424 0.000265297 -0.25912852
+8 1 399.8731348 0.234628387 0.000155369 -0.158060462
+8.5 1 641.6198938 0.273115291 9.11973E-05 -0.02616868
+9 1 1030.914709 0.245307412 5.3637E-05 0.104314575
+9.5 1 1658.453068 0.161262274 3.1602E-05 0.203179899
+10 1 2670.988321 0.04347225 1.86488E-05 0.249015424
+0.5 2 0.031906148 0.030604023 7.550183547 -5.441370834
+1 2 0.135747667 0.114903485 1.624838884 -1.650682613
+1.5 2 0.337834621 0.232087679 0.583655974 -0.932193761
+2 2 0.688948449 0.352834208 0.25375976 -0.617408098
+2.5 2 1.276466159 0.446059058 0.121460207 -0.381335848
+3 2 2.245212431 0.486091263 0.061510459 -0.160400399
+3.5 2 3.832012072 0.458629185 0.032307122 0.045371436
+4 2 6.422189499 0.364128143 0.017401426 0.215903599
+4.5 2 10.64151731 0.217848984 0.009545677 0.328481596
+5 2 17.50561501 0.046565119 0.005308944 0.367662879
+5.5 2 28.66258571 -0.117315483 0.00298437 0.330841237
+6 2 46.7870949 -0.242873213 0.001691968 0.229857908
+6.5 2 76.22054669 -0.307430389 0.000965899 0.088906659
+7 2 124.01131 -0.301417224 0.000554562 -0.0605266
+7.5 2 201.6054825 -0.230273411 0.000319924 -0.186414216
+8 2 327.5958389 -0.112993711 0.000185313 -0.263036605
+8.5 2 532.1925545 0.022323169 0.000107716 -0.276362442
+9 2 864.4962206 0.144846369 6.28006E-05 -0.226755682
+9.5 2 1404.333008 0.2278787 3.67109E-05 -0.128435911
+10 2 2281.518997 0.254630214 2.15098E-05 -0.005868083
+0.5 3 0.002645112 0.00256373 62.0579095 -42.05949428
+1 3 0.022168424 0.019563354 7.101262769 -5.821517632
+1.5 3 0.080774114 0.060963951 1.833803735 -2.073541401
+2 3 0.21273996 0.12894325 0.6473854 -1.127783765
+2.5 3 0.474370413 0.216600403 0.268227146 -0.756055495
+3 3 0.959753625 0.309062865 0.122170376 -0.538541623
+3.5 3 1.826392593 0.386770113 0.059161818 -0.358335346
+4 3 3.337275843 0.430171471 0.029884924 -0.18202211
+4.5 3 5.930096274 0.424703973 0.015563141 -0.009013679
+5 3 10.3311502 0.364831234 0.008291768 0.146267163
+5.5 3 17.7426483 0.256117865 0.00449602 0.264370052
+6 3 30.15054042 0.114768384 0.002471898 0.328248955
+6.5 3 50.83005827 -0.035346628 0.001374294 0.328803071
+7 3 85.17548647 -0.167555588 0.000771075 0.268080616
+7.5 3 142.0614436 -0.25806091 0.000435923 0.159707605
+8 3 236.0752263 -0.291125242 0.000248026 0.026542159
+8.5 3 391.1763672 -0.26261027 0.000141887 -0.103884234
+9 3 646.6942118 -0.180931248 8.15484E-05 -0.205094878
+9.5 3 1067.155009 -0.065313348 4.70593E-05 -0.257258177
+10 3 1758.380739 0.058379836 2.72527E-05 -0.251362657
diff --git a/test_files/mma.tsv b/test_files/mma.tsv
new file mode 100644
index 0000000..5bd4c8c
--- /dev/null
+++ b/test_files/mma.tsv
@@ -0,0 +1,60 @@
+0.5 1 0.2578943054 0.2422684577 1.65644112 -1.471472393
+1 1 0.565159104 0.4400505857 0.6019072302 -0.7812128213
+1.5 1 0.9816664286 0.5579365079 0.2773878005 -0.412308627
+2 1 1.590636855 0.5767248078 0.1398658818 -0.1070324315
+2.5 1 2.516716245 0.4970941025 0.07389081635 0.145918138
+3 1 3.953370217 0.3390589585 0.04015643113 0.3246744248
+3.5 1 6.205834922 0.1373775274 0.02223939293 0.4101884179
+4 1 9.759465154 -0.06604332802 0.01248349889 0.3979257106
+4.5 1 15.38922275 -0.2310604319 0.007078094909 0.3009973231
+5 1 24.33564214 -0.3275791376 0.004044613445 0.1478631434
+5.5 1 38.58816462 -0.3414382154 0.002325569009 -0.02375823896
+6 1 61.34193678 -0.2766838581 0.001343919718 -0.1750103443
+6.5 1 97.73501077 -0.1538413014 0.0007798943982 -0.274091274
+7 1 156.0390929 -0.004682823482 0.0004541824869 -0.302667237
+7.5 1 249.5843654 0.1352484276 0.0002652973901 -0.2591285105
+8 1 399.8731368 0.2346363469 0.0001553692118 -0.1580604617
+8.5 1 641.6199025 0.2731219637 0.00009119724775 -0.0261686794
+9 1 1030.914723 0.2453117866 0.00005363701638 0.1043145752
+9.5 1 1658.453078 0.1612644308 0.00003160203411 0.2031798994
+10 1 2670.988304 0.04347274617 0.00001864877345 0.2490154242
+0.5 2 0.03190614918 0.03060402346 7.550183551 -5.441370837
+1 2 0.1357476698 0.1149034849 1.624838899 -1.650682607
+1.5 2 0.3378346183 0.2320876721 0.5836559633 -0.9321937598
+2 2 0.6889484477 0.3528340286 0.2537597546 -0.6174081042
+2.5 2 1.276466148 0.4460590584 0.1214602063 -0.3813358492
+3 2 2.245212441 0.4860912606 0.06151045847 -0.1604003935
+3.5 2 3.832012048 0.4586291842 0.0323071217 0.04537143773
+4 2 6.422189375 0.3641281459 0.01740142553 0.2159035946
+4.5 2 10.6415173 0.2178489837 0.009545677203 0.3284815967
+5 2 17.50561497 0.04656511628 0.005308943712 0.3676628826
+5.5 2 28.66258529 -0.1173154816 0.002984370024 0.3308412333
+6 2 46.78709472 -0.24287321 0.001691967567 0.2298579025
+6.5 2 76.22054724 -0.3074303906 0.0009658992748 0.08890665832
+7 2 124.0113105 -0.3014172201 0.0005545621667 -0.06052660947
+7.5 2 201.6054807 -0.2302734105 0.0003199235871 -0.1864142223
+8 2 327.5958315 -0.1129917204 0.0001853130082 -0.2630366048
+8.5 2 532.1925382 0.02232473961 0.0001077157423 -0.2763624417
+9 2 864.496194 0.1448473415 0.00006280064993 -0.2267556816
+9.5 2 1404.332974 0.2278791542 0.00003671094477 -0.1284359105
+10 2 2281.518968 0.2546303137 0.00002150981701 -0.005868082442
+0.5 3 0.002645111969 0.002563729995 62.05790953 -42.0594943
+1 3 0.02216842492 0.01956335398 7.101262825 -5.821517606
+1.5 3 0.08077411302 0.06096395114 1.833803702 -2.073541399
+2 3 0.2127399592 0.1289432495 0.6473853909 -1.127783777
+2.5 3 0.4743704088 0.216600391 0.2682271464 -0.7560554968
+3 3 0.9597536295 0.3090627223 0.1221703758 -0.5385416161
+3.5 3 1.826392582 0.3867701117 0.05916181773 -0.3583353462
+4 3 3.337275778 0.4301714739 0.02988492442 -0.182022116
+4.5 3 5.930096266 0.424703973 0.01556314131 -0.009013681594
+5 3 10.33115017 0.3648312306 0.008291768415 0.1462671627
+5.5 3 17.74264804 0.2561178651 0.004496019935 0.264370045
+6 3 30.1505403 0.1147683848 0.002471898096 0.328248946
+6.5 3 50.83005863 -0.03534663129 0.001374293952 0.3288030637
+7 3 85.17548684 -0.167555588 0.0007710751536 0.268080603
+7.5 3 142.0614424 -0.2580609132 0.0004359233032 0.1597075919
+8 3 236.075221 -0.2911322071 0.0002480257159 0.02654215932
+8.5 3 391.1763552 -0.2626162039 0.0001418870088 -0.1038842343
+9 3 646.6941919 -0.1809351903 0.00008154841635 -0.2050948781
+9.5 3 1067.154983 -0.06531531322 0.00004705927401 -0.2572581775
+10 3 1758.380717 0.05837937931 0.00002725270026 -0.2513626572