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