path: root/src/main/java/org/apache/commons/math3/special/BesselJ.java

/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package org.apache.commons.math3.special;

import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.exception.ConvergenceException;
import org.apache.commons.math3.exception.MathIllegalArgumentException;
import org.apache.commons.math3.exception.util.LocalizedFormats;
import org.apache.commons.math3.util.FastMath;
import org.apache.commons.math3.util.MathArrays;

/**
 * This class provides computation methods related to Bessel functions of the first kind. Detailed
 * descriptions of these functions are available in <a
 * href="http://en.wikipedia.org/wiki/Bessel_function">Wikipedia</a>, <a
 * href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun">Abrabowitz and Stegun</a> (Ch. 9-11),
 * and <a href="http://dlmf.nist.gov/">DLMF</a> (Ch. 10).
 *
 * <p>This implementation is based on the rjbesl Fortran routine at <a
 * href="http://www.netlib.org/specfun/rjbesl">Netlib</a>.
 *
 * <p>From the Fortran code:
 *
 * <p>This program is based on a program written by David J. Sookne (2) that computes values of the
 * Bessel functions J or I of real argument and integer order. Modifications include the restriction
 * of the computation to the J Bessel function of non-negative real argument, the extension of the
 * computation to arbitrary positive order, and the elimination of most underflow.
 *
 * <p>References:
 *
 * <ul>
 *   <li>"A Note on Backward Recurrence Algorithms," Olver, F. W. J., and Sookne, D. J., Math. Comp.
 *       26, 1972, pp 941-947.
 *   <li>"Bessel Functions of Real Argument and Integer Order," Sookne, D. J., NBS Jour. of Res. B.
 *       77B, 1973, pp 125-132.
 * </ul>
 *
 * @since 3.4
 */
public class BesselJ implements UnivariateFunction {

    // ---------------------------------------------------------------------
    // Mathematical constants
    // ---------------------------------------------------------------------

    /** -2 / pi */
    private static final double PI2 = 0.636619772367581343075535;

    /** first few significant digits of 2pi */
    private static final double TOWPI1 = 6.28125;

    /** 2pi - TWOPI1 to working precision */
    private static final double TWOPI2 = 1.935307179586476925286767e-3;

    /** TOWPI1 + TWOPI2 */
    private static final double TWOPI = TOWPI1 + TWOPI2;

    // ---------------------------------------------------------------------
    // Machine-dependent parameters
    // ---------------------------------------------------------------------

    /**
     * 10.0^K, where K is the largest integer such that ENTEN is machine-representable in working
     * precision
     */
    private static final double ENTEN = 1.0e308;

    /**
     * Decimal significance desired. Should be set to (INT(log_{10}(2) * (it)+1)). Setting NSIG
     * lower will result in decreased accuracy while setting NSIG higher will increase CPU time
     * without increasing accuracy. The truncation error is limited to a relative error of
     * T=.5(10^(-NSIG)).
     */
    private static final double ENSIG = 1.0e16;

    /** 10.0 ** (-K) for the smallest integer K such that K >= NSIG/4 */
    private static final double RTNSIG = 1.0e-4;

    /** Smallest ABS(X) such that X/4 does not underflow */
    private static final double ENMTEN = 8.90e-308;

    /** Minimum acceptable value for x */
    private static final double X_MIN = 0.0;

    /**
     * Upper limit on the magnitude of x. If abs(x) = n, then at least n iterations of the backward
     * recursion will be executed. The value of 10.0 ** 4 is used on every machine.
     */
    private static final double X_MAX = 1.0e4;

    /** First 25 factorials as doubles */
    private static final double[] FACT = {
        1.0,
        1.0,
        2.0,
        6.0,
        24.0,
        120.0,
        720.0,
        5040.0,
        40320.0,
        362880.0,
        3628800.0,
        39916800.0,
        479001600.0,
        6227020800.0,
        87178291200.0,
        1.307674368e12,
        2.0922789888e13,
        3.55687428096e14,
        6.402373705728e15,
        1.21645100408832e17,
        2.43290200817664e18,
        5.109094217170944e19,
        1.12400072777760768e21,
        2.585201673888497664e22,
        6.2044840173323943936e23
    };

    /** Order of the function computed when {@link #value(double)} is used */
    private final double order;

    /**
     * Create a new BesselJ with the given order.
     *
     * @param order order of the function computed when using {@link #value(double)}.
     */
    public BesselJ(double order) {
        this.order = order;
    }

    /**
     * Returns the value of the constructed Bessel function of the first kind, for the passed
     * argument.
     *
     * @param x Argument
     * @return Value of the Bessel function at x
     * @throws MathIllegalArgumentException if {@code x} is too large relative to {@code order}
     * @throws ConvergenceException if the algorithm fails to converge
     */
    public double value(double x) throws MathIllegalArgumentException, ConvergenceException {
        return BesselJ.value(order, x);
    }

    /**
     * Returns the first Bessel function, \(J_{order}(x)\).
     *
     * @param order Order of the Bessel function
     * @param x Argument
     * @return Value of the Bessel function of the first kind, \(J_{order}(x)\)
     * @throws MathIllegalArgumentException if {@code x} is too large relative to {@code order}
     * @throws ConvergenceException if the algorithm fails to converge
     */
    public static double value(double order, double x)
            throws MathIllegalArgumentException, ConvergenceException {
        final int n = (int) order;
        final double alpha = order - n;
        final int nb = n + 1;
        final BesselJResult res = rjBesl(x, alpha, nb);

        if (res.nVals >= nb) {
            return res.vals[n];
        } else if (res.nVals < 0) {
            throw new MathIllegalArgumentException(
                    LocalizedFormats.BESSEL_FUNCTION_BAD_ARGUMENT, order, x);
        } else if (FastMath.abs(res.vals[res.nVals - 1]) < 1e-100) {
            return res.vals[n]; // underflow; return value (will be zero)
        }
        throw new ConvergenceException(
                LocalizedFormats.BESSEL_FUNCTION_FAILED_CONVERGENCE, order, x);
    }

    /**
     * Encapsulates the results returned by {@link BesselJ#rjBesl(double, double, int)}.
     *
     * <p>{@link #getVals()} returns the computed function values. {@link #getnVals()} is the number
     * of values among those returned by {@link #getnVals()} that can be considered accurate.
     *
     * <p>
     *
     * <ul>
     *   <li>nVals < 0: An argument is out of range. For example, nb <= 0, alpha < 0 or > 1, or x is
     *       too large. In this case, b(0) is set to zero, the remainder of the b-vector is not
     *       calculated, and nVals is set to MIN(nb,0) - 1 so that nVals != nb.
     *   <li>nb > nVals > 0: Not all requested function values could be calculated accurately. This
     *       usually occurs because nb is much larger than abs(x). In this case, b(n) is calculated
     *       to the desired accuracy for n < nVals, but precision is lost for nVals < n <= nb. If
     *       b(n) does not vanish for n > nVals (because it is too small to be represented), and
     *       b(n)/b(nVals) = \(10^{-k}\), then only the first NSIG-k significant figures of b(n) can
     *       be trusted.
     * </ul>
     */
    public static class BesselJResult {

        /** Bessel function values */
        private final double[] vals;

        /** Valid value count */
        private final int nVals;

        /**
         * Create a new BesselJResult with the given values and valid value count.
         *
         * @param b values
         * @param n count of valid values
         */
        public BesselJResult(double[] b, int n) {
            vals = MathArrays.copyOf(b, b.length);
            nVals = n;
        }

        /**
         * @return the computed function values
         */
        public double[] getVals() {
            return MathArrays.copyOf(vals, vals.length);
        }

        /**
         * @return the number of valid function values (normally the same as the length of the array
         *     returned by {@link #getnVals()})
         */
        public int getnVals() {
            return nVals;
        }
    }

    /**
     * Calculates Bessel functions \(J_{n+alpha}(x)\) for non-negative argument x, and non-negative
     * order n + alpha.
     *
     * <p>Before using the output vector, the user should check that nVals = nb, i.e., all orders
     * have been calculated to the desired accuracy. See BesselResult class javadoc for details on
     * return values.
     *
     * @param x non-negative real argument for which J's are to be calculated
     * @param alpha fractional part of order for which J's or exponentially scaled J's (\(J\cdot
     *     e^{x}\)) are to be calculated. 0 <= alpha < 1.0.
     * @param nb integer number of functions to be calculated, nb > 0. The first function calculated
     *     is of order alpha, and the last is of order nb - 1 + alpha.
     * @return BesselJResult a vector of the functions \(J_{alpha}(x)\) through
     *     \(J_{nb-1+alpha}(x)\), or the corresponding exponentially scaled functions and an integer
     *     output variable indicating possible errors
     */
    public static BesselJResult rjBesl(double x, double alpha, int nb) {
        final double[] b = new double[nb];

        int ncalc = 0;
        double alpem = 0;
        double alp2em = 0;

        // ---------------------------------------------------------------------
        // Check for out of range arguments.
        // ---------------------------------------------------------------------
        final int magx = (int) x;
        if ((nb > 0) && (x >= X_MIN) && (x <= X_MAX) && (alpha >= 0) && (alpha < 1)) {
            // ---------------------------------------------------------------------
            // Initialize result array to zero.
            // ---------------------------------------------------------------------
            ncalc = nb;
            for (int i = 0; i < nb; ++i) {
                b[i] = 0;
            }

            // ---------------------------------------------------------------------
            // Branch to use 2-term ascending series for small X and asymptotic
            // form for large X when NB is not too large.
            // ---------------------------------------------------------------------
            double tempa;
            double tempb;
            double tempc;
            double tover;
            if (x < RTNSIG) {
                // ---------------------------------------------------------------------
                // Two-term ascending series for small X.
                // ---------------------------------------------------------------------
                tempa = 1;
                alpem = 1 + alpha;
                double halfx = 0;
                if (x > ENMTEN) {
                    halfx = 0.5 * x;
                }
                if (alpha != 0) {
                    tempa = FastMath.pow(halfx, alpha) / (alpha * Gamma.gamma(alpha));
                }
                tempb = 0;
                if (x + 1 > 1) {
                    tempb = -halfx * halfx;
                }
                b[0] = tempa + (tempa * tempb / alpem);
                if ((x != 0) && (b[0] == 0)) {
                    ncalc = 0;
                }
                if (nb != 1) {
                    if (x <= 0) {
                        for (int n = 1; n < nb; ++n) {
                            b[n] = 0;
                        }
                    } else {
                        // ---------------------------------------------------------------------
                        // Calculate higher order functions.
                        // ---------------------------------------------------------------------
                        tempc = halfx;
                        tover = tempb != 0 ? ENMTEN / tempb : 2 * ENMTEN / x;
                        for (int n = 1; n < nb; ++n) {
                            tempa /= alpem;
                            alpem += 1;
                            tempa *= tempc;
                            if (tempa <= tover * alpem) {
                                tempa = 0;
                            }
                            b[n] = tempa + (tempa * tempb / alpem);
                            if ((b[n] == 0) && (ncalc > n)) {
                                ncalc = n;
                            }
                        }
                    }
                }
            } else if ((x > 25.0) && (nb <= magx + 1)) {
                // ---------------------------------------------------------------------
                // Asymptotic series for X > 25
                // ---------------------------------------------------------------------
                final double xc = FastMath.sqrt(PI2 / x);
                final double mul = 0.125 / x;
                final double xin = mul * mul;
                int m = 0;
                if (x >= 130.0) {
                    m = 4;
                } else if (x >= 35.0) {
                    m = 8;
                } else {
                    m = 11;
                }

                final double xm = 4.0 * m;
                // ---------------------------------------------------------------------
                // Argument reduction for SIN and COS routines.
                // ---------------------------------------------------------------------
                double t = (double) ((int) ((x / TWOPI) + 0.5));
                final double z = x - t * TOWPI1 - t * TWOPI2 - (alpha + 0.5) / PI2;
                double vsin = FastMath.sin(z);
                double vcos = FastMath.cos(z);
                double gnu = 2 * alpha;
                double capq;
                double capp;
                double s;
                double t1;
                double xk;
                for (int i = 1; i <= 2; i++) {
                    s = (xm - 1 - gnu) * (xm - 1 + gnu) * xin * 0.5;
                    t = (gnu - (xm - 3.0)) * (gnu + (xm - 3.0));
                    capp = (s * t) / FACT[2 * m];
                    t1 = (gnu - (xm + 1)) * (gnu + (xm + 1));
                    capq = (s * t1) / FACT[2 * m + 1];
                    xk = xm;
                    int k = 2 * m;
                    t1 = t;

                    for (int j = 2; j <= m; j++) {
                        xk -= 4.0;
                        s = (xk - 1 - gnu) * (xk - 1 + gnu);
                        t = (gnu - (xk - 3.0)) * (gnu + (xk - 3.0));
                        capp = (capp + 1 / FACT[k - 2]) * s * t * xin;
                        capq = (capq + 1 / FACT[k - 1]) * s * t1 * xin;
                        k -= 2;
                        t1 = t;
                    }

                    capp += 1;
                    capq = (capq + 1) * ((gnu * gnu) - 1) * (0.125 / x);
                    b[i - 1] = xc * (capp * vcos - capq * vsin);
                    if (nb == 1) {
                        return new BesselJResult(MathArrays.copyOf(b, b.length), ncalc);
                    }
                    t = vsin;
                    vsin = -vcos;
                    vcos = t;
                    gnu += 2.0;
                }

                // ---------------------------------------------------------------------
                // If NB > 2, compute J(X,ORDER+I) I = 2, NB-1
                // ---------------------------------------------------------------------
                if (nb > 2) {
                    gnu = 2 * alpha + 2.0;
                    for (int j = 2; j < nb; ++j) {
                        b[j] = gnu * b[j - 1] / x - b[j - 2];
                        gnu += 2.0;
                    }
                }
            } else {
                // ---------------------------------------------------------------------
                // Use recurrence to generate results. First initialize the
                // calculation of P*S.
                // ---------------------------------------------------------------------
                final int nbmx = nb - magx;
                int n = magx + 1;
                int nstart = 0;
                int nend = 0;
                double en = 2 * (n + alpha);
                double plast = 1;
                double p = en / x;
                double pold;
                // ---------------------------------------------------------------------
                // Calculate general significance test.
                // ---------------------------------------------------------------------
                double test = 2 * ENSIG;
                boolean readyToInitialize = false;
                if (nbmx >= 3) {
                    // ---------------------------------------------------------------------
                    // Calculate P*S until N = NB-1. Check for possible
                    // overflow.
                    // ---------------------------------------------------------------------
                    tover = ENTEN / ENSIG;
                    nstart = magx + 2;
                    nend = nb - 1;
                    en = 2 * (nstart - 1 + alpha);
                    double psave;
                    double psavel;
                    for (int k = nstart; k <= nend; k++) {
                        n = k;
                        en += 2.0;
                        pold = plast;
                        plast = p;
                        p = (en * plast / x) - pold;
                        if (p > tover) {
                            // ---------------------------------------------------------------------
                            // To avoid overflow, divide P*S by TOVER. Calculate
                            // P*S until
                            // ABS(P) > 1.
                            // ---------------------------------------------------------------------
                            tover = ENTEN;
                            p /= tover;
                            plast /= tover;
                            psave = p;
                            psavel = plast;
                            nstart = n + 1;
                            do {
                                n += 1;
                                en += 2.0;
                                pold = plast;
                                plast = p;
                                p = (en * plast / x) - pold;
                            } while (p <= 1);
                            tempb = en / x;
                            // ---------------------------------------------------------------------
                            // Calculate backward test and find NCALC, the
                            // highest N such that
                            // the test is passed.
                            // ---------------------------------------------------------------------
                            test = pold * plast * (0.5 - 0.5 / (tempb * tempb));
                            test /= ENSIG;
                            p = plast * tover;
                            n -= 1;
                            en -= 2.0;
                            nend = FastMath.min(nb, n);
                            for (int l = nstart; l <= nend; l++) {
                                pold = psavel;
                                psavel = psave;
                                psave = (en * psavel / x) - pold;
                                if (psave * psavel > test) {
                                    ncalc = l - 1;
                                    readyToInitialize = true;
                                    break;
                                }
                            }
                            ncalc = nend;
                            readyToInitialize = true;
                            break;
                        }
                    }
                    if (!readyToInitialize) {
                        n = nend;
                        en = 2 * (n + alpha);
                        // ---------------------------------------------------------------------
                        // Calculate special significance test for NBMX > 2.
                        // ---------------------------------------------------------------------
                        test =
                                FastMath.max(
                                        test, FastMath.sqrt(plast * ENSIG) * FastMath.sqrt(2 * p));
                    }
                }
                // ---------------------------------------------------------------------
                // Calculate P*S until significance test passes.
                // ---------------------------------------------------------------------
                if (!readyToInitialize) {
                    do {
                        n += 1;
                        en += 2.0;
                        pold = plast;
                        plast = p;
                        p = (en * plast / x) - pold;
                    } while (p < test);
                }
                // ---------------------------------------------------------------------
                // Initialize the backward recursion and the normalization sum.
                // ---------------------------------------------------------------------
                n += 1;
                en += 2.0;
                tempb = 0;
                tempa = 1 / p;
                int m = (2 * n) - 4 * (n / 2);
                double sum = 0;
                double em = (double) (n / 2);
                alpem = em - 1 + alpha;
                alp2em = 2 * em + alpha;
                if (m != 0) {
                    sum = tempa * alpem * alp2em / em;
                }
                nend = n - nb;

                boolean readyToNormalize = false;
                boolean calculatedB0 = false;

                // ---------------------------------------------------------------------
                // Recur backward via difference equation, calculating (but not
                // storing) B(N), until N = NB.
                // ---------------------------------------------------------------------
                for (int l = 1; l <= nend; l++) {
                    n -= 1;
                    en -= 2.0;
                    tempc = tempb;
                    tempb = tempa;
                    tempa = (en * tempb / x) - tempc;
                    m = 2 - m;
                    if (m != 0) {
                        em -= 1;
                        alp2em = 2 * em + alpha;
                        if (n == 1) {
                            break;
                        }
                        alpem = em - 1 + alpha;
                        if (alpem == 0) {
                            alpem = 1;
                        }
                        sum = (sum + tempa * alp2em) * alpem / em;
                    }
                }

                // ---------------------------------------------------------------------
                // Store B(NB).
                // ---------------------------------------------------------------------
                b[n - 1] = tempa;
                if (nend >= 0) {
                    if (nb <= 1) {
                        alp2em = alpha;
                        if (alpha + 1 == 1) {
                            alp2em = 1;
                        }
                        sum += b[0] * alp2em;
                        readyToNormalize = true;
                    } else {
                        // ---------------------------------------------------------------------
                        // Calculate and store B(NB-1).
                        // ---------------------------------------------------------------------
                        n -= 1;
                        en -= 2.0;
                        b[n - 1] = (en * tempa / x) - tempb;
                        if (n == 1) {
                            calculatedB0 = true;
                        } else {
                            m = 2 - m;
                            if (m != 0) {
                                em -= 1;
                                alp2em = 2 * em + alpha;
                                alpem = em - 1 + alpha;
                                if (alpem == 0) {
                                    alpem = 1;
                                }

                                sum = (sum + (b[n - 1] * alp2em)) * alpem / em;
                            }
                        }
                    }
                }
                if (!readyToNormalize && !calculatedB0) {
                    nend = n - 2;
                    if (nend != 0) {
                        // ---------------------------------------------------------------------
                        // Calculate via difference equation and store B(N),
                        // until N = 2.
                        // ---------------------------------------------------------------------

                        for (int l = 1; l <= nend; l++) {
                            n -= 1;
                            en -= 2.0;
                            b[n - 1] = (en * b[n] / x) - b[n + 1];
                            m = 2 - m;
                            if (m != 0) {
                                em -= 1;
                                alp2em = 2 * em + alpha;
                                alpem = em - 1 + alpha;
                                if (alpem == 0) {
                                    alpem = 1;
                                }

                                sum = (sum + b[n - 1] * alp2em) * alpem / em;
                            }
                        }
                    }
                }
                // ---------------------------------------------------------------------
                // Calculate b[0]
                // ---------------------------------------------------------------------
                if (!readyToNormalize) {
                    if (!calculatedB0) {
                        b[0] = 2.0 * (alpha + 1) * b[1] / x - b[2];
                    }
                    em -= 1;
                    alp2em = 2 * em + alpha;
                    if (alp2em == 0) {
                        alp2em = 1;
                    }
                    sum += b[0] * alp2em;
                }
                // ---------------------------------------------------------------------
                // Normalize. Divide all B(N) by sum.
                // ---------------------------------------------------------------------

                if (FastMath.abs(alpha) > 1e-16) {
                    sum *= Gamma.gamma(alpha) * FastMath.pow(x * 0.5, -alpha);
                }
                tempa = ENMTEN;
                if (sum > 1) {
                    tempa *= sum;
                }

                for (n = 0; n < nb; n++) {
                    if (FastMath.abs(b[n]) < tempa) {
                        b[n] = 0;
                    }
                    b[n] /= sum;
                }
            }
            // ---------------------------------------------------------------------
            // Error return -- X, NB, or ALPHA is out of range.
            // ---------------------------------------------------------------------
        } else {
            if (b.length > 0) {
                b[0] = 0;
            }
            ncalc = FastMath.min(nb, 0) - 1;
        }
        return new BesselJResult(MathArrays.copyOf(b, b.length), ncalc);
    }
}