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;
18
19 import java.util.function.DoubleFunction;
20 import java.util.function.Function;
21
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;
31
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
75
76 /** Private constructor for a utility class.
77 */
78 private LegendreEllipticIntegral() {
79 // nothing to do
80 }
81
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 }
95
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 }
111
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 }
139
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()) {
164
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));
167
168 } else {
169 return CarlsonEllipticIntegral.rF(zero, one.subtract(m), one);
170 }
171 }
172
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 }
200
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()) {
225
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));
228
229 } else {
230 return CarlsonEllipticIntegral.rF(zero, one.subtract(m), one);
231 }
232 }
233
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 }
255
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 }
280
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 }
302
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 }
327
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 }
348
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 }
372
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 }
395
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 }
419
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 }
438
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 }
460
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 }
479
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 }
501
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 }
526
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 }
554
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 }
579
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 }
607
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) {
627
628 // argument reduction
629 final DoubleArgumentReduction ar = new DoubleArgumentReduction(phi, m, LegendreEllipticIntegral::bigK);
630
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);
635
636 // combine complete and incomplete parts
637 return ar.negate ? ar.complete - incomplete : ar.complete + incomplete;
638
639 }
640
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) {
661
662 // argument reduction
663 final FieldArgumentReduction<T> ar = new FieldArgumentReduction<>(phi, m, LegendreEllipticIntegral::bigK);
664
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);
669
670 // combine complete and incomplete parts
671 return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);
672
673 }
674
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) {
700
701 // argument reduction
702 final FieldArgumentReduction<Complex> ar = new FieldArgumentReduction<>(phi, m, LegendreEllipticIntegral::bigK);
703
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);
708
709 // combine complete and incomplete parts
710 return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);
711
712 }
713
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 }
747
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) {
774
775 // argument reduction
776 final FieldArgumentReduction<FieldComplex<T>> ar = new FieldArgumentReduction<>(phi, m, LegendreEllipticIntegral::bigK);
777
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);
782
783 // combine complete and incomplete parts
784 return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);
785
786 }
787
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 }
823
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) {
843
844 // argument reduction
845 final DoubleArgumentReduction ar = new DoubleArgumentReduction(phi, m, LegendreEllipticIntegral::bigE);
846
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);
852
853 // combine complete and incomplete parts
854 return ar.negate ? ar.complete - incomplete : ar.complete + incomplete;
855
856 }
857
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) {
878
879 // argument reduction
880 final FieldArgumentReduction<T> ar = new FieldArgumentReduction<>(phi, m, LegendreEllipticIntegral::bigE);
881
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)));
887
888 // combine complete and incomplete parts
889 return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);
890
891 }
892
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) {
918
919 // argument reduction
920 final FieldArgumentReduction<Complex> ar = new FieldArgumentReduction<>(phi, m, LegendreEllipticIntegral::bigE);
921
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)));
927
928 // combine complete and incomplete parts
929 return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);
930
931 }
932
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 }
966
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) {
993
994 // argument reduction
995 final FieldArgumentReduction<FieldComplex<T>> ar = new FieldArgumentReduction<>(phi, m, LegendreEllipticIntegral::bigE);
996
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)));
1002
1003 // combine complete and incomplete parts
1004 return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);
1005
1006 }
1007
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 }
1043
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) {
1061
1062 // argument reduction
1063 final DoubleArgumentReduction ar = new DoubleArgumentReduction(phi, m, LegendreEllipticIntegral::bigD);
1064
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;
1069
1070 // combine complete and incomplete parts
1071 return ar.negate ? ar.complete - incomplete : ar.complete + incomplete;
1072
1073 }
1074
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) {
1093
1094 // argument reduction
1095 final FieldArgumentReduction<T> ar = new FieldArgumentReduction<>(phi, m, LegendreEllipticIntegral::bigD);
1096
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);
1101
1102 // combine complete and incomplete parts
1103 return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);
1104
1105 }
1106
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) {
1130
1131 // argument reduction
1132 final FieldArgumentReduction<Complex> ar = new FieldArgumentReduction<>(phi, m, LegendreEllipticIntegral::bigD);
1133
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);
1138
1139 // combine complete and incomplete parts
1140 return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);
1141
1142 }
1143
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) {
1168
1169 // argument reduction
1170 final FieldArgumentReduction<FieldComplex<T>> ar = new FieldArgumentReduction<>(phi, m, LegendreEllipticIntegral::bigD);
1171
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);
1176
1177 // combine complete and incomplete parts
1178 return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);
1179
1180 }
1181
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) {
1202
1203 // argument reduction
1204 final DoubleArgumentReduction ar = new DoubleArgumentReduction(phi, m, parameter -> bigPi(n, parameter));
1205
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;
1213
1214 // combine complete and incomplete parts
1215 return ar.negate ? ar.complete - incomplete : ar.complete + incomplete;
1216
1217 }
1218
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) {
1240
1241 // argument reduction
1242 final FieldArgumentReduction<T> ar = new FieldArgumentReduction<>(phi, m, parameter -> bigPi(n, parameter));
1243
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));
1251
1252 // combine complete and incomplete parts
1253 return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);
1254
1255 }
1256
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) {
1283
1284 // argument reduction
1285 final FieldArgumentReduction<Complex> ar = new FieldArgumentReduction<>(phi, m, parameter -> bigPi(n, parameter));
1286
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));
1294
1295 // combine complete and incomplete parts
1296 return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);
1297
1298 }
1299
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 }
1333
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) {
1363
1364 // argument reduction
1365 final FieldArgumentReduction<FieldComplex<T>> ar = new FieldArgumentReduction<>(phi, m, parameter -> bigPi(n, parameter));
1366
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));
1374
1375 // combine complete and incomplete parts
1376 return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);
1377
1378 }
1379
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 }
1418
1419 /** Argument reduction for an incomplete integral. */
1420 private static class DoubleArgumentReduction {
1421
1422 /** Complete part. */
1423 private final double complete;
1424
1425 /** Squared cosecant of the Jacobi amplitude. */
1426 private final double csc2;
1427
1428 /** Indicator for negated Jacobi amplitude. */
1429 private boolean negate;
1430
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 }
1443
1444 }
1445
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>> {
1450
1451 /** Complete part. */
1452 private final T complete;
1453
1454 /** Squared cosecant of the Jacobi amplitude. */
1455 private final T csc2;
1456
1457 /** Indicator for negated Jacobi amplitude. */
1458 private boolean negate;
1459
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 }
1472
1473 }
1474
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> {
1479
1480 /** Parameter. */
1481 private final T m;
1482
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 }
1489
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 }
1497
1498 }
1499
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> {
1504
1505 /** Parameter. */
1506 private final T m;
1507
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 }
1514
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 }
1522
1523 }
1524
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> {
1529
1530 /** Elliptic characteristic. */
1531 private final T n;
1532
1533 /** Parameter. */
1534 private final T m;
1535
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 }
1544
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 }
1554
1555 }
1556 }