View Javadoc
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 &lt; 0: An argument is out of range. For example, nb &lt;= 0, alpha
187      * &lt; 0 or &gt; 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 &gt; nVals &gt; 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 &lt; nVals, but
193      * precision is lost for nVals &lt; n &lt;= nb. If b(n) does not vanish for n &gt;
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 &lt;= alpha &lt; 1.0.
244      * @param nb integer number of functions to be calculated, nb &gt; 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 }