LegendreEllipticIntegral.java

  1. /*
  2.  * Licensed to the Hipparchus project 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 Hipparchus project 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.  */
  17. package org.hipparchus.special.elliptic.legendre;

  19. import java.util.function.DoubleFunction;
  20. import java.util.function.Function;

  22. import org.hipparchus.CalculusFieldElement;
  23. import org.hipparchus.analysis.CalculusFieldUnivariateFunction;
  24. import org.hipparchus.complex.Complex;
  25. import org.hipparchus.complex.ComplexUnivariateIntegrator;
  26. import org.hipparchus.complex.FieldComplex;
  27. import org.hipparchus.complex.FieldComplexUnivariateIntegrator;
  28. import org.hipparchus.special.elliptic.carlson.CarlsonEllipticIntegral;
  29. import org.hipparchus.util.FastMath;
  30. import org.hipparchus.util.MathUtils;

  32. /** Complete and incomplete elliptic integrals in Legendre form.
  33.  * <p>
  34.  * The elliptic integrals are related to Jacobi elliptic functions.
  35.  * </p>
  36.  * <p>
  37.  * <em>
  38.  * Beware that when computing elliptic integrals in the complex plane,
  39.  * many issues arise due to branch cuts. See the
  40.  * <a href="https://www.hipparchus.org/hipparchus-core/special.html#Elliptic_functions_and_integrals">user guide</a>
  41.  * for a thorough explanation.
  42.  * </em>
  43.  * </p>
  44.  * <p>
  45.  * There are different conventions to interpret the arguments of
  46.  * Legendre elliptic integrals. In mathematical texts, these conventions show
  47.  * up using the separator between arguments. So for example for the incomplete
  48.  * integral of the first kind F we have:
  49.  * </p>
  50.  * <ul>
  51.  *   <li>F(φ, k): the first argument φ is an angle and the second argument k
  52.  *       is the elliptic modulus: this is the trigonometric form of the integral</li>
  53.  *   <li>F(φ; m): the first argument φ is an angle and the second argument m=k²
  54.  *       is the parameter: this is also a trigonometric form of the integral</li>
  55.  *   <li>F(x|m): the first argument x=sin(φ) is not an angle anymore and the
  56.  *       second argument m=k² is the parameter: this is the Legendre form</li>
  57.  *   <li>F(φ\α): the first argument φ is an angle and the second argument α is the
  58.  *       modular angle</li>
  59.  * </ul>
  60.  * <p>
  61.  * As we have no separator in a method call, we have to adopt one convention
  62.  * and stick to it. In Hipparchus, we adopted the Legendre form (i.e. F(x|m),
  63.  * with x=sin(φ) and m=k². These conventions are consistent with Wolfram Alpha
  64.  * functions EllipticF, EllipticE, ElliptiPI…
  65.  * </p>
  66.  * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  67.  * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheFirstKind.html">Complete Elliptic Integrals of the First Kind (MathWorld)</a>
  68.  * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheSecondKind.html">Complete Elliptic Integrals of the Second Kind (MathWorld)</a>
  69.  * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheFirstKind.html">Elliptic Integrals of the First Kind (MathWorld)</a>
  70.  * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheSecondKind.html">Elliptic Integrals of the Second Kind (MathWorld)</a>
  71.  * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheThirdKind.html">Elliptic Integrals of the Third Kind (MathWorld)</a>
  72.  * @since 2.0
  73.  */
  74. public class LegendreEllipticIntegral { // NOPMD - this class has a high number of methods, it is normal

  76.     /** Private constructor for a utility class.
  77.      */
  78.     private LegendreEllipticIntegral() {
  79.         // nothing to do
  80.     }

  82.     /** Get the nome q.
  83.      * @param m parameter (m=k² where k is the elliptic modulus)
  84.      * @return nome q
  85.      */
  86.     public static double nome(final double m) {
  87.         if (m < 1.0e-16) {
  88.             // first terms of infinite series in Abramowitz and Stegun 17.3.21
  89.             final double m16 = m * 0.0625;
  90.             return m16 * (1 + 8 * m16);
  91.         } else {
  92.             return FastMath.exp(-FastMath.PI * bigKPrime(m) / bigK(m));
  93.         }
  94.     }

  96.     /** Get the nome q.
  97.      * @param m parameter (m=k² where k is the elliptic modulus)
  98.      * @param <T> the type of the field elements
  99.      * @return nome q
  100.      */
  101.     public static <T extends CalculusFieldElement<T>> T nome(final T m) {
  102.         final T one = m.getField().getOne();
  103.         if (m.norm() < 100 * one.ulp().getReal()) {
  104.             // first terms of infinite series in Abramowitz and Stegun 17.3.21
  105.             final T m16 = m.multiply(0.0625);
  106.             return m16.multiply(m16.multiply(8).add(1));
  107.         } else {
  108.             return FastMath.exp(bigKPrime(m).divide(bigK(m)).multiply(one.getPi().negate()));
  109.         }
  110.     }

  112.     /** Get the complete elliptic integral of the first kind K(m).
  113.      * <p>
  114.      * The complete elliptic integral of the first kind K(m) is
  115.      * \[
  116.      *    \int_0^{\frac{\pi}{2}} \frac{d\theta}{\sqrt{1-m \sin^2\theta}}
  117.      * \]
  118.      * it corresponds to the real quarter-period of Jacobi elliptic functions
  119.      * </p>
  120.      * <p>
  121.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  122.      * Carlson elliptic integrals}.
  123.      * </p>
  124.      * @param m parameter (m=k² where k is the elliptic modulus)
  125.      * @return complete elliptic integral of the first kind K(m)
  126.      * @see #bigKPrime(double)
  127.      * @see #bigF(double, double)
  128.      * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheFirstKind.html">Complete Elliptic Integrals of the First Kind (MathWorld)</a>
  129.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  130.      */
  131.     public static double bigK(final double m) {
  132.         if (m < 1.0e-8) {
  133.             // first terms of infinite series in Abramowitz and Stegun 17.3.11
  134.             return (1 + 0.25 * m) * MathUtils.SEMI_PI;
  135.         } else {
  136.             return CarlsonEllipticIntegral.rF(0, 1.0 - m, 1);
  137.         }
  138.     }

  140.     /** Get the complete elliptic integral of the first kind K(m).
  141.      * <p>
  142.      * The complete elliptic integral of the first kind K(m) is
  143.      * \[
  144.      *    \int_0^{\frac{\pi}{2}} \frac{d\theta}{\sqrt{1-m \sin^2\theta}}
  145.      * \]
  146.      * it corresponds to the real quarter-period of Jacobi elliptic functions
  147.      * </p>
  148.      * <p>
  149.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  150.      * Carlson elliptic integrals}.
  151.      * </p>
  152.      * @param m parameter (m=k² where k is the elliptic modulus)
  153.      * @param <T> the type of the field elements
  154.      * @return complete elliptic integral of the first kind K(m)
  155.      * @see #bigKPrime(CalculusFieldElement)
  156.      * @see #bigF(CalculusFieldElement, CalculusFieldElement)
  157.      * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheFirstKind.html">Complete Elliptic Integrals of the First Kind (MathWorld)</a>
  158.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  159.      */
  160.     public static <T extends CalculusFieldElement<T>> T bigK(final T m) {
  161.         final T zero = m.getField().getZero();
  162.         final T one  = m.getField().getOne();
  163.         if (m.norm() < 1.0e7 * one.ulp().getReal()) {

  165.             // first terms of infinite series in Abramowitz and Stegun 17.3.11
  166.             return one.add(m.multiply(0.25)).multiply(zero.getPi().multiply(0.5));

  168.         } else {
  169.             return CarlsonEllipticIntegral.rF(zero, one.subtract(m), one);
  170.         }
  171.     }

  173.     /** Get the complete elliptic integral of the first kind K(m).
  174.      * <p>
  175.      * The complete elliptic integral of the first kind K(m) is
  176.      * \[
  177.      *    \int_0^{\frac{\pi}{2}} \frac{d\theta}{\sqrt{1-m \sin^2\theta}}
  178.      * \]
  179.      * it corresponds to the real quarter-period of Jacobi elliptic functions
  180.      * </p>
  181.      * <p>
  182.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  183.      * Carlson elliptic integrals}.
  184.      * </p>
  185.      * @param m parameter (m=k² where k is the elliptic modulus)
  186.      * @return complete elliptic integral of the first kind K(m)
  187.      * @see #bigKPrime(Complex)
  188.      * @see #bigF(Complex, Complex)
  189.      * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheFirstKind.html">Complete Elliptic Integrals of the First Kind (MathWorld)</a>
  190.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  191.      */
  192.     public static Complex bigK(final Complex m) {
  193.         if (m.norm() < 1.0e-8) {
  194.             // first terms of infinite series in Abramowitz and Stegun 17.3.11
  195.             return Complex.ONE.add(m.multiply(0.25)).multiply(MathUtils.SEMI_PI);
  196.         } else {
  197.             return CarlsonEllipticIntegral.rF(Complex.ZERO, Complex.ONE.subtract(m), Complex.ONE);
  198.         }
  199.     }

  201.     /** Get the complete elliptic integral of the first kind K(m).
  202.      * <p>
  203.      * The complete elliptic integral of the first kind K(m) is
  204.      * \[
  205.      *    \int_0^{\frac{\pi}{2}} \frac{d\theta}{\sqrt{1-m \sin^2\theta}}
  206.      * \]
  207.      * it corresponds to the real quarter-period of Jacobi elliptic functions
  208.      * </p>
  209.      * <p>
  210.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  211.      * Carlson elliptic integrals}.
  212.      * </p>
  213.      * @param m parameter (m=k² where k is the elliptic modulus)
  214.      * @param <T> the type of the field elements
  215.      * @return complete elliptic integral of the first kind K(m)
  216.      * @see #bigKPrime(FieldComplex)
  217.      * @see #bigF(FieldComplex, FieldComplex)
  218.      * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheFirstKind.html">Complete Elliptic Integrals of the First Kind (MathWorld)</a>
  219.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  220.      */
  221.     public static <T extends CalculusFieldElement<T>> FieldComplex<T> bigK(final FieldComplex<T> m) {
  222.         final FieldComplex<T> zero = m.getField().getZero();
  223.         final FieldComplex<T> one  = m.getField().getOne();
  224.         if (m.norm() < 1.0e7 * one.ulp().getReal()) {

  226.             // first terms of infinite series in Abramowitz and Stegun 17.3.11
  227.             return one.add(m.multiply(0.25)).multiply(zero.getPi().multiply(0.5));

  229.         } else {
  230.             return CarlsonEllipticIntegral.rF(zero, one.subtract(m), one);
  231.         }
  232.     }

  234.     /** Get the complete elliptic integral of the first kind K'(m).
  235.      * <p>
  236.      * The complete elliptic integral of the first kind K'(m) is
  237.      * \[
  238.      *    \int_0^{\frac{\pi}{2}} \frac{d\theta}{\sqrt{1-(1-m) \sin^2\theta}}
  239.      * \]
  240.      * it corresponds to the imaginary quarter-period of Jacobi elliptic functions
  241.      * </p>
  242.      * <p>
  243.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  244.      * Carlson elliptic integrals}.
  245.      * </p>
  246.      * @param m parameter (m=k² where k is the elliptic modulus)
  247.      * @return complete elliptic integral of the first kind K'(m)
  248.      * @see #bigK(double)
  249.      * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheFirstKind.html">Complete Elliptic Integrals of the First Kind (MathWorld)</a>
  250.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  251.      */
  252.     public static double bigKPrime(final double m) {
  253.         return CarlsonEllipticIntegral.rF(0, m, 1);
  254.     }

  256.     /** Get the complete elliptic integral of the first kind K'(m).
  257.      * <p>
  258.      * The complete elliptic integral of the first kind K'(m) is
  259.      * \[
  260.      *    \int_0^{\frac{\pi}{2}} \frac{d\theta}{\sqrt{1-(1-m) \sin^2\theta}}
  261.      * \]
  262.      * it corresponds to the imaginary quarter-period of Jacobi elliptic functions
  263.      * </p>
  264.      * <p>
  265.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  266.      * Carlson elliptic integrals}.
  267.      * </p>
  268.      * @param m parameter (m=k² where k is the elliptic modulus)
  269.      * @param <T> the type of the field elements
  270.      * @return complete elliptic integral of the first kind K'(m)
  271.      * @see #bigK(CalculusFieldElement)
  272.      * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheFirstKind.html">Complete Elliptic Integrals of the First Kind (MathWorld)</a>
  273.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  274.      */
  275.     public static <T extends CalculusFieldElement<T>> T bigKPrime(final T m) {
  276.         final T zero = m.getField().getZero();
  277.         final T one  = m.getField().getOne();
  278.         return CarlsonEllipticIntegral.rF(zero, m, one);
  279.     }

  281.     /** Get the complete elliptic integral of the first kind K'(m).
  282.      * <p>
  283.      * The complete elliptic integral of the first kind K'(m) is
  284.      * \[
  285.      *    \int_0^{\frac{\pi}{2}} \frac{d\theta}{\sqrt{1-(1-m) \sin^2\theta}}
  286.      * \]
  287.      * it corresponds to the imaginary quarter-period of Jacobi elliptic functions
  288.      * </p>
  289.      * <p>
  290.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  291.      * Carlson elliptic integrals}.
  292.      * </p>
  293.      * @param m parameter (m=k² where k is the elliptic modulus)
  294.      * @return complete elliptic integral of the first kind K'(m)
  295.      * @see #bigK(Complex)
  296.      * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheFirstKind.html">Complete Elliptic Integrals of the First Kind (MathWorld)</a>
  297.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  298.      */
  299.     public static Complex bigKPrime(final Complex m) {
  300.         return CarlsonEllipticIntegral.rF(Complex.ZERO, m, Complex.ONE);
  301.     }

  303.     /** Get the complete elliptic integral of the first kind K'(m).
  304.      * <p>
  305.      * The complete elliptic integral of the first kind K'(m) is
  306.      * \[
  307.      *    \int_0^{\frac{\pi}{2}} \frac{d\theta}{\sqrt{1-(1-m) \sin^2\theta}}
  308.      * \]
  309.      * it corresponds to the imaginary quarter-period of Jacobi elliptic functions
  310.      * </p>
  311.      * <p>
  312.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  313.      * Carlson elliptic integrals}.
  314.      * </p>
  315.      * @param m parameter (m=k² where k is the elliptic modulus)
  316.      * @param <T> the type of the field elements
  317.      * @return complete elliptic integral of the first kind K'(m)
  318.      * @see #bigK(FieldComplex)
  319.      * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheFirstKind.html">Complete Elliptic Integrals of the First Kind (MathWorld)</a>
  320.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  321.      */
  322.     public static <T extends CalculusFieldElement<T>> FieldComplex<T> bigKPrime(final FieldComplex<T> m) {
  323.         final FieldComplex<T> zero = m.getField().getZero();
  324.         final FieldComplex<T> one  = m.getField().getOne();
  325.         return CarlsonEllipticIntegral.rF(zero, m, one);
  326.     }

  328.     /** Get the complete elliptic integral of the second kind E(m).
  329.      * <p>
  330.      * The complete elliptic integral of the second kind E(m) is
  331.      * \[
  332.      *    \int_0^{\frac{\pi}{2}} \sqrt{1-m \sin^2\theta} d\theta
  333.      * \]
  334.      * </p>
  335.      * <p>
  336.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  337.      * Carlson elliptic integrals}.
  338.      * </p>
  339.      * @param m parameter (m=k² where k is the elliptic modulus)
  340.      * @return complete elliptic integral of the second kind E(m)
  341.      * @see #bigE(double, double)
  342.      * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheSecondKind.html">Complete Elliptic Integrals of the Second Kind (MathWorld)</a>
  343.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  344.      */
  345.     public static double bigE(final double m) {
  346.         return CarlsonEllipticIntegral.rG(0, 1 - m, 1) * 2;
  347.     }

  349.     /** Get the complete elliptic integral of the second kind E(m).
  350.      * <p>
  351.      * The complete elliptic integral of the second kind E(m) is
  352.      * \[
  353.      *    \int_0^{\frac{\pi}{2}} \sqrt{1-m \sin^2\theta} d\theta
  354.      * \]
  355.      * </p>
  356.      * <p>
  357.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  358.      * Carlson elliptic integrals}.
  359.      * </p>
  360.      * @param m parameter (m=k² where k is the elliptic modulus)
  361.      * @param <T> the type of the field elements
  362.      * @return complete elliptic integral of the second kind E(m)
  363.      * @see #bigE(CalculusFieldElement, CalculusFieldElement)
  364.      * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheSecondKind.html">Complete Elliptic Integrals of the Second Kind (MathWorld)</a>
  365.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  366.      */
  367.     public static <T extends CalculusFieldElement<T>> T bigE(final T m) {
  368.         final T zero = m.getField().getZero();
  369.         final T one  = m.getField().getOne();
  370.         return CarlsonEllipticIntegral.rG(zero, one.subtract(m), one).multiply(2);
  371.     }

  373.     /** Get the complete elliptic integral of the second kind E(m).
  374.      * <p>
  375.      * The complete elliptic integral of the second kind E(m) is
  376.      * \[
  377.      *    \int_0^{\frac{\pi}{2}} \sqrt{1-m \sin^2\theta} d\theta
  378.      * \]
  379.      * </p>
  380.      * <p>
  381.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  382.      * Carlson elliptic integrals}.
  383.      * </p>
  384.      * @param m parameter (m=k² where k is the elliptic modulus)
  385.      * @return complete elliptic integral of the second kind E(m)
  386.      * @see #bigE(Complex, Complex)
  387.      * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheSecondKind.html">Complete Elliptic Integrals of the Second Kind (MathWorld)</a>
  388.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  389.      */
  390.     public static Complex bigE(final Complex m) {
  391.         return CarlsonEllipticIntegral.rG(Complex.ZERO,
  392.                                           Complex.ONE.subtract(m),
  393.                                           Complex.ONE).multiply(2);
  394.     }

  396.     /** Get the complete elliptic integral of the second kind E(m).
  397.      * <p>
  398.      * The complete elliptic integral of the second kind E(m) is
  399.      * \[
  400.      *    \int_0^{\frac{\pi}{2}} \sqrt{1-m \sin^2\theta} d\theta
  401.      * \]
  402.      * </p>
  403.      * <p>
  404.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  405.      * Carlson elliptic integrals}.
  406.      * </p>
  407.      * @param m parameter (m=k² where k is the elliptic modulus)
  408.      * @param <T> the type of the field elements
  409.      * @return complete elliptic integral of the second kind E(m)
  410.      * @see #bigE(FieldComplex, FieldComplex)
  411.      * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheSecondKind.html">Complete Elliptic Integrals of the Second Kind (MathWorld)</a>
  412.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  413.      */
  414.     public static <T extends CalculusFieldElement<T>> FieldComplex<T> bigE(final FieldComplex<T> m) {
  415.         final FieldComplex<T> zero = m.getField().getZero();
  416.         final FieldComplex<T> one  = m.getField().getOne();
  417.         return CarlsonEllipticIntegral.rG(zero, one.subtract(m), one).multiply(2);
  418.     }

  420.     /** Get the complete elliptic integral D(m) = [K(m) - E(m)]/m.
  421.      * <p>
  422.      * The complete elliptic integral D(m) is
  423.      * \[
  424.      *    \int_0^{\frac{\pi}{2}} \frac{\sin^2\theta}{\sqrt{1-m \sin^2\theta}} d\theta
  425.      * \]
  426.      * </p>
  427.      * <p>
  428.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  429.      * Carlson elliptic integrals}.
  430.      * </p>
  431.      * @param m parameter (m=k² where k is the elliptic modulus)
  432.      * @return complete elliptic integral D(m)
  433.      * @see #bigD(double, double)
  434.      */
  435.     public static double bigD(final double m) {
  436.         return CarlsonEllipticIntegral.rD(0, 1 - m, 1) / 3;
  437.     }

  439.     /** Get the complete elliptic integral D(m) = [K(m) - E(m)]/m.
  440.      * <p>
  441.      * The complete elliptic integral D(m) is
  442.      * \[
  443.      *    \int_0^{\frac{\pi}{2}} \frac{\sin^2\theta}{\sqrt{1-m \sin^2\theta}} d\theta
  444.      * \]
  445.      * </p>
  446.      * <p>
  447.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  448.      * Carlson elliptic integrals}.
  449.      * </p>
  450.      * @param m parameter (m=k² where k is the elliptic modulus)
  451.      * @param <T> the type of the field elements
  452.      * @return complete elliptic integral D(m)
  453.      * @see #bigD(CalculusFieldElement, CalculusFieldElement)
  454.      */
  455.     public static <T extends CalculusFieldElement<T>> T bigD(final T m) {
  456.         final T zero = m.getField().getZero();
  457.         final T one  = m.getField().getOne();
  458.         return CarlsonEllipticIntegral.rD(zero, one.subtract(m), one).divide(3);
  459.     }

  461.     /** Get the complete elliptic integral D(m) = [K(m) - E(m)]/m.
  462.      * <p>
  463.      * The complete elliptic integral D(m) is
  464.      * \[
  465.      *    \int_0^{\frac{\pi}{2}} \frac{\sin^2\theta}{\sqrt{1-m \sin^2\theta}} d\theta
  466.      * \]
  467.      * </p>
  468.      * <p>
  469.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  470.      * Carlson elliptic integrals}.
  471.      * </p>
  472.      * @param m parameter (m=k² where k is the elliptic modulus)
  473.      * @return complete elliptic integral D(m)
  474.      * @see #bigD(Complex, Complex)
  475.      */
  476.     public static Complex bigD(final Complex m) {
  477.         return CarlsonEllipticIntegral.rD(Complex.ZERO, Complex.ONE.subtract(m), Complex.ONE).divide(3);
  478.     }

  480.     /** Get the complete elliptic integral D(m) = [K(m) - E(m)]/m.
  481.      * <p>
  482.      * The complete elliptic integral D(m) is
  483.      * \[
  484.      *    \int_0^{\frac{\pi}{2}} \frac{\sin^2\theta}{\sqrt{1-m \sin^2\theta}} d\theta
  485.      * \]
  486.      * </p>
  487.      * <p>
  488.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  489.      * Carlson elliptic integrals}.
  490.      * </p>
  491.      * @param m parameter (m=k² where k is the elliptic modulus)
  492.      * @param <T> the type of the field elements
  493.      * @return complete elliptic integral D(m)
  494.      * @see #bigD(FieldComplex, FieldComplex)
  495.      */
  496.     public static <T extends CalculusFieldElement<T>> FieldComplex<T> bigD(final FieldComplex<T> m) {
  497.         final FieldComplex<T> zero = m.getField().getZero();
  498.         final FieldComplex<T> one  = m.getField().getOne();
  499.         return CarlsonEllipticIntegral.rD(zero, one.subtract(m), one).divide(3);
  500.     }

  502.     /** Get the complete elliptic integral of the third kind Π(n, m).
  503.      * <p>
  504.      * The complete elliptic integral of the third kind Π(n, m) is
  505.      * \[
  506.      *    \int_0^{\frac{\pi}{2}} \frac{d\theta}{\sqrt{1-m \sin^2\theta}(1-n \sin^2\theta)}
  507.      * \]
  508.      * </p>
  509.      * <p>
  510.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  511.      * Carlson elliptic integrals}.
  512.      * </p>
  513.      * @param n elliptic characteristic
  514.      * @param m parameter (m=k² where k is the elliptic modulus)
  515.      * @return complete elliptic integral of the third kind Π(n, m)
  516.      * @see #bigPi(double, double, double)
  517.      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheThirdKind.html">Elliptic Integrals of the Third Kind (MathWorld)</a>
  518.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  519.      */
  520.     public static double bigPi(final double n, final double m) {
  521.         final double kPrime2 = 1 - m;
  522.         final double delta   = n * (m - n) * (n - 1);
  523.         return CarlsonEllipticIntegral.rF(0, kPrime2, 1) +
  524.                CarlsonEllipticIntegral.rJ(0, kPrime2, 1, 1 - n, delta) * n / 3;
  525.     }

  527.     /** Get the complete elliptic integral of the third kind Π(n, m).
  528.      * <p>
  529.      * The complete elliptic integral of the third kind Π(n, m) is
  530.      * \[
  531.      *    \int_0^{\frac{\pi}{2}} \frac{d\theta}{\sqrt{1-m \sin^2\theta}(1-n \sin^2\theta)}
  532.      * \]
  533.      * </p>
  534.      * <p>
  535.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  536.      * Carlson elliptic integrals}.
  537.      * </p>
  538.      * @param n elliptic characteristic
  539.      * @param m parameter (m=k² where k is the elliptic modulus)
  540.      * @param <T> the type of the field elements
  541.      * @return complete elliptic integral of the third kind Π(n, m)
  542.      * @see #bigPi(CalculusFieldElement, CalculusFieldElement, CalculusFieldElement)
  543.      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheThirdKind.html">Elliptic Integrals of the Third Kind (MathWorld)</a>
  544.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  545.      */
  546.     public static <T extends CalculusFieldElement<T>> T bigPi(final T n, final T m) {
  547.         final T zero    = m.getField().getZero();
  548.         final T one     = m.getField().getOne();
  549.         final T kPrime2 = one.subtract(m);
  550.         final T delta   = n.multiply(m.subtract(n)).multiply(n.subtract(1));
  551.         return CarlsonEllipticIntegral.rF(zero, kPrime2, one).
  552.                add(CarlsonEllipticIntegral.rJ(zero, kPrime2, one, one.subtract(n), delta).multiply(n).divide(3));
  553.     }

  555.     /** Get the complete elliptic integral of the third kind Π(n, m).
  556.      * <p>
  557.      * The complete elliptic integral of the third kind Π(n, m) is
  558.      * \[
  559.      *    \int_0^{\frac{\pi}{2}} \frac{d\theta}{\sqrt{1-m \sin^2\theta}(1-n \sin^2\theta)}
  560.      * \]
  561.      * </p>
  562.      * <p>
  563.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  564.      * Carlson elliptic integrals}.
  565.      * </p>
  566.      * @param n elliptic characteristic
  567.      * @param m parameter (m=k² where k is the elliptic modulus)
  568.      * @return complete elliptic integral of the third kind Π(n, m)
  569.      * @see #bigPi(Complex, Complex, Complex)
  570.      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheThirdKind.html">Elliptic Integrals of the Third Kind (MathWorld)</a>
  571.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  572.      */
  573.     public static Complex bigPi(final Complex n, final Complex m) {
  574.         final Complex kPrime2 = Complex.ONE.subtract(m);
  575.         final Complex delta   = n.multiply(m.subtract(n)).multiply(n.subtract(1));
  576.         return CarlsonEllipticIntegral.rF(Complex.ZERO, kPrime2, Complex.ONE).
  577.                add(CarlsonEllipticIntegral.rJ(Complex.ZERO, kPrime2, Complex.ONE, Complex.ONE.subtract(n), delta).multiply(n).divide(3));
  578.     }

  580.     /** Get the complete elliptic integral of the third kind Π(n, m).
  581.      * <p>
  582.      * The complete elliptic integral of the third kind Π(n, m) is
  583.      * \[
  584.      *    \int_0^{\frac{\pi}{2}} \frac{d\theta}{\sqrt{1-m \sin^2\theta}(1-n \sin^2\theta)}
  585.      * \]
  586.      * </p>
  587.      * <p>
  588.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  589.      * Carlson elliptic integrals}.
  590.      * </p>
  591.      * @param n elliptic characteristic
  592.      * @param m parameter (m=k² where k is the elliptic modulus)
  593.      * @param <T> the type of the field elements
  594.      * @return complete elliptic integral of the third kind Π(n, m)
  595.      * @see #bigPi(FieldComplex, FieldComplex, FieldComplex)
  596.      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheThirdKind.html">Elliptic Integrals of the Third Kind (MathWorld)</a>
  597.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  598.      */
  599.     public static <T extends CalculusFieldElement<T>> FieldComplex<T> bigPi(final FieldComplex<T> n, final FieldComplex<T> m) {
  600.         final FieldComplex<T> zero    = m.getField().getZero();
  601.         final FieldComplex<T> one     = m.getField().getOne();
  602.         final FieldComplex<T> kPrime2 = one.subtract(m);
  603.         final FieldComplex<T> delta   = n.multiply(m.subtract(n)).multiply(n.subtract(1));
  604.         return CarlsonEllipticIntegral.rF(zero, kPrime2, one).
  605.                add(CarlsonEllipticIntegral.rJ(zero, kPrime2, one, one.subtract(n), delta).multiply(n).divide(3));
  606.     }

  608.     /** Get the incomplete elliptic integral of the first kind F(φ, m).
  609.      * <p>
  610.      * The incomplete elliptic integral of the first kind F(φ, m) is
  611.      * \[
  612.      *    \int_0^{\phi} \frac{d\theta}{\sqrt{1-m \sin^2\theta}}
  613.      * \]
  614.      * </p>
  615.      * <p>
  616.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  617.      * Carlson elliptic integrals}.
  618.      * </p>
  619.      * @param phi amplitude (i.e. upper bound of the integral)
  620.      * @param m parameter (m=k² where k is the elliptic modulus)
  621.      * @return incomplete elliptic integral of the first kind F(φ, m)
  622.      * @see #bigK(double)
  623.      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheFirstKind.html">Elliptic Integrals of the First Kind (MathWorld)</a>
  624.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  625.      */
  626.     public static double bigF(final double phi, final double m) {

  628.         // argument reduction
  629.         final DoubleArgumentReduction ar = new DoubleArgumentReduction(phi, m, LegendreEllipticIntegral::bigK);

  631.         // integrate part between 0 and π/2
  632.         final double cM1 = ar.csc2 - 1.0;
  633.         final double cMm = ar.csc2 - m;
  634.         final double incomplete =  CarlsonEllipticIntegral.rF(cM1, cMm, ar.csc2);

  636.         // combine complete and incomplete parts
  637.         return ar.negate ? ar.complete - incomplete : ar.complete + incomplete;

  639.     }

  641.     /** Get the incomplete elliptic integral of the first kind F(φ, m).
  642.      * <p>
  643.      * The incomplete elliptic integral of the first kind F(φ, m) is
  644.      * \[
  645.      *    \int_0^{\phi} \frac{d\theta}{\sqrt{1-m \sin^2\theta}}
  646.      * \]
  647.      * </p>
  648.      * <p>
  649.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  650.      * Carlson elliptic integrals}.
  651.      * </p>
  652.      * @param phi amplitude (i.e. upper bound of the integral)
  653.      * @param m parameter (m=k² where k is the elliptic modulus)
  654.      * @param <T> the type of the field elements
  655.      * @return incomplete elliptic integral of the first kind F(φ, m)
  656.      * @see #bigK(CalculusFieldElement)
  657.      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheFirstKind.html">Elliptic Integrals of the First Kind (MathWorld)</a>
  658.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  659.      */
  660.     public static <T extends CalculusFieldElement<T>> T bigF(final T phi, final T m) {

  662.         // argument reduction
  663.         final FieldArgumentReduction<T> ar = new FieldArgumentReduction<>(phi, m, LegendreEllipticIntegral::bigK);

  665.         // integrate part between 0 and π/2
  666.         final T cM1        = ar.csc2.subtract(1);
  667.         final T cMm        = ar.csc2.subtract(m);
  668.         final T incomplete = CarlsonEllipticIntegral.rF(cM1, cMm, ar.csc2);

  670.         // combine complete and incomplete parts
  671.         return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);

  673.     }

  675.     /** Get the incomplete elliptic integral of the first kind F(φ, m).
  676.      * <p>
  677.      * <em>
  678.      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
  679.      * are considered experimental for now, they have known issues.
  680.      * </em>
  681.      * </p>
  682.      * <p>
  683.      * The incomplete elliptic integral of the first kind F(φ, m) is
  684.      * \[
  685.      *    \int_0^{\phi} \frac{d\theta}{\sqrt{1-m \sin^2\theta}}
  686.      * \]
  687.      * </p>
  688.      * <p>
  689.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  690.      * Carlson elliptic integrals}.
  691.      * </p>
  692.      * @param phi amplitude (i.e. upper bound of the integral)
  693.      * @param m parameter (m=k² where k is the elliptic modulus)
  694.      * @return incomplete elliptic integral of the first kind F(φ, m)
  695.      * @see #bigK(Complex)
  696.      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheFirstKind.html">Elliptic Integrals of the First Kind (MathWorld)</a>
  697.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  698.      */
  699.     public static Complex bigF(final Complex phi, final Complex m) {

  701.         // argument reduction
  702.         final FieldArgumentReduction<Complex> ar = new FieldArgumentReduction<>(phi, m, LegendreEllipticIntegral::bigK);

  704.         // integrate part between 0 and π/2
  705.         final Complex cM1        = ar.csc2.subtract(1);
  706.         final Complex cMm        = ar.csc2.subtract(m);
  707.         final Complex incomplete = CarlsonEllipticIntegral.rF(cM1, cMm, ar.csc2);

  709.         // combine complete and incomplete parts
  710.         return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);

  712.     }

  714.     /** Get the incomplete elliptic integral of the first kind F(φ, m) using numerical integration.
  715.      * <p>
  716.      * <em>
  717.      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
  718.      * are considered experimental for now, they have known issues.
  719.      * </em>
  720.      * </p>
  721.      * <p>
  722.      * The incomplete elliptic integral of the first kind F(φ, m) is
  723.      * \[
  724.      *    \int_0^{\phi} \frac{d\theta}{\sqrt{1-m \sin^2\theta}}
  725.      * \]
  726.      * </p>
  727.      * <p>
  728.      * The algorithm for evaluating the functions is based on numerical integration.
  729.      * If integration path comes too close to a pole of the integrand, then integration will fail
  730.      * with a {@link org.hipparchus.exception.MathIllegalStateException MathIllegalStateException}
  731.      * even for very large {@code maxEval}. This is normal behavior.
  732.      * </p>
  733.      * @param phi amplitude (i.e. upper bound of the integral)
  734.      * @param m parameter (m=k² where k is the elliptic modulus)
  735.      * @param integrator integrator to use
  736.      * @param maxEval maximum number of evaluations (real and imaginary
  737.      * parts are evaluated separately, so up to twice this number may be used)
  738.      * @return incomplete elliptic integral of the first kind F(φ, m)
  739.      * @see #bigK(Complex)
  740.      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheFirstKind.html">Elliptic Integrals of the First Kind (MathWorld)</a>
  741.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  742.      */
  743.     public static Complex bigF(final Complex phi, final Complex m,
  744.                                final ComplexUnivariateIntegrator integrator, final int maxEval) {
  745.         return integrator.integrate(maxEval, new First<>(m), phi.getField().getZero(), phi);
  746.     }

  748.     /** Get the incomplete elliptic integral of the first kind F(φ, m).
  749.      * <p>
  750.      * <em>
  751.      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
  752.      * are considered experimental for now, they have known issues.
  753.      * </em>
  754.      * </p>
  755.      * <p>
  756.      * The incomplete elliptic integral of the first kind F(φ, m) is
  757.      * \[
  758.      *    \int_0^{\phi} \frac{d\theta}{\sqrt{1-m \sin^2\theta}}
  759.      * \]
  760.      * </p>
  761.      * <p>
  762.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  763.      * Carlson elliptic integrals}.
  764.      * </p>
  765.      * @param phi amplitude (i.e. upper bound of the integral)
  766.      * @param m parameter (m=k² where k is the elliptic modulus)
  767.      * @param <T> the type of the field elements
  768.      * @return incomplete elliptic integral of the first kind F(φ, m)
  769.      * @see #bigK(CalculusFieldElement)
  770.      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheFirstKind.html">Elliptic Integrals of the First Kind (MathWorld)</a>
  771.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  772.      */
  773.     public static <T extends CalculusFieldElement<T>> FieldComplex<T> bigF(final FieldComplex<T> phi, final FieldComplex<T> m) {

  775.         // argument reduction
  776.         final FieldArgumentReduction<FieldComplex<T>> ar = new FieldArgumentReduction<>(phi, m, LegendreEllipticIntegral::bigK);

  778.         // integrate part between 0 and π/2
  779.         final FieldComplex<T> cM1        = ar.csc2.subtract(1);
  780.         final FieldComplex<T> cMm        = ar.csc2.subtract(m);
  781.         final FieldComplex<T> incomplete = CarlsonEllipticIntegral.rF(cM1, cMm, ar.csc2);

  783.         // combine complete and incomplete parts
  784.         return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);

  786.     }

  788.     /** Get the incomplete elliptic integral of the first kind F(φ, m).
  789.      * <p>
  790.      * <em>
  791.      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
  792.      * are considered experimental for now, they have known issues.
  793.      * </em>
  794.      * </p>
  795.      * <p>
  796.      * The incomplete elliptic integral of the first kind F(φ, m) is
  797.      * \[
  798.      *    \int_0^{\phi} \frac{d\theta}{\sqrt{1-m \sin^2\theta}}
  799.      * \]
  800.      * </p>
  801.      * <p>
  802.      * The algorithm for evaluating the functions is based on numerical integration.
  803.      * If integration path comes too close to a pole of the integrand, then integration will fail
  804.      * with a {@link org.hipparchus.exception.MathIllegalStateException MathIllegalStateException}
  805.      * even for very large {@code maxEval}. This is normal behavior.
  806.      * </p>
  807.      * @param phi amplitude (i.e. upper bound of the integral)
  808.      * @param m parameter (m=k² where k is the elliptic modulus)
  809.      * @param integrator integrator to use
  810.      * @param maxEval maximum number of evaluations (real and imaginary
  811.      * parts are evaluated separately, so up to twice this number may be used)
  812.      * @param <T> the type of the field elements
  813.      * @return incomplete elliptic integral of the first kind F(φ, m)
  814.      * @see #bigK(CalculusFieldElement)
  815.      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheFirstKind.html">Elliptic Integrals of the First Kind (MathWorld)</a>
  816.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  817.      */
  818.     public static <T extends CalculusFieldElement<T>> FieldComplex<T> bigF(final FieldComplex<T> phi, final FieldComplex<T> m,
  819.                                                                            final FieldComplexUnivariateIntegrator<T> integrator,
  820.                                                                            final int maxEval) {
  821.         return integrator.integrate(maxEval, new First<>(m), phi.getField().getZero(), phi);
  822.     }

  824.     /** Get the incomplete elliptic integral of the second kind E(φ, m).
  825.      * <p>
  826.      * The incomplete elliptic integral of the second kind E(φ, m) is
  827.      * \[
  828.      *    \int_0^{\phi} \sqrt{1-m \sin^2\theta} d\theta
  829.      * \]
  830.      * </p>
  831.      * <p>
  832.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  833.      * Carlson elliptic integrals}.
  834.      * </p>
  835.      * @param phi amplitude (i.e. upper bound of the integral)
  836.      * @param m parameter (m=k² where k is the elliptic modulus)
  837.      * @return incomplete elliptic integral of the second kind E(φ, m)
  838.      * @see #bigE(double)
  839.      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheSecondKind.html">Elliptic Integrals of the Second Kind (MathWorld)</a>
  840.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  841.      */
  842.     public static double bigE(final double phi, final double m) {

  844.         // argument reduction
  845.         final DoubleArgumentReduction ar = new DoubleArgumentReduction(phi, m, LegendreEllipticIntegral::bigE);

  847.         // integrate part between 0 and π/2
  848.         final double cM1        = ar.csc2 - 1.0;
  849.         final double cMm        = ar.csc2 - m;
  850.         final double incomplete = CarlsonEllipticIntegral.rF(cM1, cMm, ar.csc2) -
  851.                                   CarlsonEllipticIntegral.rD(cM1, cMm, ar.csc2) * (m / 3);

  853.         // combine complete and incomplete parts
  854.         return ar.negate ? ar.complete - incomplete : ar.complete + incomplete;

  856.     }

  858.     /** Get the incomplete elliptic integral of the second kind E(φ, m).
  859.      * <p>
  860.      * The incomplete elliptic integral of the second kind E(φ, m) is
  861.      * \[
  862.      *    \int_0^{\phi} \sqrt{1-m \sin^2\theta} d\theta
  863.      * \]
  864.      * </p>
  865.      * <p>
  866.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  867.      * Carlson elliptic integrals}.
  868.      * </p>
  869.      * @param phi amplitude (i.e. upper bound of the integral)
  870.      * @param m parameter (m=k² where k is the elliptic modulus)
  871.      * @param <T> the type of the field elements
  872.      * @return incomplete elliptic integral of the second kind E(φ, m)
  873.      * @see #bigE(CalculusFieldElement)
  874.      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheSecondKind.html">Elliptic Integrals of the Second Kind (MathWorld)</a>
  875.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  876.      */
  877.     public static <T extends CalculusFieldElement<T>> T bigE(final T phi, final T m) {

  879.         // argument reduction
  880.         final FieldArgumentReduction<T> ar = new FieldArgumentReduction<>(phi, m, LegendreEllipticIntegral::bigE);

  882.         // integrate part between 0 and π/2
  883.         final T cM1        = ar.csc2.subtract(1);
  884.         final T cMm        = ar.csc2.subtract(m);
  885.         final T incomplete = CarlsonEllipticIntegral.rF(cM1, cMm, ar.csc2).
  886.                              subtract(CarlsonEllipticIntegral.rD(cM1, cMm, ar.csc2).multiply(m.divide(3)));

  888.         // combine complete and incomplete parts
  889.         return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);

  891.     }

  893.     /** Get the incomplete elliptic integral of the second kind E(φ, m).
  894.      * <p>
  895.      * <em>
  896.      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
  897.      * are considered experimental for now, they have known issues.
  898.      * </em>
  899.      * </p>
  900.      * <p>
  901.      * The incomplete elliptic integral of the second kind E(φ, m) is
  902.      * \[
  903.      *    \int_0^{\phi} \sqrt{1-m \sin^2\theta} d\theta
  904.      * \]
  905.      * </p>
  906.      * <p>
  907.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  908.      * Carlson elliptic integrals}.
  909.      * </p>
  910.      * @param phi amplitude (i.e. upper bound of the integral)
  911.      * @param m parameter (m=k² where k is the elliptic modulus)
  912.      * @return incomplete elliptic integral of the second kind E(φ, m)
  913.      * @see #bigE(Complex)
  914.      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheSecondKind.html">Elliptic Integrals of the Second Kind (MathWorld)</a>
  915.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  916.      */
  917.     public static Complex bigE(final Complex phi, final Complex m) {

  919.         // argument reduction
  920.         final FieldArgumentReduction<Complex> ar = new FieldArgumentReduction<>(phi, m, LegendreEllipticIntegral::bigE);

  922.         // integrate part between 0 and π/2
  923.         final Complex cM1        = ar.csc2.subtract(1);
  924.         final Complex cMm        = ar.csc2.subtract(m);
  925.         final Complex incomplete = CarlsonEllipticIntegral.rF(cM1, cMm, ar.csc2).
  926.                                    subtract(CarlsonEllipticIntegral.rD(cM1, cMm, ar.csc2).multiply(m.divide(3)));

  928.         // combine complete and incomplete parts
  929.         return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);

  931.     }

  933.     /** Get the incomplete elliptic integral of the second kind E(φ, m) using numerical integration.
  934.      * <p>
  935.      * <em>
  936.      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
  937.      * are considered experimental for now, they have known issues.
  938.      * </em>
  939.      * </p>
  940.      * <p>
  941.      * The incomplete elliptic integral of the second kind E(φ, m) is
  942.      * \[
  943.      *    \int_0^{\phi} \sqrt{1-m \sin^2\theta} d\theta
  944.      * \]
  945.      * </p>
  946.      * <p>
  947.      * The algorithm for evaluating the functions is based on numerical integration.
  948.      * If integration path comes too close to a pole of the integrand, then integration will fail
  949.      * with a {@link org.hipparchus.exception.MathIllegalStateException MathIllegalStateException}
  950.      * even for very large {@code maxEval}. This is normal behavior.
  951.      * </p>
  952.      * @param phi amplitude (i.e. upper bound of the integral)
  953.      * @param m parameter (m=k² where k is the elliptic modulus)
  954.      * @param integrator integrator to use
  955.      * @param maxEval maximum number of evaluations (real and imaginary
  956.      * parts are evaluated separately, so up to twice this number may be used)
  957.      * @return incomplete elliptic integral of the second kind E(φ, m)
  958.      * @see #bigE(Complex)
  959.      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheSecondKind.html">Elliptic Integrals of the Second Kind (MathWorld)</a>
  960.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  961.      */
  962.     public static Complex bigE(final Complex phi, final Complex m,
  963.                                final ComplexUnivariateIntegrator integrator, final int maxEval) {
  964.         return integrator.integrate(maxEval, new Second<>(m), phi.getField().getZero(), phi);
  965.     }

  967.     /** Get the incomplete elliptic integral of the second kind E(φ, m).
  968.      * <p>
  969.      * <em>
  970.      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
  971.      * are considered experimental for now, they have known issues.
  972.      * </em>
  973.      * </p>
  974.      * <p>
  975.      * The incomplete elliptic integral of the second kind E(φ, m) is
  976.      * \[
  977.      *    \int_0^{\phi} \sqrt{1-m \sin^2\theta} d\theta
  978.      * \]
  979.      * </p>
  980.      * <p>
  981.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  982.      * Carlson elliptic integrals}.
  983.      * </p>
  984.      * @param phi amplitude (i.e. upper bound of the integral)
  985.      * @param m parameter (m=k² where k is the elliptic modulus)
  986.      * @param <T> the type of the field elements
  987.      * @return incomplete elliptic integral of the second kind E(φ, m)
  988.      * @see #bigE(FieldComplex)
  989.      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheSecondKind.html">Elliptic Integrals of the Second Kind (MathWorld)</a>
  990.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  991.      */
  992.     public static <T extends CalculusFieldElement<T>> FieldComplex<T> bigE(final FieldComplex<T> phi, final FieldComplex<T> m) {

  994.         // argument reduction
  995.         final FieldArgumentReduction<FieldComplex<T>> ar = new FieldArgumentReduction<>(phi, m, LegendreEllipticIntegral::bigE);

  997.         // integrate part between 0 and π/2
  998.         final FieldComplex<T> cM1        = ar.csc2.subtract(1);
  999.         final FieldComplex<T> cMm        = ar.csc2.subtract(m);
  1000.         final FieldComplex<T> incomplete = CarlsonEllipticIntegral.rF(cM1, cMm, ar.csc2).
  1001.                                            subtract(CarlsonEllipticIntegral.rD(cM1, cMm, ar.csc2).multiply(m.divide(3)));

  1003.         // combine complete and incomplete parts
  1004.         return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);

  1006.     }

  1008.     /** Get the incomplete elliptic integral of the second kind E(φ, m).
  1009.      * <p>
  1010.      * <em>
  1011.      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
  1012.      * are considered experimental for now, they have known issues.
  1013.      * </em>
  1014.      * </p>
  1015.      * <p>
  1016.      * The incomplete elliptic integral of the second kind E(φ, m) is
  1017.      * \[
  1018.      *    \int_0^{\phi} \sqrt{1-m \sin^2\theta} d\theta
  1019.      * \]
  1020.      * </p>
  1021.      * <p>
  1022.      * The algorithm for evaluating the functions is based on numerical integration.
  1023.      * If integration path comes too close to a pole of the integrand, then integration will fail
  1024.      * with a {@link org.hipparchus.exception.MathIllegalStateException MathIllegalStateException}
  1025.      * even for very large {@code maxEval}. This is normal behavior.
  1026.      * </p>
  1027.      * @param phi amplitude (i.e. upper bound of the integral)
  1028.      * @param m parameter (m=k² where k is the elliptic modulus)
  1029.      * @param integrator integrator to use
  1030.      * @param maxEval maximum number of evaluations (real and imaginary
  1031.      * parts are evaluated separately, so up to twice this number may be used)
  1032.      * @param <T> the type of the field elements
  1033.      * @return incomplete elliptic integral of the second kind E(φ, m)
  1034.      * @see #bigE(FieldComplex)
  1035.      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheSecondKind.html">Elliptic Integrals of the Second Kind (MathWorld)</a>
  1036.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  1037.      */
  1038.     public static <T extends CalculusFieldElement<T>> FieldComplex<T> bigE(final FieldComplex<T> phi, final FieldComplex<T> m,
  1039.                                                                            final FieldComplexUnivariateIntegrator<T> integrator,
  1040.                                                                            final int maxEval) {
  1041.         return integrator.integrate(maxEval, new Second<>(m), phi.getField().getZero(), phi);
  1042.     }

  1044.     /** Get the incomplete elliptic integral D(φ, m) = [F(φ, m) - E(φ, m)]/m.
  1045.      * <p>
  1046.      * The incomplete elliptic integral D(φ, m) is
  1047.      * \[
  1048.      *    \int_0^{\phi} \frac{\sin^2\theta}{\sqrt{1-m \sin^2\theta}} d\theta
  1049.      * \]
  1050.      * </p>
  1051.      * <p>
  1052.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  1053.      * Carlson elliptic integrals}.
  1054.      * </p>
  1055.      * @param phi amplitude (i.e. upper bound of the integral)
  1056.      * @param m parameter (m=k² where k is the elliptic modulus)
  1057.      * @return incomplete elliptic integral D(φ, m)
  1058.      * @see #bigD(double)
  1059.      */
  1060.     public static double bigD(final double phi, final double m) {

  1062.         // argument reduction
  1063.         final DoubleArgumentReduction ar = new DoubleArgumentReduction(phi, m, LegendreEllipticIntegral::bigD);

  1065.         // integrate part between 0 and π/2
  1066.         final double cM1        = ar.csc2 - 1.0;
  1067.         final double cMm        = ar.csc2 - m;
  1068.         final double incomplete = CarlsonEllipticIntegral.rD(cM1, cMm, ar.csc2) / 3;

  1070.         // combine complete and incomplete parts
  1071.         return ar.negate ? ar.complete - incomplete : ar.complete + incomplete;

  1073.     }

  1075.     /** Get the incomplete elliptic integral D(φ, m) = [F(φ, m) - E(φ, m)]/m.
  1076.      * <p>
  1077.      * The incomplete elliptic integral D(φ, m) is
  1078.      * \[
  1079.      *    \int_0^{\phi} \frac{\sin^2\theta}{\sqrt{1-m \sin^2\theta}} d\theta
  1080.      * \]
  1081.      * </p>
  1082.      * <p>
  1083.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  1084.      * Carlson elliptic integrals}.
  1085.      * </p>
  1086.      * @param phi amplitude (i.e. upper bound of the integral)
  1087.      * @param m parameter (m=k² where k is the elliptic modulus)
  1088.      * @param <T> the type of the field elements
  1089.      * @return incomplete elliptic integral D(φ, m)
  1090.      * @see #bigD(CalculusFieldElement)
  1091.      */
  1092.     public static <T extends CalculusFieldElement<T>> T bigD(final T phi, final T m) {

  1094.         // argument reduction
  1095.         final FieldArgumentReduction<T> ar = new FieldArgumentReduction<>(phi, m, LegendreEllipticIntegral::bigD);

  1097.         // integrate part between 0 and π/2
  1098.         final T cM1        = ar.csc2.subtract(1);
  1099.         final T cMm        = ar.csc2.subtract(m);
  1100.         final T incomplete = CarlsonEllipticIntegral.rD(cM1, cMm, ar.csc2).divide(3);

  1102.         // combine complete and incomplete parts
  1103.         return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);

  1105.     }

  1107.     /** Get the incomplete elliptic integral D(φ, m) = [F(φ, m) - E(φ, m)]/m.
  1108.      * <p>
  1109.      * <em>
  1110.      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
  1111.      * are considered experimental for now, they have known issues.
  1112.      * </em>
  1113.      * </p>
  1114.      * <p>
  1115.      * The incomplete elliptic integral D(φ, m) is
  1116.      * \[
  1117.      *    \int_0^{\phi} \frac{\sin^2\theta}{\sqrt{1-m \sin^2\theta}} d\theta
  1118.      * \]
  1119.      * </p>
  1120.      * <p>
  1121.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  1122.      * Carlson elliptic integrals}.
  1123.      * </p>
  1124.      * @param phi amplitude (i.e. upper bound of the integral)
  1125.      * @param m parameter (m=k² where k is the elliptic modulus)
  1126.      * @return incomplete elliptic integral D(φ, m)
  1127.      * @see #bigD(Complex)
  1128.      */
  1129.     public static Complex bigD(final Complex phi, final Complex m) {

  1131.         // argument reduction
  1132.         final FieldArgumentReduction<Complex> ar = new FieldArgumentReduction<>(phi, m, LegendreEllipticIntegral::bigD);

  1134.         // integrate part between 0 and π/2
  1135.         final Complex cM1        = ar.csc2.subtract(1);
  1136.         final Complex cMm        = ar.csc2.subtract(m);
  1137.         final Complex incomplete = CarlsonEllipticIntegral.rD(cM1, cMm, ar.csc2).divide(3);

  1139.         // combine complete and incomplete parts
  1140.         return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);

  1142.     }

  1144.     /** Get the incomplete elliptic integral D(φ, m) = [F(φ, m) - E(φ, m)]/m.
  1145.      * <p>
  1146.      * <em>
  1147.      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
  1148.      * are considered experimental for now, they have known issues.
  1149.      * </em>
  1150.      * </p>
  1151.      * <p>
  1152.      * The incomplete elliptic integral D(φ, m) is
  1153.      * \[
  1154.      *    \int_0^{\phi} \frac{\sin^2\theta}{\sqrt{1-m \sin^2\theta}} d\theta
  1155.      * \]
  1156.      * </p>
  1157.      * <p>
  1158.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  1159.      * Carlson elliptic integrals}.
  1160.      * </p>
  1161.      * @param phi amplitude (i.e. upper bound of the integral)
  1162.      * @param m parameter (m=k² where k is the elliptic modulus)
  1163.      * @param <T> the type of the field elements
  1164.      * @return incomplete elliptic integral D(φ, m)
  1165.      * @see #bigD(CalculusFieldElement)
  1166.      */
  1167.     public static <T extends CalculusFieldElement<T>> FieldComplex<T> bigD(final FieldComplex<T> phi, final FieldComplex<T> m) {

  1169.         // argument reduction
  1170.         final FieldArgumentReduction<FieldComplex<T>> ar = new FieldArgumentReduction<>(phi, m, LegendreEllipticIntegral::bigD);

  1172.         // integrate part between 0 and π/2
  1173.         final FieldComplex<T> cM1        = ar.csc2.subtract(1);
  1174.         final FieldComplex<T> cMm        = ar.csc2.subtract(m);
  1175.         final FieldComplex<T> incomplete = CarlsonEllipticIntegral.rD(cM1, cMm, ar.csc2).divide(3);

  1177.         // combine complete and incomplete parts
  1178.         return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);

  1180.     }

  1182.     /** Get the incomplete elliptic integral of the third kind Π(n, φ, m).
  1183.      * <p>
  1184.      * The incomplete elliptic integral of the third kind Π(n, φ, m) is
  1185.      * \[
  1186.      *    \int_0^{\phi} \frac{d\theta}{\sqrt{1-m \sin^2\theta}(1-n \sin^2\theta)}
  1187.      * \]
  1188.      * </p>
  1189.      * <p>
  1190.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  1191.      * Carlson elliptic integrals}.
  1192.      * </p>
  1193.      * @param n elliptic characteristic
  1194.      * @param phi amplitude (i.e. upper bound of the integral)
  1195.      * @param m parameter (m=k² where k is the elliptic modulus)
  1196.      * @return incomplete elliptic integral of the third kind Π(n, φ, m)
  1197.      * @see #bigPi(double, double)
  1198.      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheThirdKind.html">Elliptic Integrals of the Third Kind (MathWorld)</a>
  1199.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  1200.      */
  1201.     public static double bigPi(final double n, final double phi, final double m) {

  1203.         // argument reduction
  1204.         final DoubleArgumentReduction ar = new DoubleArgumentReduction(phi, m, parameter -> bigPi(n, parameter));

  1206.         // integrate part between 0 and π/2
  1207.         final double cM1        = ar.csc2 - 1.0;
  1208.         final double cMm        = ar.csc2 - m;
  1209.         final double cMn        = ar.csc2 - n;
  1210.         final double delta      = n * (m - n) * (n - 1);
  1211.         final double incomplete = CarlsonEllipticIntegral.rF(cM1, cMm, ar.csc2) +
  1212.                                   CarlsonEllipticIntegral.rJ(cM1, cMm, ar.csc2, cMn, delta) * n / 3;

  1214.         // combine complete and incomplete parts
  1215.         return ar.negate ? ar.complete - incomplete : ar.complete + incomplete;

  1217.     }

  1219.     /** Get the incomplete elliptic integral of the third kind Π(n, φ, m).
  1220.      * <p>
  1221.      * The incomplete elliptic integral of the third kind Π(n, φ, m) is
  1222.      * \[
  1223.      *    \int_0^{\phi} \frac{d\theta}{\sqrt{1-m \sin^2\theta}(1-n \sin^2\theta)}
  1224.      * \]
  1225.      * </p>
  1226.      * <p>
  1227.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  1228.      * Carlson elliptic integrals}.
  1229.      * </p>
  1230.      * @param n elliptic characteristic
  1231.      * @param phi amplitude (i.e. upper bound of the integral)
  1232.      * @param m parameter (m=k² where k is the elliptic modulus)
  1233.      * @param <T> the type of the field elements
  1234.      * @return incomplete elliptic integral of the third kind Π(n, φ, m)
  1235.      * @see #bigPi(CalculusFieldElement, CalculusFieldElement)
  1236.      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheThirdKind.html">Elliptic Integrals of the Third Kind (MathWorld)</a>
  1237.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  1238.      */
  1239.     public static <T extends CalculusFieldElement<T>> T bigPi(final T n, final T phi, final T m) {

  1241.         // argument reduction
  1242.         final FieldArgumentReduction<T> ar = new FieldArgumentReduction<>(phi, m, parameter -> bigPi(n, parameter));

  1244.         // integrate part between 0 and π/2
  1245.         final T cM1        = ar.csc2.subtract(1);
  1246.         final T cMm        = ar.csc2.subtract(m);
  1247.         final T cMn        = ar.csc2.subtract(n);
  1248.         final T delta      = n.multiply(m.subtract(n)).multiply(n.subtract(1));
  1249.         final T incomplete = CarlsonEllipticIntegral.rF(cM1, cMm, ar.csc2).
  1250.                              add(CarlsonEllipticIntegral.rJ(cM1, cMm, ar.csc2, cMn, delta).multiply(n).divide(3));

  1252.         // combine complete and incomplete parts
  1253.         return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);

  1255.     }

  1257.     /** Get the incomplete elliptic integral of the third kind Π(n, φ, m).
  1258.      * <p>
  1259.      * <em>
  1260.      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
  1261.      * are considered experimental for now, they have known issues.
  1262.      * </em>
  1263.      * </p>
  1264.      * <p>
  1265.      * The incomplete elliptic integral of the third kind Π(n, φ, m) is
  1266.      * \[
  1267.      *    \int_0^{\phi} \frac{d\theta}{\sqrt{1-m \sin^2\theta}(1-n \sin^2\theta)}
  1268.      * \]
  1269.      * </p>
  1270.      * <p>
  1271.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  1272.      * Carlson elliptic integrals}.
  1273.      * </p>
  1274.      * @param n elliptic characteristic
  1275.      * @param phi amplitude (i.e. upper bound of the integral)
  1276.      * @param m parameter (m=k² where k is the elliptic modulus)
  1277.      * @return incomplete elliptic integral of the third kind Π(n, φ, m)
  1278.      * @see #bigPi(Complex, Complex)
  1279.      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheThirdKind.html">Elliptic Integrals of the Third Kind (MathWorld)</a>
  1280.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  1281.      */
  1282.     public static Complex bigPi(final Complex n, final Complex phi, final Complex m) {

  1284.         // argument reduction
  1285.         final FieldArgumentReduction<Complex> ar = new FieldArgumentReduction<>(phi, m, parameter -> bigPi(n, parameter));

  1287.         // integrate part between 0 and π/2
  1288.         final Complex cM1        = ar.csc2.subtract(1);
  1289.         final Complex cMm        = ar.csc2.subtract(m);
  1290.         final Complex cMn        = ar.csc2.subtract(n);
  1291.         final Complex delta      = n.multiply(m.subtract(n)).multiply(n.subtract(1));
  1292.         final Complex incomplete = CarlsonEllipticIntegral.rF(cM1, cMm, ar.csc2).
  1293.                                    add(CarlsonEllipticIntegral.rJ(cM1, cMm, ar.csc2, cMn, delta).multiply(n).divide(3));

  1295.         // combine complete and incomplete parts
  1296.         return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);

  1298.     }

  1300.     /** Get the incomplete elliptic integral of the third kind Π(n, φ, m) using numerical integration.
  1301.      * <p>
  1302.      * <em>
  1303.      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
  1304.      * are considered experimental for now, they have known issues.
  1305.      * </em>
  1306.      * </p>
  1307.      * <p>
  1308.      * The incomplete elliptic integral of the third kind Π(n, φ, m) is
  1309.      * \[
  1310.      *    \int_0^{\phi} \frac{d\theta}{\sqrt{1-m \sin^2\theta}(1-n \sin^2\theta)}
  1311.      * \]
  1312.      * </p>
  1313.      * <p>
  1314.      * The algorithm for evaluating the functions is based on numerical integration.
  1315.      * If integration path comes too close to a pole of the integrand, then integration will fail
  1316.      * with a {@link org.hipparchus.exception.MathIllegalStateException MathIllegalStateException}
  1317.      * even for very large {@code maxEval}. This is normal behavior.
  1318.      * </p>
  1319.      * @param n elliptic characteristic
  1320.      * @param phi amplitude (i.e. upper bound of the integral)
  1321.      * @param m parameter (m=k² where k is the elliptic modulus)
  1322.      * @param integrator integrator to use
  1323.      * @param maxEval maximum number of evaluations (real and imaginary
  1324.      * @return incomplete elliptic integral of the third kind Π(n, φ, m)
  1325.      * @see #bigPi(Complex, Complex)
  1326.      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheThirdKind.html">Elliptic Integrals of the Third Kind (MathWorld)</a>
  1327.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  1328.      */
  1329.     public static Complex bigPi(final Complex n, final Complex phi, final Complex m,
  1330.                                 final ComplexUnivariateIntegrator integrator, final int maxEval) {
  1331.          return integrator.integrate(maxEval, new Third<>(n, m), phi.getField().getZero(), phi);
  1332.     }

  1334.     /** Get the incomplete elliptic integral of the third kind Π(n, φ, m).
  1335.      * <p>
  1336.      * <em>
  1337.      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
  1338.      * are considered experimental for now, they have known issues.
  1339.      * </em>
  1340.      * </p>
  1341.      * <p>
  1342.      * The incomplete elliptic integral of the third kind Π(n, φ, m) is
  1343.      * \[
  1344.      *    \int_0^{\phi} \frac{d\theta}{\sqrt{1-m \sin^2\theta}(1-n \sin^2\theta)}
  1345.      * \]
  1346.      * </p>
  1347.      * <p>
  1348.      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
  1349.      * Carlson elliptic integrals}.
  1350.      * </p>
  1351.      * @param n elliptic characteristic
  1352.      * @param phi amplitude (i.e. upper bound of the integral)
  1353.      * @param m parameter (m=k² where k is the elliptic modulus)
  1354.      * @param <T> the type of the field elements
  1355.      * @return incomplete elliptic integral of the third kind Π(n, φ, m)
  1356.      * @see #bigPi(FieldComplex, FieldComplex)
  1357.      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheThirdKind.html">Elliptic Integrals of the Third Kind (MathWorld)</a>
  1358.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  1359.      */
  1360.     public static <T extends CalculusFieldElement<T>> FieldComplex<T> bigPi(final FieldComplex<T> n,
  1361.                                                                             final FieldComplex<T> phi,
  1362.                                                                             final FieldComplex<T> m) {

  1364.         // argument reduction
  1365.         final FieldArgumentReduction<FieldComplex<T>> ar = new FieldArgumentReduction<>(phi, m, parameter -> bigPi(n, parameter));

  1367.         // integrate part between 0 and π/2
  1368.         final FieldComplex<T> cM1        = ar.csc2.subtract(1);
  1369.         final FieldComplex<T> cMm        = ar.csc2.subtract(m);
  1370.         final FieldComplex<T> cMn        = ar.csc2.subtract(n);
  1371.         final FieldComplex<T> delta      = n.multiply(m.subtract(n)).multiply(n.subtract(1));
  1372.         final FieldComplex<T> incomplete = CarlsonEllipticIntegral.rF(cM1, cMm, ar.csc2).
  1373.                                            add(CarlsonEllipticIntegral.rJ(cM1, cMm, ar.csc2, cMn, delta).multiply(n).divide(3));

  1375.         // combine complete and incomplete parts
  1376.         return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);

  1378.     }

  1380.     /** Get the incomplete elliptic integral of the third kind Π(n, φ, m).
  1381.      * <p>
  1382.      * <em>
  1383.      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
  1384.      * are considered experimental for now, they have known issues.
  1385.      * </em>
  1386.      * </p>
  1387.      * <p>
  1388.      * The incomplete elliptic integral of the third kind Π(n, φ, m) is
  1389.      * \[
  1390.      *    \int_0^{\phi} \frac{d\theta}{\sqrt{1-m \sin^2\theta}(1-n \sin^2\theta)}
  1391.      * \]
  1392.      * </p>
  1393.      * <p>
  1394.      * The algorithm for evaluating the functions is based on numerical integration.
  1395.      * If integration path comes too close to a pole of the integrand, then integration will fail
  1396.      * with a {@link org.hipparchus.exception.MathIllegalStateException MathIllegalStateException}
  1397.      * even for very large {@code maxEval}. This is normal behavior.
  1398.      * </p>
  1399.      * @param n elliptic characteristic
  1400.      * @param phi amplitude (i.e. upper bound of the integral)
  1401.      * @param m parameter (m=k² where k is the elliptic modulus)
  1402.      * @param integrator integrator to use
  1403.      * @param maxEval maximum number of evaluations (real and imaginary
  1404.      * parts are evaluated separately, so up to twice this number may be used)
  1405.      * @param <T> the type of the field elements
  1406.      * @return incomplete elliptic integral of the third kind Π(n, φ, m)
  1407.      * @see #bigPi(FieldComplex, FieldComplex)
  1408.      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheThirdKind.html">Elliptic Integrals of the Third Kind (MathWorld)</a>
  1409.      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
  1410.      */
  1411.     public static <T extends CalculusFieldElement<T>> FieldComplex<T> bigPi(final FieldComplex<T> n,
  1412.                                                                             final FieldComplex<T> phi,
  1413.                                                                             final FieldComplex<T> m,
  1414.                                                                             final FieldComplexUnivariateIntegrator<T> integrator,
  1415.                                                                             final int maxEval) {
  1416.         return integrator.integrate(maxEval, new Third<>(n, m), phi.getField().getZero(), phi);
  1417.     }

  1419.     /** Argument reduction for an incomplete integral. */
  1420.     private static class DoubleArgumentReduction {

  1422.         /** Complete part. */
  1423.         private final double complete;

  1425.         /** Squared cosecant of the Jacobi amplitude. */
  1426.         private final double csc2;

  1428.         /** Indicator for negated Jacobi amplitude. */
  1429.         private boolean negate;

  1431.         /** Simple constructor.
  1432.          * @param phi amplitude (i.e. upper bound of the integral)
  1433.          * @param m parameter (m=k² where k is the elliptic modulus)
  1434.          * @param integral provider for complete integral
  1435.          */
  1436.         DoubleArgumentReduction(final double phi, final double m, final DoubleFunction<Double> integral) {
  1437.             final double sin = FastMath.sin(phi);
  1438.             final int    p   = (int) FastMath.rint(phi / FastMath.PI);
  1439.             complete         = p == 0 ? 0 : integral.apply(m) * 2 * p;
  1440.             negate           = sin < 0 ^ (p & 0x1) == 1;
  1441.             csc2             = 1.0 / (sin * sin);
  1442.         }

  1444.     }

  1446.     /** Argument reduction for an incomplete integral.
  1447.      * @param <T> type fo the field elements
  1448.      */
  1449.     private static class FieldArgumentReduction<T extends CalculusFieldElement<T>> {

  1451.         /** Complete part. */
  1452.         private final T complete;

  1454.         /** Squared cosecant of the Jacobi amplitude. */
  1455.         private final T csc2;

  1457.         /** Indicator for negated Jacobi amplitude. */
  1458.         private boolean negate;

  1460.         /** Simple constructor.
  1461.          * @param phi amplitude (i.e. upper bound of the integral)
  1462.          * @param m parameter (m=k² where k is the elliptic modulus)
  1463.          * @param integral provider for complete integral
  1464.          */
  1465.         FieldArgumentReduction(final T phi, final T m, final Function<T, T> integral) {
  1466.             final T   sin = FastMath.sin(phi);
  1467.             final int p   = (int) FastMath.rint(phi.getReal() / FastMath.PI);
  1468.             complete      = p == 0 ? phi.getField().getZero() : integral.apply(m).multiply(2 * p);
  1469.             negate        = sin.getReal() < 0 ^ (p & 0x1) == 1;
  1470.             csc2          = sin.multiply(sin).reciprocal();
  1471.         }

  1473.     }

  1475.     /** Integrand for elliptic integrals of the first kind.
  1476.      * @param <T> type of the field elements
  1477.      */
  1478.     private static class First<T extends CalculusFieldElement<T>> implements CalculusFieldUnivariateFunction<T> {

  1480.         /** Parameter. */
  1481.         private final T m;

  1483.         /** Simple constructor.
  1484.          * @param m parameter (m=k² where k is the elliptic modulus)
  1485.          */
  1486.         First(final T m) {
  1487.             this.m = m;
  1488.         }

  1490.         /** {@inheritDoc} */
  1491.         @Override
  1492.         public T value(final T theta) {
  1493.             final T sin  = theta.sin();
  1494.             final T sin2 = sin.multiply(sin);
  1495.             return sin2.multiply(m).negate().add(1).sqrt().reciprocal();
  1496.         }

  1498.     }

  1500.     /** Integrand for elliptic integrals of the second kind.
  1501.      * @param <T> type of the field elements
  1502.      */
  1503.     private static class Second<T extends CalculusFieldElement<T>> implements CalculusFieldUnivariateFunction<T> {

  1505.         /** Parameter. */
  1506.         private final T m;

  1508.         /** Simple constructor.
  1509.          * @param m parameter (m=k² where k is the elliptic modulus)
  1510.          */
  1511.         Second(final T m) {
  1512.             this.m = m;
  1513.         }

  1515.         /** {@inheritDoc} */
  1516.         @Override
  1517.         public T value(final T theta) {
  1518.             final T sin = theta.sin();
  1519.             final T sin2 = sin.multiply(sin);
  1520.             return sin2.multiply(m).negate().add(1).sqrt();
  1521.         }

  1523.     }

  1525.     /** Integrand for elliptic integrals of the third kind.
  1526.      * @param <T> type of the field elements
  1527.      */
  1528.     private static class Third<T extends CalculusFieldElement<T>> implements CalculusFieldUnivariateFunction<T> {

  1530.         /** Elliptic characteristic. */
  1531.         private final T n;

  1533.         /** Parameter. */
  1534.         private final T m;

  1536.         /** Simple constructor.
  1537.          * @param n elliptic characteristic
  1538.          * @param m parameter (m=k² where k is the elliptic modulus)
  1539.          */
  1540.         Third(final T n, final T m) {
  1541.             this.n = n;
  1542.             this.m = m;
  1543.         }

  1545.         /** {@inheritDoc} */
  1546.         @Override
  1547.         public T value(final T theta) {
  1548.             final T sin  = theta.sin();
  1549.             final T sin2 = sin.multiply(sin);
  1550.             final T d1   = sin2.multiply(m).negate().add(1).sqrt();
  1551.             final T da   = sin2.multiply(n).negate().add(1);
  1552.             return d1.multiply(da).reciprocal();
  1553.         }

  1555.     }
  1556. }
Created with JaCoCo 0.8.12.202403310830