| /* |
| * 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.math4.legacy.special; |
| |
| import java.util.Arrays; |
| import org.apache.commons.numbers.gamma.Gamma; |
| import org.apache.commons.math4.legacy.analysis.UnivariateFunction; |
| import org.apache.commons.math4.legacy.exception.ConvergenceException; |
| import org.apache.commons.math4.legacy.exception.MathIllegalArgumentException; |
| import org.apache.commons.math4.legacy.exception.util.LocalizedFormats; |
| import org.apache.commons.math4.core.jdkmath.JdkMath; |
| |
| /** |
| * 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> |
| * <p> |
| * From the Fortran code: </p> |
| * <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> |
| * <p> |
| * References:</p> |
| * <ul> |
| * <li>"A Note on Backward Recurrence Algorithms," Olver, F. W. J., and Sookne, |
| * D. J., Math. Comp. 26, 1972, pp 941-947.</li> |
| * <li>"Bessel Functions of Real Argument and Integer Order," Sookne, D. J., NBS |
| * Jour. of Res. B. 77B, 1973, pp 125-132.</li> |
| * </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 |
| */ |
| @Override |
| 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 (JdkMath.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>{@code nVals < 0}: An argument is out of range. For example, {@code nb <= 0}, |
| * {@code 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> |
| * <li>{@code 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 {@code n < nVals}, but |
| * precision is lost for {@code nVals < n <= nb}. If b(n) does not vanish for |
| * {@code 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.</li></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 = Arrays.copyOf(b, b.length); |
| nVals = n; |
| } |
| |
| /** |
| * @return the computed function values |
| */ |
| public double[] getVals() { |
| return Arrays.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. |
| * </p> |
| * @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. {@code 0 <= alpha < 1.0} |
| * @param nb integer number of functions to be calculated, {@code nb > 0}. The first |
| * function calculated is of order alpha, and the last is of order |
| * {@code 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 = JdkMath.pow(halfx, alpha) / |
| (alpha * Gamma.value(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 = JdkMath.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 = JdkMath.sin(z); |
| double vcos = JdkMath.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(Arrays.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 = JdkMath.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 = JdkMath.max(test, JdkMath.sqrt(plast * ENSIG) * |
| JdkMath.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 (JdkMath.abs(alpha) > 1e-16) { |
| sum *= Gamma.value(alpha) * JdkMath.pow(x * 0.5, -alpha); |
| } |
| tempa = ENMTEN; |
| if (sum > 1) { |
| tempa *= sum; |
| } |
| |
| for (n = 0; n < nb; n++) { |
| if (JdkMath.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 = JdkMath.min(nb, 0) - 1; |
| } |
| return new BesselJResult(Arrays.copyOf(b, b.length), ncalc); |
| } |
| } |
