diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/special/BesselJ.java')
| -rw-r--r-- | src/main/java/org/apache/commons/math3/special/BesselJ.java | 661 |
1 files changed, 661 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/special/BesselJ.java b/src/main/java/org/apache/commons/math3/special/BesselJ.java new file mode 100644 index 0000000..61aa4ff --- /dev/null +++ b/src/main/java/org/apache/commons/math3/special/BesselJ.java @@ -0,0 +1,661 @@ +/* + * 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); + } +} |