001/* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 018package org.apache.commons.math3.special; 019 020import org.apache.commons.math3.analysis.UnivariateFunction; 021import org.apache.commons.math3.exception.ConvergenceException; 022import org.apache.commons.math3.exception.MathIllegalArgumentException; 023import org.apache.commons.math3.exception.util.LocalizedFormats; 024import org.apache.commons.math3.util.FastMath; 025import org.apache.commons.math3.util.MathArrays; 026 027/** 028 * This class provides computation methods related to Bessel 029 * functions of the first kind. Detailed descriptions of these functions are 030 * available in <a 031 * href="http://en.wikipedia.org/wiki/Bessel_function">Wikipedia</a>, <a 032 * href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun">Abrabowitz and 033 * Stegun</a> (Ch. 9-11), and <a href="http://dlmf.nist.gov/">DLMF</a> (Ch. 10). 034 * <p> 035 * This implementation is based on the rjbesl Fortran routine at 036 * <a href="http://www.netlib.org/specfun/rjbesl">Netlib</a>.</p> 037 * <p> 038 * From the Fortran code: </p> 039 * <p> 040 * This program is based on a program written by David J. Sookne (2) that 041 * computes values of the Bessel functions J or I of real argument and integer 042 * order. Modifications include the restriction of the computation to the J 043 * Bessel function of non-negative real argument, the extension of the 044 * computation to arbitrary positive order, and the elimination of most 045 * underflow.</p> 046 * <p> 047 * References:</p> 048 * <ul> 049 * <li>"A Note on Backward Recurrence Algorithms," Olver, F. W. J., and Sookne, 050 * D. J., Math. Comp. 26, 1972, pp 941-947.</li> 051 * <li>"Bessel Functions of Real Argument and Integer Order," Sookne, D. J., NBS 052 * Jour. of Res. B. 77B, 1973, pp 125-132.</li> 053 * </ul> </p> 054 * @since 3.4 055 */ 056public class BesselJ 057 implements UnivariateFunction { 058 059 // --------------------------------------------------------------------- 060 // Mathematical constants 061 // --------------------------------------------------------------------- 062 063 /** -2 / pi */ 064 private static final double PI2 = 0.636619772367581343075535; 065 066 /** first few significant digits of 2pi */ 067 private static final double TOWPI1 = 6.28125; 068 069 /** 2pi - TWOPI1 to working precision */ 070 private static final double TWOPI2 = 1.935307179586476925286767e-3; 071 072 /** TOWPI1 + TWOPI2 */ 073 private static final double TWOPI = TOWPI1 + TWOPI2; 074 075 // --------------------------------------------------------------------- 076 // Machine-dependent parameters 077 // --------------------------------------------------------------------- 078 079 /** 080 * 10.0^K, where K is the largest integer such that ENTEN is 081 * machine-representable in working precision 082 */ 083 private static final double ENTEN = 1.0e308; 084 085 /** 086 * Decimal significance desired. Should be set to (INT(log_{10}(2) * (it)+1)). 087 * Setting NSIG lower will result in decreased accuracy while setting 088 * NSIG higher will increase CPU time without increasing accuracy. 089 * The truncation error is limited to a relative error of 090 * T=.5(10^(-NSIG)). 091 */ 092 private static final double ENSIG = 1.0e16; 093 094 /** 10.0 ** (-K) for the smallest integer K such that K >= NSIG/4 */ 095 private static final double RTNSIG = 1.0e-4; 096 097 /** Smallest ABS(X) such that X/4 does not underflow */ 098 private static final double ENMTEN = 8.90e-308; 099 100 /** Minimum acceptable value for x */ 101 private static final double X_MIN = 0.0; 102 103 /** 104 * Upper limit on the magnitude of x. If abs(x) = n, then at least 105 * n iterations of the backward recursion will be executed. The value of 106 * 10.0 ** 4 is used on every machine. 107 */ 108 private static final double X_MAX = 1.0e4; 109 110 /** First 25 factorials as doubles */ 111 private static final double[] FACT = { 112 1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0, 113 3628800.0, 39916800.0, 479001600.0, 6227020800.0, 87178291200.0, 114 1.307674368e12, 2.0922789888e13, 3.55687428096e14, 6.402373705728e15, 115 1.21645100408832e17, 2.43290200817664e18, 5.109094217170944e19, 116 1.12400072777760768e21, 2.585201673888497664e22, 117 6.2044840173323943936e23 118 }; 119 120 /** Order of the function computed when {@link #value(double)} is used */ 121 private final double order; 122 123 /** 124 * Create a new BesselJ with the given order. 125 * 126 * @param order order of the function computed when using {@link #value(double)}. 127 */ 128 public BesselJ(double order) { 129 this.order = order; 130 } 131 132 /** 133 * Returns the value of the constructed Bessel function of the first kind, 134 * for the passed argument. 135 * 136 * @param x Argument 137 * @return Value of the Bessel function at x 138 * @throws MathIllegalArgumentException if {@code x} is too large relative to {@code order} 139 * @throws ConvergenceException if the algorithm fails to converge 140 */ 141 public double value(double x) 142 throws MathIllegalArgumentException, ConvergenceException { 143 return BesselJ.value(order, x); 144 } 145 146 /** 147 * Returns the first Bessel function, \(J_{order}(x)\). 148 * 149 * @param order Order of the Bessel function 150 * @param x Argument 151 * @return Value of the Bessel function of the first kind, \(J_{order}(x)\) 152 * @throws MathIllegalArgumentException if {@code x} is too large relative to {@code order} 153 * @throws ConvergenceException if the algorithm fails to converge 154 */ 155 public static double value(double order, double x) 156 throws MathIllegalArgumentException, ConvergenceException { 157 final int n = (int) order; 158 final double alpha = order - n; 159 final int nb = n + 1; 160 final BesselJResult res = rjBesl(x, alpha, nb); 161 162 if (res.nVals >= nb) { 163 return res.vals[n]; 164 } else if (res.nVals < 0) { 165 throw new MathIllegalArgumentException(LocalizedFormats.BESSEL_FUNCTION_BAD_ARGUMENT,order, x); 166 } else if (FastMath.abs(res.vals[res.nVals - 1]) < 1e-100) { 167 return res.vals[n]; // underflow; return value (will be zero) 168 } 169 throw new ConvergenceException(LocalizedFormats.BESSEL_FUNCTION_FAILED_CONVERGENCE, order, x); 170 } 171 172 /** 173 * Encapsulates the results returned by {@link BesselJ#rjBesl(double, double, int)}. 174 * <p> 175 * {@link #getVals()} returns the computed function values. 176 * {@link #getnVals()} is the number of values among those returned by {@link #getnVals()} 177 * that can be considered accurate. 178 * </p><p> 179 * <ul> 180 * <li>nVals < 0: An argument is out of range. For example, nb <= 0, alpha 181 * < 0 or > 1, or x is too large. In this case, b(0) is set to zero, the 182 * remainder of the b-vector is not calculated, and nVals is set to 183 * MIN(nb,0) - 1 so that nVals != nb.</li> 184 * <li>nb > nVals > 0: Not all requested function values could be calculated 185 * accurately. This usually occurs because nb is much larger than abs(x). In 186 * this case, b(n) is calculated to the desired accuracy for n < nVals, but 187 * precision is lost for nVals < n <= nb. If b(n) does not vanish for n > 188 * nVals (because it is too small to be represented), and b(n)/b(nVals) = 189 * \(10^{-k}\), then only the first NSIG-k significant figures of b(n) can be 190 * trusted.</li></ul></p> 191 */ 192 public static class BesselJResult { 193 194 /** Bessel function values */ 195 private final double[] vals; 196 197 /** Valid value count */ 198 private final int nVals; 199 200 /** 201 * Create a new BesselJResult with the given values and valid value count. 202 * 203 * @param b values 204 * @param n count of valid values 205 */ 206 public BesselJResult(double[] b, int n) { 207 vals = MathArrays.copyOf(b, b.length); 208 nVals = n; 209 } 210 211 /** 212 * @return the computed function values 213 */ 214 public double[] getVals() { 215 return MathArrays.copyOf(vals, vals.length); 216 } 217 218 /** 219 * @return the number of valid function values (normally the same as the 220 * length of the array returned by {@link #getnVals()}) 221 */ 222 public int getnVals() { 223 return nVals; 224 } 225 } 226 227 /** 228 * Calculates Bessel functions \(J_{n+alpha}(x)\) for 229 * non-negative argument x, and non-negative order n + alpha. 230 * <p> 231 * Before using the output vector, the user should check that 232 * nVals = nb, i.e., all orders have been calculated to the desired accuracy. 233 * See BesselResult class javadoc for details on return values. 234 * </p> 235 * @param x non-negative real argument for which J's are to be calculated 236 * @param alpha fractional part of order for which J's or exponentially 237 * scaled J's (\(J\cdot e^{x}\)) are to be calculated. 0 <= alpha < 1.0. 238 * @param nb integer number of functions to be calculated, nb > 0. The first 239 * function calculated is of order alpha, and the last is of order 240 * nb - 1 + alpha. 241 * @return BesselJResult a vector of the functions 242 * \(J_{alpha}(x)\) through \(J_{nb-1+alpha}(x)\), or the corresponding exponentially 243 * scaled functions and an integer output variable indicating possible errors 244 */ 245 public static BesselJResult rjBesl(double x, double alpha, int nb) { 246 final double[] b = new double[nb]; 247 248 int ncalc = 0; 249 double alpem = 0; 250 double alp2em = 0; 251 252 // --------------------------------------------------------------------- 253 // Check for out of range arguments. 254 // --------------------------------------------------------------------- 255 final int magx = (int) x; 256 if ((nb > 0) && (x >= X_MIN) && (x <= X_MAX) && (alpha >= 0) && 257 (alpha < 1)) { 258 // --------------------------------------------------------------------- 259 // Initialize result array to zero. 260 // --------------------------------------------------------------------- 261 ncalc = nb; 262 for (int i = 0; i < nb; ++i) { 263 b[i] = 0; 264 } 265 266 // --------------------------------------------------------------------- 267 // Branch to use 2-term ascending series for small X and asymptotic 268 // form for large X when NB is not too large. 269 // --------------------------------------------------------------------- 270 double tempa; 271 double tempb; 272 double tempc; 273 double tover; 274 if (x < RTNSIG) { 275 // --------------------------------------------------------------------- 276 // Two-term ascending series for small X. 277 // --------------------------------------------------------------------- 278 tempa = 1; 279 alpem = 1 + alpha; 280 double halfx = 0; 281 if (x > ENMTEN) { 282 halfx = 0.5 * x; 283 } 284 if (alpha != 0) { 285 tempa = FastMath.pow(halfx, alpha) / 286 (alpha * Gamma.gamma(alpha)); 287 } 288 tempb = 0; 289 if (x + 1 > 1) { 290 tempb = -halfx * halfx; 291 } 292 b[0] = tempa + (tempa * tempb / alpem); 293 if ((x != 0) && (b[0] == 0)) { 294 ncalc = 0; 295 } 296 if (nb != 1) { 297 if (x <= 0) { 298 for (int n = 1; n < nb; ++n) { 299 b[n] = 0; 300 } 301 } else { 302 // --------------------------------------------------------------------- 303 // Calculate higher order functions. 304 // --------------------------------------------------------------------- 305 tempc = halfx; 306 tover = tempb != 0 ? ENMTEN / tempb : 2 * ENMTEN / x; 307 for (int n = 1; n < nb; ++n) { 308 tempa /= alpem; 309 alpem += 1; 310 tempa *= tempc; 311 if (tempa <= tover * alpem) { 312 tempa = 0; 313 } 314 b[n] = tempa + (tempa * tempb / alpem); 315 if ((b[n] == 0) && (ncalc > n)) { 316 ncalc = n; 317 } 318 } 319 } 320 } 321 } else if ((x > 25.0) && (nb <= magx + 1)) { 322 // --------------------------------------------------------------------- 323 // Asymptotic series for X > 25 324 // --------------------------------------------------------------------- 325 final double xc = FastMath.sqrt(PI2 / x); 326 final double mul = 0.125 / x; 327 final double xin = mul * mul; 328 int m = 0; 329 if (x >= 130.0) { 330 m = 4; 331 } else if (x >= 35.0) { 332 m = 8; 333 } else { 334 m = 11; 335 } 336 337 final double xm = 4.0 * m; 338 // --------------------------------------------------------------------- 339 // Argument reduction for SIN and COS routines. 340 // --------------------------------------------------------------------- 341 double t = (double) ((int) ((x / TWOPI) + 0.5)); 342 final double z = x - t * TOWPI1 - t * TWOPI2 - (alpha + 0.5) / PI2; 343 double vsin = FastMath.sin(z); 344 double vcos = FastMath.cos(z); 345 double gnu = 2 * alpha; 346 double capq; 347 double capp; 348 double s; 349 double t1; 350 double xk; 351 for (int i = 1; i <= 2; i++) { 352 s = (xm - 1 - gnu) * (xm - 1 + gnu) * xin * 0.5; 353 t = (gnu - (xm - 3.0)) * (gnu + (xm - 3.0)); 354 capp = (s * t) / FACT[2 * m]; 355 t1 = (gnu - (xm + 1)) * (gnu + (xm + 1)); 356 capq = (s * t1) / FACT[2 * m + 1]; 357 xk = xm; 358 int k = 2 * m; 359 t1 = t; 360 361 for (int j = 2; j <= m; j++) { 362 xk -= 4.0; 363 s = (xk - 1 - gnu) * (xk - 1 + gnu); 364 t = (gnu - (xk - 3.0)) * (gnu + (xk - 3.0)); 365 capp = (capp + 1 / FACT[k - 2]) * s * t * xin; 366 capq = (capq + 1 / FACT[k - 1]) * s * t1 * xin; 367 k -= 2; 368 t1 = t; 369 } 370 371 capp += 1; 372 capq = (capq + 1) * ((gnu * gnu) - 1) * (0.125 / x); 373 b[i - 1] = xc * (capp * vcos - capq * vsin); 374 if (nb == 1) { 375 return new BesselJResult(MathArrays.copyOf(b, b.length), 376 ncalc); 377 } 378 t = vsin; 379 vsin = -vcos; 380 vcos = t; 381 gnu += 2.0; 382 } 383 384 // --------------------------------------------------------------------- 385 // If NB > 2, compute J(X,ORDER+I) I = 2, NB-1 386 // --------------------------------------------------------------------- 387 if (nb > 2) { 388 gnu = 2 * alpha + 2.0; 389 for (int j = 2; j < nb; ++j) { 390 b[j] = gnu * b[j - 1] / x - b[j - 2]; 391 gnu += 2.0; 392 } 393 } 394 } else { 395 // --------------------------------------------------------------------- 396 // Use recurrence to generate results. First initialize the 397 // calculation of P*S. 398 // --------------------------------------------------------------------- 399 final int nbmx = nb - magx; 400 int n = magx + 1; 401 int nstart = 0; 402 int nend = 0; 403 double en = 2 * (n + alpha); 404 double plast = 1; 405 double p = en / x; 406 double pold; 407 // --------------------------------------------------------------------- 408 // Calculate general significance test. 409 // --------------------------------------------------------------------- 410 double test = 2 * ENSIG; 411 boolean readyToInitialize = false; 412 if (nbmx >= 3) { 413 // --------------------------------------------------------------------- 414 // Calculate P*S until N = NB-1. Check for possible 415 // overflow. 416 // --------------------------------------------------------------------- 417 tover = ENTEN / ENSIG; 418 nstart = magx + 2; 419 nend = nb - 1; 420 en = 2 * (nstart - 1 + alpha); 421 double psave; 422 double psavel; 423 for (int k = nstart; k <= nend; k++) { 424 n = k; 425 en += 2.0; 426 pold = plast; 427 plast = p; 428 p = (en * plast / x) - pold; 429 if (p > tover) { 430 // --------------------------------------------------------------------- 431 // To avoid overflow, divide P*S by TOVER. Calculate 432 // P*S until 433 // ABS(P) > 1. 434 // --------------------------------------------------------------------- 435 tover = ENTEN; 436 p /= tover; 437 plast /= tover; 438 psave = p; 439 psavel = plast; 440 nstart = n + 1; 441 do { 442 n += 1; 443 en += 2.0; 444 pold = plast; 445 plast = p; 446 p = (en * plast / x) - pold; 447 } while (p <= 1); 448 tempb = en / x; 449 // --------------------------------------------------------------------- 450 // Calculate backward test and find NCALC, the 451 // highest N such that 452 // the test is passed. 453 // --------------------------------------------------------------------- 454 test = pold * plast * (0.5 - 0.5 / (tempb * tempb)); 455 test /= ENSIG; 456 p = plast * tover; 457 n -= 1; 458 en -= 2.0; 459 nend = FastMath.min(nb, n); 460 for (int l = nstart; l <= nend; l++) { 461 pold = psavel; 462 psavel = psave; 463 psave = (en * psavel / x) - pold; 464 if (psave * psavel > test) { 465 ncalc = l - 1; 466 readyToInitialize = true; 467 break; 468 } 469 } 470 ncalc = nend; 471 readyToInitialize = true; 472 break; 473 } 474 } 475 if (!readyToInitialize) { 476 n = nend; 477 en = 2 * (n + alpha); 478 // --------------------------------------------------------------------- 479 // Calculate special significance test for NBMX > 2. 480 // --------------------------------------------------------------------- 481 test = FastMath.max(test, FastMath.sqrt(plast * ENSIG) * 482 FastMath.sqrt(2 * p)); 483 } 484 } 485 // --------------------------------------------------------------------- 486 // Calculate P*S until significance test passes. 487 // --------------------------------------------------------------------- 488 if (!readyToInitialize) { 489 do { 490 n += 1; 491 en += 2.0; 492 pold = plast; 493 plast = p; 494 p = (en * plast / x) - pold; 495 } while (p < test); 496 } 497 // --------------------------------------------------------------------- 498 // Initialize the backward recursion and the normalization sum. 499 // --------------------------------------------------------------------- 500 n += 1; 501 en += 2.0; 502 tempb = 0; 503 tempa = 1 / p; 504 int m = (2 * n) - 4 * (n / 2); 505 double sum = 0; 506 double em = (double) (n / 2); 507 alpem = em - 1 + alpha; 508 alp2em = 2 * em + alpha; 509 if (m != 0) { 510 sum = tempa * alpem * alp2em / em; 511 } 512 nend = n - nb; 513 514 boolean readyToNormalize = false; 515 boolean calculatedB0 = false; 516 517 // --------------------------------------------------------------------- 518 // Recur backward via difference equation, calculating (but not 519 // storing) B(N), until N = NB. 520 // --------------------------------------------------------------------- 521 for (int l = 1; l <= nend; l++) { 522 n -= 1; 523 en -= 2.0; 524 tempc = tempb; 525 tempb = tempa; 526 tempa = (en * tempb / x) - tempc; 527 m = 2 - m; 528 if (m != 0) { 529 em -= 1; 530 alp2em = 2 * em + alpha; 531 if (n == 1) { 532 break; 533 } 534 alpem = em - 1 + alpha; 535 if (alpem == 0) { 536 alpem = 1; 537 } 538 sum = (sum + tempa * alp2em) * alpem / em; 539 } 540 } 541 542 // --------------------------------------------------------------------- 543 // Store B(NB). 544 // --------------------------------------------------------------------- 545 b[n - 1] = tempa; 546 if (nend >= 0) { 547 if (nb <= 1) { 548 alp2em = alpha; 549 if (alpha + 1 == 1) { 550 alp2em = 1; 551 } 552 sum += b[0] * alp2em; 553 readyToNormalize = true; 554 } else { 555 // --------------------------------------------------------------------- 556 // Calculate and store B(NB-1). 557 // --------------------------------------------------------------------- 558 n -= 1; 559 en -= 2.0; 560 b[n - 1] = (en * tempa / x) - tempb; 561 if (n == 1) { 562 calculatedB0 = true; 563 } else { 564 m = 2 - m; 565 if (m != 0) { 566 em -= 1; 567 alp2em = 2 * em + alpha; 568 alpem = em - 1 + alpha; 569 if (alpem == 0) { 570 alpem = 1; 571 } 572 573 sum = (sum + (b[n - 1] * alp2em)) * alpem / em; 574 } 575 } 576 } 577 } 578 if (!readyToNormalize && !calculatedB0) { 579 nend = n - 2; 580 if (nend != 0) { 581 // --------------------------------------------------------------------- 582 // Calculate via difference equation and store B(N), 583 // until N = 2. 584 // --------------------------------------------------------------------- 585 586 for (int l = 1; l <= nend; l++) { 587 n -= 1; 588 en -= 2.0; 589 b[n - 1] = (en * b[n] / x) - b[n + 1]; 590 m = 2 - m; 591 if (m != 0) { 592 em -= 1; 593 alp2em = 2 * em + alpha; 594 alpem = em - 1 + alpha; 595 if (alpem == 0) { 596 alpem = 1; 597 } 598 599 sum = (sum + b[n - 1] * alp2em) * alpem / em; 600 } 601 } 602 } 603 } 604 // --------------------------------------------------------------------- 605 // Calculate b[0] 606 // --------------------------------------------------------------------- 607 if (!readyToNormalize) { 608 if (!calculatedB0) { 609 b[0] = 2.0 * (alpha + 1) * b[1] / x - b[2]; 610 } 611 em -= 1; 612 alp2em = 2 * em + alpha; 613 if (alp2em == 0) { 614 alp2em = 1; 615 } 616 sum += b[0] * alp2em; 617 } 618 // --------------------------------------------------------------------- 619 // Normalize. Divide all B(N) by sum. 620 // --------------------------------------------------------------------- 621 622 if (FastMath.abs(alpha) > 1e-16) { 623 sum *= Gamma.gamma(alpha) * FastMath.pow(x * 0.5, -alpha); 624 } 625 tempa = ENMTEN; 626 if (sum > 1) { 627 tempa *= sum; 628 } 629 630 for (n = 0; n < nb; n++) { 631 if (FastMath.abs(b[n]) < tempa) { 632 b[n] = 0; 633 } 634 b[n] /= sum; 635 } 636 } 637 // --------------------------------------------------------------------- 638 // Error return -- X, NB, or ALPHA is out of range. 639 // --------------------------------------------------------------------- 640 } else { 641 if (b.length > 0) { 642 b[0] = 0; 643 } 644 ncalc = FastMath.min(nb, 0) - 1; 645 } 646 return new BesselJResult(MathArrays.copyOf(b, b.length), ncalc); 647 } 648}