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