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.special;
019
020import java.util.Arrays;
021import org.apache.commons.numbers.gamma.Gamma;
022import org.apache.commons.math4.analysis.UnivariateFunction;
023import org.apache.commons.math4.exception.ConvergenceException;
024import org.apache.commons.math4.exception.MathIllegalArgumentException;
025import org.apache.commons.math4.exception.util.LocalizedFormats;
026import org.apache.commons.math4.util.FastMath;
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 (FastMath.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) &&
258            (alpha < 1)) {
259            // ---------------------------------------------------------------------
260            // Initialize result array to zero.
261            // ---------------------------------------------------------------------
262            ncalc = nb;
263            for (int i = 0; i < nb; ++i) {
264                b[i] = 0;
265            }
266
267            // ---------------------------------------------------------------------
268            // Branch to use 2-term ascending series for small X and asymptotic
269            // form for large X when NB is not too large.
270            // ---------------------------------------------------------------------
271            double tempa;
272            double tempb;
273            double tempc;
274            double tover;
275            if (x < RTNSIG) {
276                // ---------------------------------------------------------------------
277                // Two-term ascending series for small X.
278                // ---------------------------------------------------------------------
279                tempa = 1;
280                alpem = 1 + alpha;
281                double halfx = 0;
282                if (x > ENMTEN) {
283                    halfx = 0.5 * x;
284                }
285                if (alpha != 0) {
286                    tempa = FastMath.pow(halfx, alpha) /
287                            (alpha * Gamma.value(alpha));
288                }
289                tempb = 0;
290                if (x + 1 > 1) {
291                    tempb = -halfx * halfx;
292                }
293                b[0] = tempa + (tempa * tempb / alpem);
294                if ((x != 0) && (b[0] == 0)) {
295                    ncalc = 0;
296                }
297                if (nb != 1) {
298                    if (x <= 0) {
299                        for (int n = 1; n < nb; ++n) {
300                            b[n] = 0;
301                        }
302                    } else {
303                        // ---------------------------------------------------------------------
304                        // Calculate higher order functions.
305                        // ---------------------------------------------------------------------
306                        tempc = halfx;
307                        tover = tempb != 0 ? ENMTEN / tempb :  2 * ENMTEN / x;
308                        for (int n = 1; n < nb; ++n) {
309                            tempa /= alpem;
310                            alpem += 1;
311                            tempa *= tempc;
312                            if (tempa <= tover * alpem) {
313                                tempa = 0;
314                            }
315                            b[n] = tempa + (tempa * tempb / alpem);
316                            if ((b[n] == 0) && (ncalc > n)) {
317                                ncalc = n;
318                            }
319                        }
320                    }
321                }
322            } else if ((x > 25.0) && (nb <= magx + 1)) {
323                // ---------------------------------------------------------------------
324                // Asymptotic series for X > 25
325                // ---------------------------------------------------------------------
326                final double xc = FastMath.sqrt(PI2 / x);
327                final double mul = 0.125 / x;
328                final double xin = mul * mul;
329                int m = 0;
330                if (x >= 130.0) {
331                    m = 4;
332                } else if (x >= 35.0) {
333                    m = 8;
334                } else {
335                    m = 11;
336                }
337
338                final double xm = 4.0 * m;
339                // ---------------------------------------------------------------------
340                // Argument reduction for SIN and COS routines.
341                // ---------------------------------------------------------------------
342                double t = (double) ((int) ((x / TWOPI) + 0.5));
343                final double z = x - t * TOWPI1 - t * TWOPI2 - (alpha + 0.5) / PI2;
344                double vsin = FastMath.sin(z);
345                double vcos = FastMath.cos(z);
346                double gnu = 2 * alpha;
347                double capq;
348                double capp;
349                double s;
350                double t1;
351                double xk;
352                for (int i = 1; i <= 2; i++) {
353                    s = (xm - 1 - gnu) * (xm - 1 + gnu) * xin * 0.5;
354                    t = (gnu - (xm - 3.0)) * (gnu + (xm - 3.0));
355                    capp = (s * t) / FACT[2 * m];
356                    t1 = (gnu - (xm + 1)) * (gnu + (xm + 1));
357                    capq = (s * t1) / FACT[2 * m + 1];
358                    xk = xm;
359                    int k = 2 * m;
360                    t1 = t;
361
362                    for (int j = 2; j <= m; j++) {
363                        xk -= 4.0;
364                        s = (xk - 1 - gnu) * (xk - 1 + gnu);
365                        t = (gnu - (xk - 3.0)) * (gnu + (xk - 3.0));
366                        capp = (capp + 1 / FACT[k - 2]) * s * t * xin;
367                        capq = (capq + 1 / FACT[k - 1]) * s * t1 * xin;
368                        k -= 2;
369                        t1 = t;
370                    }
371
372                    capp += 1;
373                    capq = (capq + 1) * ((gnu * gnu) - 1) * (0.125 / x);
374                    b[i - 1] = xc * (capp * vcos - capq * vsin);
375                    if (nb == 1) {
376                        return new BesselJResult(Arrays.copyOf(b, b.length),
377                                                 ncalc);
378                    }
379                    t = vsin;
380                    vsin = -vcos;
381                    vcos = t;
382                    gnu += 2.0;
383                }
384
385                // ---------------------------------------------------------------------
386                // If NB > 2, compute J(X,ORDER+I) I = 2, NB-1
387                // ---------------------------------------------------------------------
388                if (nb > 2) {
389                    gnu = 2 * alpha + 2.0;
390                    for (int j = 2; j < nb; ++j) {
391                        b[j] = gnu * b[j - 1] / x - b[j - 2];
392                        gnu += 2.0;
393                    }
394                }
395            } else {
396                // ---------------------------------------------------------------------
397                // Use recurrence to generate results. First initialize the
398                // calculation of P*S.
399                // ---------------------------------------------------------------------
400                final int nbmx = nb - magx;
401                int n = magx + 1;
402                int nstart = 0;
403                int nend = 0;
404                double en = 2 * (n + alpha);
405                double plast = 1;
406                double p = en / x;
407                double pold;
408                // ---------------------------------------------------------------------
409                // Calculate general significance test.
410                // ---------------------------------------------------------------------
411                double test = 2 * ENSIG;
412                boolean readyToInitialize = false;
413                if (nbmx >= 3) {
414                    // ---------------------------------------------------------------------
415                    // Calculate P*S until N = NB-1. Check for possible
416                    // overflow.
417                    // ---------------------------------------------------------------------
418                    tover = ENTEN / ENSIG;
419                    nstart = magx + 2;
420                    nend = nb - 1;
421                    en = 2 * (nstart - 1 + alpha);
422                    double psave;
423                    double psavel;
424                    for (int k = nstart; k <= nend; k++) {
425                        n = k;
426                        en += 2.0;
427                        pold = plast;
428                        plast = p;
429                        p = (en * plast / x) - pold;
430                        if (p > tover) {
431                            // ---------------------------------------------------------------------
432                            // To avoid overflow, divide P*S by TOVER. Calculate
433                            // P*S until
434                            // ABS(P) > 1.
435                            // ---------------------------------------------------------------------
436                            tover = ENTEN;
437                            p /= tover;
438                            plast /= tover;
439                            psave = p;
440                            psavel = plast;
441                            nstart = n + 1;
442                            do {
443                                n += 1;
444                                en += 2.0;
445                                pold = plast;
446                                plast = p;
447                                p = (en * plast / x) - pold;
448                            } while (p <= 1);
449                            tempb = en / x;
450                            // ---------------------------------------------------------------------
451                            // Calculate backward test and find NCALC, the
452                            // highest N such that
453                            // the test is passed.
454                            // ---------------------------------------------------------------------
455                            test = pold * plast * (0.5 - 0.5 / (tempb * tempb));
456                            test /= ENSIG;
457                            p = plast * tover;
458                            n -= 1;
459                            en -= 2.0;
460                            nend = FastMath.min(nb, n);
461                            for (int l = nstart; l <= nend; l++) {
462                                pold = psavel;
463                                psavel = psave;
464                                psave = (en * psavel / x) - pold;
465                                if (psave * psavel > test) {
466                                    ncalc = l - 1;
467                                    readyToInitialize = true;
468                                    break;
469                                }
470                            }
471                            ncalc = nend;
472                            readyToInitialize = true;
473                            break;
474                        }
475                    }
476                    if (!readyToInitialize) {
477                        n = nend;
478                        en = 2 * (n + alpha);
479                        // ---------------------------------------------------------------------
480                        // Calculate special significance test for NBMX > 2.
481                        // ---------------------------------------------------------------------
482                        test = FastMath.max(test, FastMath.sqrt(plast * ENSIG) *
483                                                  FastMath.sqrt(2 * p));
484                    }
485                }
486                // ---------------------------------------------------------------------
487                // Calculate P*S until significance test passes.
488                // ---------------------------------------------------------------------
489                if (!readyToInitialize) {
490                    do {
491                        n += 1;
492                        en += 2.0;
493                        pold = plast;
494                        plast = p;
495                        p = (en * plast / x) - pold;
496                    } while (p < test);
497                }
498                // ---------------------------------------------------------------------
499                // Initialize the backward recursion and the normalization sum.
500                // ---------------------------------------------------------------------
501                n += 1;
502                en += 2.0;
503                tempb = 0;
504                tempa = 1 / p;
505                int m = (2 * n) - 4 * (n / 2);
506                double sum = 0;
507                double em = (double) (n / 2);
508                alpem = em - 1 + alpha;
509                alp2em = 2 * em + alpha;
510                if (m != 0) {
511                    sum = tempa * alpem * alp2em / em;
512                }
513                nend = n - nb;
514
515                boolean readyToNormalize = false;
516                boolean calculatedB0 = false;
517
518                // ---------------------------------------------------------------------
519                // Recur backward via difference equation, calculating (but not
520                // storing) B(N), until N = NB.
521                // ---------------------------------------------------------------------
522                for (int l = 1; l <= nend; l++) {
523                    n -= 1;
524                    en -= 2.0;
525                    tempc = tempb;
526                    tempb = tempa;
527                    tempa = (en * tempb / x) - tempc;
528                    m = 2 - m;
529                    if (m != 0) {
530                        em -= 1;
531                        alp2em = 2 * em + alpha;
532                        if (n == 1) {
533                            break;
534                        }
535                        alpem = em - 1 + alpha;
536                        if (alpem == 0) {
537                            alpem = 1;
538                        }
539                        sum = (sum + tempa * alp2em) * alpem / em;
540                    }
541                }
542
543                // ---------------------------------------------------------------------
544                // Store B(NB).
545                // ---------------------------------------------------------------------
546                b[n - 1] = tempa;
547                if (nend >= 0) {
548                    if (nb <= 1) {
549                        alp2em = alpha;
550                        if (alpha + 1 == 1) {
551                            alp2em = 1;
552                        }
553                        sum += b[0] * alp2em;
554                        readyToNormalize = true;
555                    } else {
556                        // ---------------------------------------------------------------------
557                        // Calculate and store B(NB-1).
558                        // ---------------------------------------------------------------------
559                        n -= 1;
560                        en -= 2.0;
561                        b[n - 1] = (en * tempa / x) - tempb;
562                        if (n == 1) {
563                            calculatedB0 = true;
564                        } else {
565                            m = 2 - m;
566                            if (m != 0) {
567                                em -= 1;
568                                alp2em = 2 * em + alpha;
569                                alpem = em - 1 + alpha;
570                                if (alpem == 0) {
571                                    alpem = 1;
572                                }
573
574                                sum = (sum + (b[n - 1] * alp2em)) * alpem / em;
575                            }
576                        }
577                    }
578                }
579                if (!readyToNormalize && !calculatedB0) {
580                    nend = n - 2;
581                    if (nend != 0) {
582                        // ---------------------------------------------------------------------
583                        // Calculate via difference equation and store B(N),
584                        // until N = 2.
585                        // ---------------------------------------------------------------------
586
587                        for (int l = 1; l <= nend; l++) {
588                            n -= 1;
589                            en -= 2.0;
590                            b[n - 1] = (en * b[n] / x) - b[n + 1];
591                            m = 2 - m;
592                            if (m != 0) {
593                                em -= 1;
594                                alp2em = 2 * em + alpha;
595                                alpem = em - 1 + alpha;
596                                if (alpem == 0) {
597                                    alpem = 1;
598                                }
599
600                                sum = (sum + b[n - 1] * alp2em) * alpem / em;
601                            }
602                        }
603                    }
604                }
605                // ---------------------------------------------------------------------
606                // Calculate b[0]
607                // ---------------------------------------------------------------------
608                if (!readyToNormalize) {
609                    if (!calculatedB0) {
610                        b[0] = 2.0 * (alpha + 1) * b[1] / x - b[2];
611                    }
612                    em -= 1;
613                    alp2em = 2 * em + alpha;
614                    if (alp2em == 0) {
615                        alp2em = 1;
616                    }
617                    sum += b[0] * alp2em;
618                }
619                // ---------------------------------------------------------------------
620                // Normalize. Divide all B(N) by sum.
621                // ---------------------------------------------------------------------
622
623                if (FastMath.abs(alpha) > 1e-16) {
624                    sum *= Gamma.value(alpha) * FastMath.pow(x * 0.5, -alpha);
625                }
626                tempa = ENMTEN;
627                if (sum > 1) {
628                    tempa *= sum;
629                }
630
631                for (n = 0; n < nb; n++) {
632                    if (FastMath.abs(b[n]) < tempa) {
633                        b[n] = 0;
634                    }
635                    b[n] /= sum;
636                }
637            }
638            // ---------------------------------------------------------------------
639            // Error return -- X, NB, or ALPHA is out of range.
640            // ---------------------------------------------------------------------
641        } else {
642            if (b.length > 0) {
643                b[0] = 0;
644            }
645            ncalc = FastMath.min(nb, 0) - 1;
646        }
647        return new BesselJResult(Arrays.copyOf(b, b.length), ncalc);
648    }
649}