1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18 package org.apache.commons.math4.legacy.special;
19
20 import java.util.Arrays;
21 import org.apache.commons.numbers.gamma.Gamma;
22 import org.apache.commons.math4.legacy.analysis.UnivariateFunction;
23 import org.apache.commons.math4.legacy.exception.ConvergenceException;
24 import org.apache.commons.math4.legacy.exception.MathIllegalArgumentException;
25 import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
26 import org.apache.commons.math4.core.jdkmath.JdkMath;
27
28 /**
29 * This class provides computation methods related to Bessel
30 * functions of the first kind. Detailed descriptions of these functions are
31 * available in <a
32 * href="http://en.wikipedia.org/wiki/Bessel_function">Wikipedia</a>, <a
33 * href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun">Abrabowitz and
34 * Stegun</a> (Ch. 9-11), and <a href="http://dlmf.nist.gov/">DLMF</a> (Ch. 10).
35 * <p>
36 * This implementation is based on the rjbesl Fortran routine at
37 * <a href="http://www.netlib.org/specfun/rjbesl">Netlib</a>.</p>
38 * <p>
39 * From the Fortran code: </p>
40 * <p>
41 * This program is based on a program written by David J. Sookne (2) that
42 * computes values of the Bessel functions J or I of real argument and integer
43 * order. Modifications include the restriction of the computation to the J
44 * Bessel function of non-negative real argument, the extension of the
45 * computation to arbitrary positive order, and the elimination of most
46 * underflow.</p>
47 * <p>
48 * References:</p>
49 * <ul>
50 * <li>"A Note on Backward Recurrence Algorithms," Olver, F. W. J., and Sookne,
51 * D. J., Math. Comp. 26, 1972, pp 941-947.</li>
52 * <li>"Bessel Functions of Real Argument and Integer Order," Sookne, D. J., NBS
53 * Jour. of Res. B. 77B, 1973, pp 125-132.</li>
54 * </ul>
55 * @since 3.4
56 */
57 public class BesselJ
58 implements UnivariateFunction {
59
60 // ---------------------------------------------------------------------
61 // Mathematical constants
62 // ---------------------------------------------------------------------
63
64 /** -2 / pi. */
65 private static final double PI2 = 0.636619772367581343075535;
66
67 /** first few significant digits of 2pi. */
68 private static final double TOWPI1 = 6.28125;
69
70 /** 2pi - TWOPI1 to working precision. */
71 private static final double TWOPI2 = 1.935307179586476925286767e-3;
72
73 /** TOWPI1 + TWOPI2. */
74 private static final double TWOPI = TOWPI1 + TWOPI2;
75
76 // ---------------------------------------------------------------------
77 // Machine-dependent parameters
78 // ---------------------------------------------------------------------
79
80 /**
81 * 10.0^K, where K is the largest integer such that ENTEN is
82 * machine-representable in working precision.
83 */
84 private static final double ENTEN = 1.0e308;
85
86 /**
87 * Decimal significance desired. Should be set to (INT(log_{10}(2) * (it)+1)).
88 * Setting NSIG lower will result in decreased accuracy while setting
89 * NSIG higher will increase CPU time without increasing accuracy.
90 * The truncation error is limited to a relative error of
91 * T=.5(10^(-NSIG)).
92 */
93 private static final double ENSIG = 1.0e16;
94
95 /** 10.0 ** (-K) for the smallest integer K such that K >= NSIG/4. */
96 private static final double RTNSIG = 1.0e-4;
97
98 /** Smallest ABS(X) such that X/4 does not underflow. */
99 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 }