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.math4.legacy.special; 019 020import java.util.Arrays; 021import org.apache.commons.numbers.gamma.Gamma; 022import org.apache.commons.math4.legacy.analysis.UnivariateFunction; 023import org.apache.commons.math4.legacy.exception.ConvergenceException; 024import org.apache.commons.math4.legacy.exception.MathIllegalArgumentException; 025import org.apache.commons.math4.legacy.exception.util.LocalizedFormats; 026import org.apache.commons.math4.core.jdkmath.JdkMath; 027 028/** 029 * This class provides computation methods related to Bessel 030 * functions of the first kind. Detailed descriptions of these functions are 031 * available in <a 032 * href="http://en.wikipedia.org/wiki/Bessel_function">Wikipedia</a>, <a 033 * href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun">Abrabowitz and 034 * Stegun</a> (Ch. 9-11), and <a href="http://dlmf.nist.gov/">DLMF</a> (Ch. 10). 035 * <p> 036 * This implementation is based on the rjbesl Fortran routine at 037 * <a href="http://www.netlib.org/specfun/rjbesl">Netlib</a>.</p> 038 * <p> 039 * From the Fortran code: </p> 040 * <p> 041 * This program is based on a program written by David J. Sookne (2) that 042 * computes values of the Bessel functions J or I of real argument and integer 043 * order. Modifications include the restriction of the computation to the J 044 * Bessel function of non-negative real argument, the extension of the 045 * computation to arbitrary positive order, and the elimination of most 046 * underflow.</p> 047 * <p> 048 * References:</p> 049 * <ul> 050 * <li>"A Note on Backward Recurrence Algorithms," Olver, F. W. J., and Sookne, 051 * D. J., Math. Comp. 26, 1972, pp 941-947.</li> 052 * <li>"Bessel Functions of Real Argument and Integer Order," Sookne, D. J., NBS 053 * Jour. of Res. B. 77B, 1973, pp 125-132.</li> 054 * </ul> 055 * @since 3.4 056 */ 057public class BesselJ 058 implements UnivariateFunction { 059 060 // --------------------------------------------------------------------- 061 // Mathematical constants 062 // --------------------------------------------------------------------- 063 064 /** -2 / pi. */ 065 private static final double PI2 = 0.636619772367581343075535; 066 067 /** first few significant digits of 2pi. */ 068 private static final double TOWPI1 = 6.28125; 069 070 /** 2pi - TWOPI1 to working precision. */ 071 private static final double TWOPI2 = 1.935307179586476925286767e-3; 072 073 /** TOWPI1 + TWOPI2. */ 074 private static final double TWOPI = TOWPI1 + TWOPI2; 075 076 // --------------------------------------------------------------------- 077 // Machine-dependent parameters 078 // --------------------------------------------------------------------- 079 080 /** 081 * 10.0^K, where K is the largest integer such that ENTEN is 082 * machine-representable in working precision. 083 */ 084 private static final double ENTEN = 1.0e308; 085 086 /** 087 * Decimal significance desired. Should be set to (INT(log_{10}(2) * (it)+1)). 088 * Setting NSIG lower will result in decreased accuracy while setting 089 * NSIG higher will increase CPU time without increasing accuracy. 090 * The truncation error is limited to a relative error of 091 * T=.5(10^(-NSIG)). 092 */ 093 private static final double ENSIG = 1.0e16; 094 095 /** 10.0 ** (-K) for the smallest integer K such that K >= NSIG/4. */ 096 private static final double RTNSIG = 1.0e-4; 097 098 /** Smallest ABS(X) such that X/4 does not underflow. */ 099 private static final double ENMTEN = 8.90e-308; 100 101 /** Minimum acceptable value for x. */ 102 private static final double X_MIN = 0.0; 103 104 /** 105 * Upper limit on the magnitude of x. If abs(x) = n, then at least 106 * n iterations of the backward recursion will be executed. The value of 107 * 10.0 ** 4 is used on every machine. 108 */ 109 private static final double X_MAX = 1.0e4; 110 111 /** First 25 factorials as doubles. */ 112 private static final double[] FACT = { 113 1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0, 114 3628800.0, 39916800.0, 479001600.0, 6227020800.0, 87178291200.0, 115 1.307674368e12, 2.0922789888e13, 3.55687428096e14, 6.402373705728e15, 116 1.21645100408832e17, 2.43290200817664e18, 5.109094217170944e19, 117 1.12400072777760768e21, 2.585201673888497664e22, 118 6.2044840173323943936e23 119 }; 120 121 /** Order of the function computed when {@link #value(double)} is used. */ 122 private final double order; 123 124 /** 125 * Create a new BesselJ with the given order. 126 * 127 * @param order order of the function computed when using {@link #value(double)}. 128 */ 129 public BesselJ(double order) { 130 this.order = order; 131 } 132 133 /** 134 * Returns the value of the constructed Bessel function of the first kind, 135 * for the passed argument. 136 * 137 * @param x Argument 138 * @return Value of the Bessel function at x 139 * @throws MathIllegalArgumentException if {@code x} is too large relative to {@code order} 140 * @throws ConvergenceException if the algorithm fails to converge 141 */ 142 @Override 143 public double value(double x) 144 throws MathIllegalArgumentException, ConvergenceException { 145 return BesselJ.value(order, x); 146 } 147 148 /** 149 * Returns the first Bessel function, \(J_{order}(x)\). 150 * 151 * @param order Order of the Bessel function 152 * @param x Argument 153 * @return Value of the Bessel function of the first kind, \(J_{order}(x)\) 154 * @throws MathIllegalArgumentException if {@code x} is too large relative to {@code order} 155 * @throws ConvergenceException if the algorithm fails to converge 156 */ 157 public static double value(double order, double x) 158 throws MathIllegalArgumentException, ConvergenceException { 159 final int n = (int) order; 160 final double alpha = order - n; 161 final int nb = n + 1; 162 final BesselJResult res = rjBesl(x, alpha, nb); 163 164 if (res.nVals >= nb) { 165 return res.vals[n]; 166 } else if (res.nVals < 0) { 167 throw new MathIllegalArgumentException(LocalizedFormats.BESSEL_FUNCTION_BAD_ARGUMENT,order, x); 168 } else if (JdkMath.abs(res.vals[res.nVals - 1]) < 1e-100) { 169 return res.vals[n]; // underflow; return value (will be zero) 170 } 171 throw new ConvergenceException(LocalizedFormats.BESSEL_FUNCTION_FAILED_CONVERGENCE, order, x); 172 } 173 174 /** 175 * Encapsulates the results returned by {@link BesselJ#rjBesl(double, double, int)}. 176 * <p> 177 * {@link #getVals()} returns the computed function values. 178 * {@link #getnVals()} is the number of values among those returned by {@link #getnVals()} 179 * that can be considered accurate. 180 * </p><ul> 181 * <li>{@code nVals < 0}: An argument is out of range. For example, {@code nb <= 0}, 182 * {@code alpha < 0 or > 1}, or x is too large. In this case, b(0) is set to zero, the 183 * remainder of the b-vector is not calculated, and nVals is set to 184 * MIN(nb,0) - 1 so that nVals != nb.</li> 185 * <li>{@code nb > nVals > 0}: Not all requested function values could be calculated 186 * accurately. This usually occurs because nb is much larger than abs(x). In 187 * this case, b(n) is calculated to the desired accuracy for {@code n < nVals}, but 188 * precision is lost for {@code nVals < n <= nb}. If b(n) does not vanish for 189 * {@code n > nVals} (because it is too small to be represented), and b(n)/b(nVals) = 190 * \(10^{-k}\), then only the first NSIG-k significant figures of b(n) can be 191 * trusted.</li></ul> 192 */ 193 public static class BesselJResult { 194 195 /** Bessel function values. */ 196 private final double[] vals; 197 198 /** Valid value count. */ 199 private final int nVals; 200 201 /** 202 * Create a new BesselJResult with the given values and valid value count. 203 * 204 * @param b values 205 * @param n count of valid values 206 */ 207 public BesselJResult(double[] b, int n) { 208 vals = Arrays.copyOf(b, b.length); 209 nVals = n; 210 } 211 212 /** 213 * @return the computed function values 214 */ 215 public double[] getVals() { 216 return Arrays.copyOf(vals, vals.length); 217 } 218 219 /** 220 * @return the number of valid function values (normally the same as the 221 * length of the array returned by {@link #getnVals()}) 222 */ 223 public int getnVals() { 224 return nVals; 225 } 226 } 227 228 /** 229 * Calculates Bessel functions \(J_{n+alpha}(x)\) for 230 * non-negative argument x, and non-negative order n + alpha. 231 * <p> 232 * Before using the output vector, the user should check that 233 * nVals = nb, i.e., all orders have been calculated to the desired accuracy. 234 * See BesselResult class javadoc for details on return values. 235 * </p> 236 * @param x non-negative real argument for which J's are to be calculated 237 * @param alpha fractional part of order for which J's or exponentially 238 * scaled J's (\(J\cdot e^{x}\)) are to be calculated. {@code 0 <= alpha < 1.0} 239 * @param nb integer number of functions to be calculated, {@code nb > 0}. The first 240 * function calculated is of order alpha, and the last is of order 241 * {@code nb - 1 + alpha}. 242 * @return BesselJResult a vector of the functions 243 * \(J_{alpha}(x)\) through \(J_{nb-1+alpha}(x)\), or the corresponding exponentially 244 * scaled functions and an integer output variable indicating possible errors 245 */ 246 public static BesselJResult rjBesl(double x, double alpha, int nb) { 247 final double[] b = new double[nb]; 248 249 int ncalc = 0; 250 double alpem = 0; 251 double alp2em = 0; 252 253 // --------------------------------------------------------------------- 254 // Check for out of range arguments. 255 // --------------------------------------------------------------------- 256 final int magx = (int) x; 257 if (nb > 0 && x >= X_MIN && x <= X_MAX && alpha >= 0 && 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 = JdkMath.pow(halfx, alpha) / 286 (alpha * Gamma.value(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 = JdkMath.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 = JdkMath.sin(z); 344 double vcos = JdkMath.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(Arrays.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 = JdkMath.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 = JdkMath.max(test, JdkMath.sqrt(plast * ENSIG) * 482 JdkMath.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 (JdkMath.abs(alpha) > 1e-16) { 623 sum *= Gamma.value(alpha) * JdkMath.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 (JdkMath.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 = JdkMath.min(nb, 0) - 1; 645 } 646 return new BesselJResult(Arrays.copyOf(b, b.length), ncalc); 647 } 648}