BesselJ.java

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

  18. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */

  23. package org.hipparchus.special;

  25. import org.hipparchus.analysis.UnivariateFunction;
  26. import org.hipparchus.exception.LocalizedCoreFormats;
  27. import org.hipparchus.exception.MathIllegalArgumentException;
  28. import org.hipparchus.exception.MathIllegalStateException;
  29. import org.hipparchus.util.FastMath;
  30. import org.hipparchus.util.SinCos;

  32. /**
  33.  * This class provides computation methods related to Bessel
  34.  * functions of the first kind. Detailed descriptions of these functions are
  35.  * available in <a
  36.  * href="http://en.wikipedia.org/wiki/Bessel_function">Wikipedia</a>, <a
  37.  * href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun">Abramowitz and
  38.  * Stegun</a> (Ch. 9-11), and <a href="http://dlmf.nist.gov/">DLMF</a> (Ch. 10).
  39.  * <p>
  40.  * This implementation is based on the rjbesl Fortran routine at
  41.  * <a href="http://www.netlib.org/specfun/rjbesl">Netlib</a>.
  42.  * </p>
  43.  * <p>
  44.  * From the Fortran code: </p>
  45.  * <p>
  46.  * This program is based on a program written by David J. Sookne (2) that
  47.  * computes values of the Bessel functions J or I of real argument and integer
  48.  * order. Modifications include the restriction of the computation to the J
  49.  * Bessel function of non-negative real argument, the extension of the
  50.  * computation to arbitrary positive order, and the elimination of most
  51.  * underflow.</p>
  52.  * <p>
  53.  * References:</p>
  54.  * <ul>
  55.  * <li>"A Note on Backward Recurrence Algorithms," Olver, F. W. J., and Sookne,
  56.  * D. J., Math. Comp. 26, 1972, pp 941-947.</li>
  57.  * <li>"Bessel Functions of Real Argument and Integer Order," Sookne, D. J., NBS
  58.  * Jour. of Res. B. 77B, 1973, pp 125-132.</li>
  59.  * </ul>
  60.  */
  61. public class BesselJ
  62.     implements UnivariateFunction {

  64.     // ---------------------------------------------------------------------
  65.     // Mathematical constants
  66.     // ---------------------------------------------------------------------

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

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

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

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

  80.     // ---------------------------------------------------------------------
  81.     // Machine-dependent parameters
  82.     // ---------------------------------------------------------------------

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

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

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

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

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

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

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

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

  128.     /**
  129.      * Create a new BesselJ with the given order.
  130.      *
  131.      * @param order order of the function computed when using {@link #value(double)}.
  132.      */
  133.     public BesselJ(double order) {
  134.         this.order = order;
  135.     }

  137.     /**
  138.      * Returns the value of the constructed Bessel function of the first kind,
  139.      * for the passed argument.
  140.      *
  141.      * @param x Argument
  142.      * @return Value of the Bessel function at x
  143.      * @throws MathIllegalArgumentException if {@code x} is too large relative to {@code order}
  144.      * @throws MathIllegalStateException if the algorithm fails to converge
  145.      */
  146.     @Override
  147.     public double value(double x)
  148.         throws MathIllegalArgumentException, MathIllegalStateException {
  149.         return value(order, x);
  150.     }

  152.     /**
  153.      * Returns the first Bessel function, \(J_{order}(x)\).
  154.      *
  155.      * @param order Order of the Bessel function
  156.      * @param x Argument
  157.      * @return Value of the Bessel function of the first kind, \(J_{order}(x)\)
  158.      * @throws MathIllegalArgumentException if {@code x} is too large relative to {@code order}
  159.      * @throws MathIllegalStateException if the algorithm fails to converge
  160.      */
  161.     public static double value(double order, double x)
  162.         throws MathIllegalArgumentException, MathIllegalStateException {
  163.         final int n = (int) order;
  164.         final double alpha = order - n;
  165.         final int nb = n + 1;
  166.         final BesselJResult res = rjBesl(x, alpha, nb);

  168.         if (res.nVals >= nb) {
  169.             return res.vals[n];
  170.         } else if (res.nVals < 0) {
  171.             throw new MathIllegalArgumentException(LocalizedCoreFormats.BESSEL_FUNCTION_BAD_ARGUMENT,order, x);
  172.         } else if (FastMath.abs(res.vals[res.nVals - 1]) < 1e-100) {
  173.             return res.vals[n]; // underflow; return value (will be zero)
  174.         }
  175.         throw new MathIllegalStateException(LocalizedCoreFormats.BESSEL_FUNCTION_FAILED_CONVERGENCE, order, x);
  176.     }

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

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

  203.         /** Valid value count */
  204.         private final int nVals;

  206.         /**
  207.          * Create a new BesselJResult with the given values and valid value count.
  208.          *
  209.          * @param b values
  210.          * @param n count of valid values
  211.          */
  212.         public BesselJResult(double[] b, int n) {
  213.             vals = b.clone();
  214.             nVals = n;
  215.         }

  217.         /** Get computed function values.
  218.          * @return the computed function values
  219.          */
  220.         public double[] getVals() {
  221.             return vals.clone();
  222.         }

  224.         /** Get number of valid function values.
  225.          * @return the number of valid function values (normally the same as the
  226.          *         length of the array returned by {@link #getnVals()})
  227.          */
  228.         public int getnVals() {
  229.             return nVals;
  230.         }
  231.     }

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

  254.         int ncalc;
  255.         double alpem;
  256.         double alp2em;

  258.         // ---------------------------------------------------------------------
  259.         // Check for out of range arguments.
  260.         // ---------------------------------------------------------------------
  261.         final int magx = (int) x;
  262.         if ((nb > 0) && (x >= X_MIN) && (x <= X_MAX) && (alpha >= 0) &&
  263.             (alpha < 1)) {
  264.             // ---------------------------------------------------------------------
  265.             // Initialize result array to zero.
  266.             // ---------------------------------------------------------------------
  267.             ncalc = nb;
  268.             for (int i = 0; i < nb; ++i) {
  269.                 b[i] = 0;
  270.             }

  272.             // ---------------------------------------------------------------------
  273.             // Branch to use 2-term ascending series for small X and asymptotic
  274.             // form for large X when NB is not too large.
  275.             // ---------------------------------------------------------------------
  276.             double tempa;
  277.             double tempb;
  278.             double tempc;
  279.             double tover;
  280.             if (x < RTNSIG) {
  281.                 // ---------------------------------------------------------------------
  282.                 // Two-term ascending series for small X.
  283.                 // ---------------------------------------------------------------------
  284.                 tempa = 1;
  285.                 alpem = 1 + alpha;
  286.                 double halfx = 0;
  287.                 if (x > ENMTEN) {
  288.                     halfx = 0.5 * x;
  289.                 }
  290.                 if (alpha != 0) {
  291.                     tempa = FastMath.pow(halfx, alpha) /
  292.                             (alpha * Gamma.gamma(alpha));
  293.                 }
  294.                 tempb = 0;
  295.                 if (x + 1 > 1) {
  296.                     tempb = -halfx * halfx;
  297.                 }
  298.                 b[0] = tempa + (tempa * tempb / alpem);
  299.                 if ((x != 0) && (b[0] == 0)) {
  300.                     ncalc = 0;
  301.                 }
  302.                 if (nb != 1) {
  303.                     if (x <= 0) {
  304.                         for (int n = 1; n < nb; ++n) {
  305.                             b[n] = 0;
  306.                         }
  307.                     } else {
  308.                         // ---------------------------------------------------------------------
  309.                         // Calculate higher order functions.
  310.                         // ---------------------------------------------------------------------
  311.                         tempc = halfx;
  312.                         tover = tempb != 0 ? ENMTEN / tempb :  2 * ENMTEN / x;
  313.                         for (int n = 1; n < nb; ++n) {
  314.                             tempa /= alpem;
  315.                             alpem += 1;
  316.                             tempa *= tempc;
  317.                             if (tempa <= tover * alpem) {
  318.                                 tempa = 0;
  319.                             }
  320.                             b[n] = tempa + (tempa * tempb / alpem);
  321.                             if ((b[n] == 0) && (ncalc > n)) {
  322.                                 ncalc = n;
  323.                             }
  324.                         }
  325.                     }
  326.                 }
  327.             } else if ((x > 25.0) && (nb <= magx + 1)) {
  328.                 // ---------------------------------------------------------------------
  329.                 // Asymptotic series for X > 25
  330.                 // ---------------------------------------------------------------------
  331.                 final double xc = FastMath.sqrt(PI2 / x);
  332.                 final double mul = 0.125 / x;
  333.                 final double xin = mul * mul;
  334.                 final int m;
  335.                 if (x >= 130.0) {
  336.                     m = 4;
  337.                 } else if (x >= 35.0) {
  338.                     m = 8;
  339.                 } else {
  340.                     m = 11;
  341.                 }

  343.                 final double xm = 4.0 * m;
  344.                 // ---------------------------------------------------------------------
  345.                 // Argument reduction for SIN and COS routines.
  346.                 // ---------------------------------------------------------------------
  347.                 double t = ((int) ((x / TWOPI) + 0.5));
  348.                 final double z = x - t * TOWPI1 - t * TWOPI2 - (alpha + 0.5) / PI2;
  349.                 final SinCos vsc = FastMath.sinCos(z);
  350.                 double vsin = vsc.sin();
  351.                 double vcos = vsc.cos();
  352.                 double gnu = 2 * alpha;
  353.                 double capq;
  354.                 double capp;
  355.                 double s;
  356.                 double t1;
  357.                 double xk;
  358.                 for (int i = 1; i <= 2; i++) {
  359.                     s = (xm - 1 - gnu) * (xm - 1 + gnu) * xin * 0.5;
  360.                     t = (gnu - (xm - 3.0)) * (gnu + (xm - 3.0));
  361.                     capp = (s * t) / FACT[2 * m];
  362.                     t1 = (gnu - (xm + 1)) * (gnu + (xm + 1));
  363.                     capq = (s * t1) / FACT[2 * m + 1];
  364.                     xk = xm;
  365.                     int k = 2 * m;
  366.                     t1 = t;

  368.                     for (int j = 2; j <= m; j++) {
  369.                         xk -= 4.0;
  370.                         s = (xk - 1 - gnu) * (xk - 1 + gnu);
  371.                         t = (gnu - (xk - 3.0)) * (gnu + (xk - 3.0));
  372.                         capp = (capp + 1 / FACT[k - 2]) * s * t * xin;
  373.                         capq = (capq + 1 / FACT[k - 1]) * s * t1 * xin;
  374.                         k -= 2;
  375.                         t1 = t;
  376.                     }

  378.                     capp += 1;
  379.                     capq = (capq + 1) * ((gnu * gnu) - 1) * (0.125 / x);
  380.                     b[i - 1] = xc * (capp * vcos - capq * vsin);
  381.                     if (nb == 1) {
  382.                         return new BesselJResult(b.clone(), ncalc);
  383.                     }
  384.                     t = vsin;
  385.                     vsin = -vcos;
  386.                     vcos = t;
  387.                     gnu += 2.0;
  388.                 }

  390.                 // ---------------------------------------------------------------------
  391.                 // If NB > 2, compute J(X,ORDER+I) I = 2, NB-1
  392.                 // ---------------------------------------------------------------------
  393.                 if (nb > 2) {
  394.                     gnu = 2 * alpha + 2.0;
  395.                     for (int j = 2; j < nb; ++j) {
  396.                         b[j] = gnu * b[j - 1] / x - b[j - 2];
  397.                         gnu += 2.0;
  398.                     }
  399.                 }
  400.             } else {
  401.                 // ---------------------------------------------------------------------
  402.                 // Use recurrence to generate results. First initialize the
  403.                 // calculation of P*S.
  404.                 // ---------------------------------------------------------------------
  405.                 final int nbmx = nb - magx;
  406.                 int n = magx + 1;
  407.                 int nstart;
  408.                 int nend;
  409.                 double en = 2 * (n + alpha);
  410.                 double plast = 1;
  411.                 double p = en / x;
  412.                 double pold;
  413.                 // ---------------------------------------------------------------------
  414.                 // Calculate general significance test.
  415.                 // ---------------------------------------------------------------------
  416.                 double test = 2 * ENSIG;
  417.                 boolean readyToInitialize = false;
  418.                 if (nbmx >= 3) {
  419.                     // ---------------------------------------------------------------------
  420.                     // Calculate P*S until N = NB-1. Check for possible
  421.                     // overflow.
  422.                     // ---------------------------------------------------------------------
  423.                     tover = ENTEN / ENSIG;
  424.                     nstart = magx + 2;
  425.                     nend = nb - 1;
  426.                     en = 2 * (nstart - 1 + alpha);
  427.                     double psave;
  428.                     double psavel;
  429.                     for (int k = nstart; k <= nend; k++) {
  430.                         n = k;
  431.                         en += 2.0;
  432.                         pold = plast;
  433.                         plast = p;
  434.                         p = (en * plast / x) - pold;
  435.                         if (p > tover) {
  436.                             // ---------------------------------------------------------------------
  437.                             // To avoid overflow, divide P*S by TOVER. Calculate
  438.                             // P*S until
  439.                             // ABS(P) > 1.
  440.                             // ---------------------------------------------------------------------
  441.                             tover = ENTEN;
  442.                             p /= tover;
  443.                             plast /= tover;
  444.                             psave = p;
  445.                             psavel = plast;
  446.                             nstart = n + 1;
  447.                             do {
  448.                                 n += 1;
  449.                                 en += 2.0;
  450.                                 pold = plast;
  451.                                 plast = p;
  452.                                 p = (en * plast / x) - pold;
  453.                             } while (p <= 1);
  454.                             tempb = en / x;
  455.                             // ---------------------------------------------------------------------
  456.                             // Calculate backward test and find NCALC, the
  457.                             // highest N such that
  458.                             // the test is passed.
  459.                             // ---------------------------------------------------------------------
  460.                             test = pold * plast * (0.5 - 0.5 / (tempb * tempb));
  461.                             test /= ENSIG;
  462.                             p = plast * tover;
  463.                             n -= 1;
  464.                             en -= 2.0;
  465.                             nend = FastMath.min(nb, n);
  466.                             for (int l = nstart; l <= nend; l++) {
  467.                                 pold = psavel;
  468.                                 psavel = psave;
  469.                                 psave = (en * psavel / x) - pold;
  470.                                 if (psave * psavel > test) {
  471.                                     break;
  472.                                 }
  473.                             }
  474.                             ncalc = nend;
  475.                             readyToInitialize = true;
  476.                             break;
  477.                         }
  478.                     }
  479.                     if (!readyToInitialize) {
  480.                         n = nend;
  481.                         en = 2 * (n + alpha);
  482.                         // ---------------------------------------------------------------------
  483.                         // Calculate special significance test for NBMX > 2.
  484.                         // ---------------------------------------------------------------------
  485.                         test = FastMath.max(test, FastMath.sqrt(plast * ENSIG) *
  486.                                                   FastMath.sqrt(2 * p));
  487.                     }
  488.                 }
  489.                 // ---------------------------------------------------------------------
  490.                 // Calculate P*S until significance test passes.
  491.                 // ---------------------------------------------------------------------
  492.                 if (!readyToInitialize) {
  493.                     do {
  494.                         n += 1;
  495.                         en += 2.0;
  496.                         pold = plast;
  497.                         plast = p;
  498.                         p = (en * plast / x) - pold;
  499.                     } while (p < test);
  500.                 }
  501.                 // ---------------------------------------------------------------------
  502.                 // Initialize the backward recursion and the normalization sum.
  503.                 // ---------------------------------------------------------------------
  504.                 n += 1;
  505.                 en += 2.0;
  506.                 tempb = 0;
  507.                 tempa = 1 / p;
  508.                 int m = (2 * n) - 4 * (n / 2);
  509.                 double sum = 0;
  510.                 int em = n / 2;
  511.                 alpem = em - 1 + alpha;
  512.                 alp2em = 2 * em + alpha;
  513.                 if (m != 0) {
  514.                     sum = tempa * alpem * alp2em / em;
  515.                 }
  516.                 nend = n - nb;

  518.                 boolean readyToNormalize = false;
  519.                 boolean calculatedB0 = false;

  521.                 // ---------------------------------------------------------------------
  522.                 // Recur backward via difference equation, calculating (but not
  523.                 // storing) B(N), until N = NB.
  524.                 // ---------------------------------------------------------------------
  525.                 for (int l = 1; l <= nend; l++) {
  526.                     n -= 1;
  527.                     en -= 2.0;
  528.                     tempc = tempb;
  529.                     tempb = tempa;
  530.                     tempa = (en * tempb / x) - tempc;
  531.                     m = 2 - m;
  532.                     if (m != 0) {
  533.                         em -= 1;
  534.                         alp2em = 2 * em + alpha;
  535.                         if (n == 1) {
  536.                             break;
  537.                         }
  538.                         alpem = em - 1 + alpha;
  539.                         if (alpem == 0) {
  540.                             alpem = 1;
  541.                         }
  542.                         sum = (sum + tempa * alp2em) * alpem / em;
  543.                     }
  544.                 }

  546.                 // ---------------------------------------------------------------------
  547.                 // Store B(NB).
  548.                 // ---------------------------------------------------------------------
  549.                 b[n - 1] = tempa;
  550.                 if (nend >= 0) {
  551.                     if (nb <= 1) {
  552.                         alp2em = alpha;
  553.                         if (alpha + 1 == 1) {
  554.                             alp2em = 1;
  555.                         }
  556.                         sum += b[0] * alp2em;
  557.                         readyToNormalize = true;
  558.                     } else {
  559.                         // ---------------------------------------------------------------------
  560.                         // Calculate and store B(NB-1).
  561.                         // ---------------------------------------------------------------------
  562.                         n -= 1;
  563.                         en -= 2.0;
  564.                         b[n - 1] = (en * tempa / x) - tempb;
  565.                         if (n == 1) {
  566.                             calculatedB0 = true;
  567.                         } else {
  568.                             m = 2 - m;
  569.                             if (m != 0) {
  570.                                 em -= 1;
  571.                                 alp2em = 2 * em + alpha;
  572.                                 alpem = em - 1 + alpha;
  573.                                 if (alpem == 0) {
  574.                                     alpem = 1;
  575.                                 }

  577.                                 sum = (sum + (b[n - 1] * alp2em)) * alpem / em;
  578.                             }
  579.                         }
  580.                     }
  581.                 }
  582.                 if (!readyToNormalize && !calculatedB0) {
  583.                     nend = n - 2;
  584.                     if (nend != 0) {
  585.                         // ---------------------------------------------------------------------
  586.                         // Calculate via difference equation and store B(N),
  587.                         // until N = 2.
  588.                         // ---------------------------------------------------------------------

  590.                         for (int l = 1; l <= nend; l++) {
  591.                             n -= 1;
  592.                             en -= 2.0;
  593.                             b[n - 1] = (en * b[n] / x) - b[n + 1];
  594.                             m = 2 - m;
  595.                             if (m != 0) {
  596.                                 em -= 1;
  597.                                 alp2em = 2 * em + alpha;
  598.                                 alpem = em - 1 + alpha;
  599.                                 if (alpem == 0) {
  600.                                     alpem = 1;
  601.                                 }

  603.                                 sum = (sum + b[n - 1] * alp2em) * alpem / em;
  604.                             }
  605.                         }
  606.                     }
  607.                 }
  608.                 // ---------------------------------------------------------------------
  609.                 // Calculate b[0]
  610.                 // ---------------------------------------------------------------------
  611.                 if (!readyToNormalize) {
  612.                     if (!calculatedB0) {
  613.                         b[0] = 2.0 * (alpha + 1) * b[1] / x - b[2];
  614.                     }
  615.                     em -= 1;
  616.                     alp2em = 2 * em + alpha;
  617.                     if (alp2em == 0) {
  618.                         alp2em = 1;
  619.                     }
  620.                     sum += b[0] * alp2em;
  621.                 }
  622.                 // ---------------------------------------------------------------------
  623.                 // Normalize. Divide all B(N) by sum.
  624.                 // ---------------------------------------------------------------------

  626.                 if (FastMath.abs(alpha) > 1e-16) {
  627.                     sum *= Gamma.gamma(alpha) * FastMath.pow(x * 0.5, -alpha);
  628.                 }
  629.                 tempa = ENMTEN;
  630.                 if (sum > 1) {
  631.                     tempa *= sum;
  632.                 }

  634.                 for (n = 0; n < nb; n++) {
  635.                     if (FastMath.abs(b[n]) < tempa) {
  636.                         b[n] = 0;
  637.                     }
  638.                     b[n] /= sum;
  639.                 }
  640.             }
  641.             // ---------------------------------------------------------------------
  642.             // Error return -- X, NB, or ALPHA is out of range.
  643.             // ---------------------------------------------------------------------
  644.         } else {
  645.             if (b.length > 0) {
  646.                 b[0] = 0;
  647.             }
  648.             ncalc = FastMath.min(nb, 0) - 1;
  649.         }
  650.         return new BesselJResult(b.clone(), ncalc);
  651.     }
  652. }
Created with JaCoCo 0.8.12.202403310830