j1.js (3610B)
1 /** 2 * @license Apache-2.0 3 * 4 * Copyright (c) 2018 The Stdlib Authors. 5 * 6 * Licensed under the Apache License, Version 2.0 (the "License"); 7 * you may not use this file except in compliance with the License. 8 * You may obtain a copy of the License at 9 * 10 * http://www.apache.org/licenses/LICENSE-2.0 11 * 12 * Unless required by applicable law or agreed to in writing, software 13 * distributed under the License is distributed on an "AS IS" BASIS, 14 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 15 * See the License for the specific language governing permissions and 16 * limitations under the License. 17 * 18 * 19 * ## Notice 20 * 21 * The original C++ code and copyright notice are from the [Boost library]{@link https://github.com/boostorg/math/blob/develop/include/boost/math/special_functions/detail/bessel_j1.hpp}. The implementation has been modified for JavaScript. 22 * 23 * ```text 24 * Copyright Xiaogang Zhang, 2006. 25 * 26 * Use, modification and distribution are subject to the 27 * Boost Software License, Version 1.0. (See accompanying file 28 * LICENSE or copy at http://www.boost.org/LICENSE_1_0.txt) 29 * ``` 30 */ 31 32 'use strict'; 33 34 // MODULES // 35 36 var sqrt = require( './../../../../base/special/sqrt' ); 37 var abs = require( './../../../../base/special/abs' ); 38 var sincos = require( './../../../../base/special/sincos' ); 39 var PINF = require( '@stdlib/constants/float64/pinf' ); 40 var SQRT_PI = require( '@stdlib/constants/float64/sqrt-pi' ); 41 var poly1 = require( './rational_p1q1.js' ); 42 var poly2 = require( './rational_p2q2.js' ); 43 var polyC = require( './rational_pcqc.js' ); 44 var polyS = require( './rational_psqs.js' ); 45 46 47 // VARIABLES // 48 49 var x1 = 3.8317059702075123156e+00; 50 var x2 = 7.0155866698156187535e+00; 51 var x11 = 9.810e+02; 52 var x12 = -3.2527979248768438556e-04; 53 var x21 = 1.7960e+03; 54 var x22 = -3.8330184381246462950e-05; 55 56 // `sincos` workspace: 57 var sc = [ 0.0, 0.0 ]; // WARNING: not thread safe 58 59 60 // MAIN // 61 62 /** 63 * Computes the Bessel function of the first kind of order one. 64 * 65 * ## Notes 66 * 67 * - Accuracy for subnormal `x` is very poor. Full accuracy is achieved at `1.0e-308` but trends progressively to zero at `5e-324`. This suggests that underflow (or overflow, perhaps due to a reciprocal) is effectively cutting off digits of precision until the computation loses all accuracy at `5e-324`. 68 * 69 * @param {number} x - input value 70 * @returns {number} evaluated Bessel function 71 * 72 * @example 73 * var v = j1( 0.0 ); 74 * // returns 0.0 75 * 76 * v = j1( 1.0 ); 77 * // returns ~0.440 78 * 79 * v = j1( Infinity ); 80 * // returns 0.0 81 * 82 * v = j1( -Infinity ); 83 * // returns 0.0 84 * 85 * v = j1( NaN ); 86 * // returns NaN 87 */ 88 function j1( x ) { 89 var value; 90 var rc; 91 var rs; 92 var y2; 93 var r; 94 var y; 95 var f; 96 var w; 97 98 w = abs( x ); 99 if ( x === 0.0 ) { 100 return 0.0; 101 } 102 if ( w === PINF ) { 103 return 0.0; 104 } 105 if ( w <= 4.0 ) { 106 y = x * x; 107 r = poly1( y ); 108 f = w * ( w+x1 ) * ( ( w - (x11/256.0) ) - x12 ); 109 value = f * r; 110 } else if ( w <= 8.0 ) { 111 y = x * x; 112 r = poly2( y ); 113 f = w * ( w+x2 ) * ( ( w - (x21/256.0) ) - x22 ); 114 value = f * r; 115 } else { 116 y = 8.0 / w; 117 y2 = y * y; 118 rc = polyC( y2 ); 119 rs = polyS( y2 ); 120 f = 1.0 / ( sqrt( w ) * SQRT_PI ); 121 122 /* 123 * What follows is really just: 124 * 125 * ``` 126 * z = w - 0.75 * pi; 127 * value = f * ( rc * cos( z ) - y * rs * sin( z ) ); 128 * ``` 129 * 130 * but using the sin/cos addition rules plus constants for the values of sin/cos of `3π/4` which then cancel out with corresponding terms in "f". 131 */ 132 sincos( sc, w ); 133 value = f * ( ( rc * (sc[0]-sc[1]) ) + ( (y*rs) * (sc[0]+sc[1]) ) ); 134 } 135 if ( x < 0.0 ) { 136 value *= -1.0; 137 } 138 return value; 139 } 140 141 142 // EXPORTS // 143 144 module.exports = j1;