1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.hipparchus.special.elliptic.legendre;
18
19 import org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.analysis.CalculusFieldUnivariateFunction;
21 import org.hipparchus.exception.MathIllegalStateException;
22 import org.hipparchus.special.elliptic.carlson.CarlsonEllipticIntegral;
23 import org.hipparchus.util.FastMath;
24 import org.hipparchus.util.MathArrays;
25 import org.hipparchus.util.MathUtils;
26 import org.junit.jupiter.api.Test;
27
28 import static org.junit.jupiter.api.Assertions.assertEquals;
29 import static org.junit.jupiter.api.Assertions.assertTrue;
30
31 public abstract class LegendreEllipticIntegralAbstractComplexTest<T extends CalculusFieldElement<T>> {
32
33 protected abstract T buildComplex(double realPart);
34 protected abstract T buildComplex(double realPart, double imaginaryPart);
35 protected abstract T K(T m);
36 protected abstract T Kprime(T m);
37 protected abstract T F(T phi, T m);
38 protected abstract T integratedF(T phi, T m);
39 protected abstract T E(T m);
40 protected abstract T E(T phi, T m);
41 protected abstract T integratedE(T phi, T m);
42 protected abstract T D(T m);
43 protected abstract T D(T phi, T m);
44 protected abstract T Pi(T n, T m);
45 protected abstract T Pi(T n, T phi, T m);
46 protected abstract T integratedPi(T n, T phi, T m);
47 protected abstract T integrate(int maxEval, CalculusFieldUnivariateFunction<T> f, T start, T end);
48 protected abstract T integrate(int maxEval, CalculusFieldUnivariateFunction<T> f, T start, T middle, T end);
49
50 private void check(double expectedReal, double expectedImaginary, T result, double tol) {
51 assertEquals(0, buildComplex(expectedReal, expectedImaginary).subtract(result).norm(), tol);
52 }
53
54 @Test
55 public void testNoConvergence() {
56 assertTrue(K(buildComplex(Double.NaN)).isNaN());
57 }
58
59 @Test
60 public void testComplementary() {
61 for (double m = 0.01; m < 1; m += 0.01) {
62 T k1 = K(buildComplex(m));
63 T k2 = Kprime(buildComplex(1 - m));
64 assertEquals(k1.getReal(), k2.getReal(), FastMath.ulp(k1).getReal());
65 }
66 }
67
68 @Test
69 public void testAbramowitzStegunExample3() {
70 T k = K(buildComplex(80.0 / 81.0));
71 assertEquals(3.591545001, k.getReal(), 2.0e-9);
72 }
73
74 public void testAbramowitzStegunExample4() {
75 check(1.019106060, 0.0, E(buildComplex(80.0 / 81.0)), 2.0e-8);
76 }
77
78 @Test
79 public void testAbramowitzStegunExample8() {
80 final T m = buildComplex(1.0 / 5.0);
81 check(1.115921, 0.0, F(buildComplex(FastMath.acos(FastMath.sqrt(2) / 3.0)), m), 1.0e-6);
82 check(0.800380, 0.0, F(buildComplex(FastMath.acos(FastMath.sqrt(2) / 2.0)), m), 1.0e-6);
83 }
84
85 @Test
86 public void testAbramowitzStegunExample9() {
87 final T m = buildComplex(1.0 / 2.0);
88 check(1.854075, 0.0, F(buildComplex(MathUtils.SEMI_PI), m), 1.0e-6);
89 check(0.535623, 0.0, F(buildComplex(FastMath.PI / 6.0), m), 1.0e-6);
90 }
91
92 @Test
93 public void testAbramowitzStegunExample10() {
94 final T m = buildComplex(4.0 / 5.0);
95 check(0.543604, 0.0, F(buildComplex(FastMath.PI / 6.0), m), 1.0e-6);
96 }
97
98 @Test
99 public void testAbramowitzStegunExample14() {
100 final T k = buildComplex(3.0 / 5.0);
101 check(0.80904, 0.0, E(buildComplex(FastMath.asin(FastMath.sqrt(5.0) / 3.0)), k.multiply(k)), 1.0e-5);
102 check(0.41192, 0.0, E(buildComplex(FastMath.asin(5.0 / (3.0 * FastMath.sqrt(17.0)))), k.multiply(k)), 1.0e-5);
103 }
104
105 @Test
106 public void testAbramowitzStegunTable175() {
107 final double sinAlpha1 = FastMath.sin(FastMath.toRadians(32));
108 check(0.26263487, 0.0, F(buildComplex(FastMath.toRadians(15)), buildComplex(sinAlpha1 * sinAlpha1)), 1.0e-8);
109 final double sinAlpha2 = FastMath.sin(FastMath.toRadians(46));
110 check(1.61923762, 0.0, F(buildComplex(FastMath.toRadians(80)), buildComplex(sinAlpha2 * sinAlpha2)), 1.0e-8);
111 }
112
113 @Test
114 public void testAbramowitzStegunTable176() {
115 final double sinAlpha1 = FastMath.sin(FastMath.toRadians(64));
116 check(0.42531712, 0.0, E(buildComplex(FastMath.toRadians(25)), buildComplex(sinAlpha1 * sinAlpha1)), 1.0e-8);
117 final double sinAlpha2 = FastMath.sin(FastMath.toRadians(76));
118 check(0.96208074, 0.0, E(buildComplex(FastMath.toRadians(70)), buildComplex(sinAlpha2 * sinAlpha2)), 1.0e-8);
119 }
120
121 @Test
122 public void testAbramowitzStegunTable179() {
123 final double sinAlpha1 = FastMath.sin(FastMath.toRadians(15));
124 check(1.62298, 0.0, Pi(buildComplex(0.4), buildComplex(FastMath.toRadians(75)), buildComplex(sinAlpha1 * sinAlpha1)), 1.0e-5);
125 final double sinAlpha2 = FastMath.sin(FastMath.toRadians(60));
126 check(1.03076, 0.0, Pi(buildComplex(0.8), buildComplex(FastMath.toRadians(45)), buildComplex(sinAlpha2 * sinAlpha2)), 1.0e-5);
127 final double sinAlpha3 = FastMath.sin(FastMath.toRadians(15));
128 check(2.79990, 0.0, Pi(buildComplex(0.9), buildComplex(FastMath.toRadians(75)), buildComplex(sinAlpha3 * sinAlpha3)), 1.0e-5);
129 }
130
131 @Test
132 public void testCompleteVsIncompleteF() {
133 for (double m = 0.01; m < 1; m += 0.01) {
134 double complete = K(buildComplex(m)).getReal();
135 double incomplete = F(buildComplex(MathUtils.SEMI_PI), buildComplex(m)).getReal();
136 assertEquals(complete, incomplete, FastMath.ulp(complete));
137 }
138 }
139
140 @Test
141 public void testCompleteVsIncompleteE() {
142 for (double m = 0.01; m < 1; m += 0.01) {
143 double complete = E(buildComplex(m)).getReal();
144 double incomplete = E(buildComplex(MathUtils.SEMI_PI), buildComplex(m)).getReal();
145 assertEquals(complete, incomplete, 4 * FastMath.ulp(complete));
146 }
147 }
148
149 @Test
150 public void testCompleteVsIncompleteD() {
151 for (double m = 0.01; m < 1; m += 0.01) {
152 double complete = D(buildComplex(m)).getReal();
153 double incomplete = D(buildComplex(MathUtils.SEMI_PI), buildComplex(m)).getReal();
154 assertEquals(complete, incomplete, FastMath.ulp(complete));
155 }
156 }
157
158 @Test
159 public void testCompleteVsIncompletePi() {
160 for (double alpha2 = 0.01; alpha2 < 1; alpha2 += 0.01) {
161 for (double m = 0.01; m < 1; m += 0.01) {
162 double complete = Pi(buildComplex(alpha2), buildComplex(m)).getReal();
163 double incomplete = Pi(buildComplex(alpha2), buildComplex(MathUtils.SEMI_PI), buildComplex(m)).getReal();
164 assertEquals(complete, incomplete, FastMath.ulp(complete));
165 }
166 }
167 }
168
169 @Test
170 public void testNomeSmallParameter() {
171 assertEquals(5.9375e-18, LegendreEllipticIntegral.nome(buildComplex(0.95e-16)).getReal(), 1.0e-50);
172 }
173
174 @Test
175 public void testIntegralsSmallParameter() {
176 assertEquals(7.8539816428e-10, K(buildComplex(2.0e-9)).getReal() - MathUtils.SEMI_PI, 1.0e-15);
177 }
178
179 @Test
180 public void testPrecomputedDelta() {
181
182 T n = buildComplex(3.4, -1.3);
183 T m = buildComplex(0.2, 0.6);
184 T phi = buildComplex(1.2, 0.08);
185 T ref = buildComplex(0.48657872913710487, -0.9325310177613093);
186 assertEquals(0.0, Pi(n, phi, m).subtract(ref).getReal(), 1.0e-15);
187
188
189 final T csc = phi.sin().reciprocal();
190 final T csc2 = csc.multiply(csc);
191 final T cM1 = csc2.subtract(1);
192 final T cMm = csc2.subtract(m);
193 final T cMn = csc2.subtract(n);
194 final T pinphim = CarlsonEllipticIntegral.rF(cM1, cMm, csc2).
195 add(CarlsonEllipticIntegral.rJ(cM1, cMm, csc2, cMn).multiply(n).divide(3));
196 assertEquals(0.0, pinphim.subtract(ref).getReal(), 1.0e-15);
197
198 }
199
200 @Test
201 public void testIncompleteDifferenceA() {
202 final T phi = buildComplex(1.2, 0.75);
203 final T m = buildComplex(0.2, 0.6);
204 final T ref = F(phi, m).
205 subtract(E(phi, m)).
206 divide(m);
207 final T integrated = integrate(100000, new Difference<>(m), buildComplex(1.0e-10, 1.0e-10), phi);
208 assertEquals(0.0, integrated.subtract(ref).norm(), 2.0e-10);
209 assertEquals(0.0, D(phi, m).subtract(ref).norm(), 1.0e-10);
210 }
211
212 @Test
213 public void testIncompleteDifferenceB() {
214 final T phi = buildComplex(1.2, 0.0);
215 final T m = buildComplex(2.3, -1.5);
216 final T ref = F(phi, m).
217 subtract(E(phi, m)).
218 divide(m);
219 final T integrated = integrate(100000, new Difference<>(m), buildComplex(1.0e-10, 1.0e-10), phi);
220 assertEquals(0.0, integrated.subtract(ref).norm(), 2.0e-10);
221 assertEquals(0.0, D(phi, m).subtract(ref).norm(), 1.0e-10);
222 }
223
224 @Test
225 public void testIncompleteDifferenceC() {
226 final T phi = buildComplex(3, 2.5);
227 final T m = buildComplex(2.3, -1.5);
228 final T ref = F(phi, m).subtract(E(phi, m)).divide(m);
229
230
231 final T integrated = integrate(100000, new Difference<>(m), buildComplex(1.0e-12, 1.0e-12), buildComplex(0, -1.5), phi);
232 assertEquals(0.0, integrated.subtract(ref).norm(), 2.0e-10);
233 assertEquals(0.0, D(phi, m).subtract(ref).norm(), 1.0e-10);
234 }
235
236 @Test
237 public void testIncompleteDifferenceD() {
238 final T phi = buildComplex(-0.4, 2.5);
239 final T m = buildComplex(2.3, -1.5);
240 final T ref = F(phi, m).
241 subtract(E(phi, m)).
242 divide(m);
243 final T integrated = integrate(100000, new Difference<>(m), buildComplex(1.0e-10, 1.0e-10), phi);
244 assertEquals(0.0, integrated.subtract(ref).norm(), 2.0e-10);
245 assertEquals(0.0, D(phi, m).subtract(ref).norm(), 1.0e-10);
246 }
247
248 @Test
249 public void testIncompleteFirstKindA() {
250 final T phi = buildComplex(1.2, 0.75);
251 final T m = buildComplex(0.2, 0.6);
252 final T ref = buildComplex(1.00265860821563927579252866, 0.80128721521822408811217);
253 assertEquals(0.0, integratedF(phi, m).subtract(ref).norm(), 2.0e-10);
254 assertEquals(0.0, F(phi, m).subtract(ref).norm(), 1.0e-10);
255 }
256
257 @Test
258 public void testIncompleteFirstKindB() {
259 final T phi = buildComplex(1.2, 0.0);
260 final T m = buildComplex(2.3, -1.5);
261 final T ref = buildComplex(1.04335840461807753156026488, -0.5872679121672512828049797);
262 assertEquals(0.0, integratedF(phi, m).subtract(ref).norm(), 2.0e-10);
263 assertEquals(0.0, F(phi, m).subtract(ref).norm(), 1.0e-10);
264 }
265
266 @Test
267 public void testIncompleteFirstKindC() {
268 final T phi = buildComplex(-0.4, 2.5);
269 final T m = buildComplex(2.3, -1.5);
270 final T ref = buildComplex(-0.20646268947416273887690961, 1.0927692344374984107332330625089);
271 assertEquals(0.0, integratedF(phi, m).subtract(ref).norm(), 2.0e-10);
272 assertEquals(0.0, F(phi, m).subtract(ref).norm(), 1.0e-10);
273 }
274
275 @Test
276 public void testIncompleteSecondKindA() {
277 final T phi = buildComplex(1.2, 0.75);
278 final T m = buildComplex(0.2, 0.6);
279 final T ref = buildComplex(1.4103674846223375296500, 0.644849758860533700396);
280 assertEquals(0.0, integratedE(phi, m).subtract(ref).norm(), 2.0e-10);
281 assertEquals(0.0, E(phi, m).subtract(ref).norm(), 1.0e-10);
282 }
283
284 @Test
285 public void testIncompleteSecondKindB() {
286 final T phi = buildComplex(1.2, 0.0);
287 final T m = buildComplex(2.3, -1.5);
288 final T ref = buildComplex(0.8591316843513079270009549421, 0.55423174445992167002660);
289 assertEquals(0.0, integratedE(phi, m).subtract(ref).norm(), 2.0e-10);
290 assertEquals(0.0, E(phi, m).subtract(ref).norm(), 1.0e-10);
291 }
292
293 @Test
294 public void testIncompleteSecondKindC() {
295 final T phi = buildComplex(-0.4, 2.5);
296 final T m = buildComplex(2.3, -1.5);
297 final T ref = buildComplex(-1.68645030068870706703580773597, 9.176675281683098106653799);
298 assertEquals(0.0, integratedE(phi, m).subtract(ref).norm(), 2.0e-10);
299 assertEquals(0.0, E(phi, m).subtract(ref).norm(), 1.0e-10);
300 }
301
302 @Test
303 public void testIncompleteThirdKind() {
304 final T n = buildComplex(3.4, -1.3);
305 final T m = buildComplex(0.2, 0.6);
306 final T[][] references = MathArrays.buildArray(n.getField(), 38, 2);
307 references[ 0][0] = buildComplex(1.2, -1.5); references[ 0][1] = buildComplex( 0.03351171124667249063, -0.57566536173018225078);
308 references[ 1][0] = buildComplex(1.2, -1.4); references[ 1][1] = buildComplex( 0.03644476655784750591, -0.57331323414589059064);
309 references[ 2][0] = buildComplex(1.2, -1.389190765); references[ 2][1] = buildComplex( 0.03682595624642736804, -0.57302208430696377063);
310 references[ 3][0] = buildComplex(1.2, -1.3); references[ 3][1] = buildComplex( 0.04057582076289887060, -0.57031938416517331851);
311 references[ 4][0] = buildComplex(1.2, -1.2); references[ 4][1] = buildComplex( 0.04640620250501646778, -0.56663563199204550096);
312 references[ 5][0] = buildComplex(1.2, -1.066681968); references[ 5][1] = buildComplex( 0.05798102095188363268, -0.56085091376439940662);
313 references[ 6][0] = buildComplex(1.2, -1.066681967); references[ 6][1] = buildComplex( 0.16640600620018822283, -0.80852909755613891276);
314 references[ 7][0] = buildComplex(1.2, -1.0); references[ 7][1] = buildComplex( 0.17433162130249770262, -0.80552550782381317481);
315 references[ 8][0] = buildComplex(1.2, -0.75); references[ 8][1] = buildComplex( 0.21971749684893289198, -0.79853199575576671312);
316 references[ 9][0] = buildComplex(1.2, -0.5); references[ 9][1] = buildComplex( 0.28863330914614968227, -0.80811912506287693026);
317 references[10][0] = buildComplex(1.2, 0.00); references[10][1] = buildComplex( 0.46199936298130610116, -0.90668901748978345243);
318 references[11][0] = buildComplex(1.2, 0.05); references[11][1] = buildComplex( 0.47768146729143941350, -0.92262785201399518938);
319 references[12][0] = buildComplex(1.2, 0.07); references[12][1] = buildComplex( 0.48365870488290366257, -0.92920571368375179065);
320 references[13][0] = buildComplex(1.2, 0.08); references[13][1] = buildComplex( 0.48657872913710538725, -0.93253101776130872023);
321 references[14][0] = buildComplex(1.2, 0.085); references[14][1] = buildComplex( 0.48802108187307861659, -0.93420195155547102730);
322 references[15][0] = buildComplex(1.2, 0.085181); references[15][1] = buildComplex( 0.48807307162394529967, -0.93426253841598734629);
323 references[16][0] = buildComplex(1.2, 0.085182); references[16][1] = buildComplex( 1.10372915320771558982, -2.71126044767925836757);
324 references[17][0] = buildComplex(1.2, 0.08524491); references[17][1] = buildComplex( 1.10374721953636884577, -2.71128150740577828433);
325 references[18][0] = buildComplex(1.2, 0.08524492); references[18][1] = buildComplex(-1.35887595515643636904, 4.39670878728780118558);
326 references[19][0] = buildComplex(1.2, 0.08524501); references[19][1] = buildComplex(-1.35887592931182948115, 4.39670875715884785695);
327 references[20][0] = buildComplex(1.2, 0.08524502); references[20][1] = buildComplex(-0.12756433765799251204, 0.84271360479056584862);
328 references[21][0] = buildComplex(1.2, 0.08524505); references[21][1] = buildComplex(-0.12756432904312455421, 0.84271359474758096952);
329 references[22][0] = buildComplex(1.2, 0.0852451); references[22][1] = buildComplex(-0.12756431468501224811, 0.84271357800927242222);
330 references[23][0] = buildComplex(1.2, 0.086); references[23][1] = buildComplex(-0.12734767232640510383, 0.84246080397106417910);
331 references[24][0] = buildComplex(1.2, 0.087); references[24][1] = buildComplex(-0.12706111141088970319, 0.84212577870702519585);
332 references[25][0] = buildComplex(1.2, 0.09); references[25][1] = buildComplex(-0.12620431480082688859, 0.84111948445051909467);
333 references[26][0] = buildComplex(1.2, 0.10); references[26][1] = buildComplex(-0.12337990144324015420, 0.83775255734168295690);
334 references[27][0] = buildComplex(1.2, 0.20); references[27][1] = buildComplex(-0.09796886081364286844, 0.80343164778772696735);
335 references[28][0] = buildComplex(1.2, 0.2049); references[28][1] = buildComplex(-0.09686124781335778609, 0.80173956525364697562);
336 references[29][0] = buildComplex(1.2, 0.2051631601); references[29][1] = buildComplex(-0.096802132146685100037, 0.80164871118350302421);
337 references[30][0] = buildComplex(1.2, 0.2051631602); references[30][1] = buildComplex(-0.096802132124228501156, 0.80164871114897922197);
338 references[31][0] = buildComplex(1.2, 0.24); references[31][1] = buildComplex(-0.08930914015165779431, 0.78965671512935519039);
339 references[32][0] = buildComplex(1.2, 0.2462); references[32][1] = buildComplex(-0.08804459384315873622, 0.78753318617257610859);
340 references[33][0] = buildComplex(1.2, 0.24675781599); references[33][1] = buildComplex(-0.087931839058365500429, 0.78734234011803621580);
341 references[34][0] = buildComplex(1.2, 0.2475); references[34][1] = buildComplex(-0.08778207674401363351, 0.78708847169559542169);
342 references[35][0] = buildComplex(1.2, 0.25); references[35][1] = buildComplex(-0.08727979375768656571, 0.78623380867526956107);
343 references[36][0] = buildComplex(1.2, 0.45); references[36][1] = buildComplex(-0.05716147875408999398, 0.72239458160027437391);
344 references[37][0] = buildComplex(1.2, 0.75); references[37][1] = buildComplex(-0.03776232345596223316, 0.65347724469972942256);
345
346 for (int i = 0; i < references.length; ++i) {
347 final T[] ref = references[i];
348 T integrated;
349 try {
350 integrated = integratedPi(n, ref[0], m);
351 } catch (MathIllegalStateException mise) {
352 integrated = buildComplex(Double.NaN);
353 }
354 T carlson = Pi(n, ref[0], m);
355 if (i < 2) {
356
357 assertTrue(true);
358 } else if (i == 2) {
359
360 assertTrue(integrated.isNaN());
361 } else if (i < 6) {
362
363
364 assertEquals(0.0, carlson.subtract(integrated).norm(), 4.1e-7);
365 } else if (i < 16) {
366
367 assertEquals(0.0, integrated.subtract(ref[1]).norm(), 1.0e-10);
368 assertEquals(0.0, carlson.subtract(ref[1]).norm(), 1.0e-10);
369 } else if (i < 30) {
370
371
372 assertEquals(0.0, carlson.subtract(integrated).norm(), 2.0e-6);
373 } else if (i < 35) {
374
375 assertTrue(true);
376 } else {
377
378
379 assertEquals(0.0, integrated.subtract(ref[1]).norm(), 2.5e-7);
380 }
381 }
382 }
383
384 private static class Difference<T extends CalculusFieldElement<T>> implements CalculusFieldUnivariateFunction<T> {
385
386 final T m;
387
388 Difference(final T m) {
389 this.m = m;
390 }
391
392 public T value(final T theta) {
393 final T sin = theta.sin();
394 final T sin2 = sin.multiply(sin);
395 return sin2.divide(sin2.multiply(m).negate().add(1).sqrt());
396 }
397
398 }
399
400 }