View Javadoc
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 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         // no argument reduction and no precomputed delta
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         // we have to use a specific path to get the correct result
228         // integrating over a single straight line gives a completely wrong result
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                 // integration, Carlson and Wolfram Alpha all give different results
355                 Assert.assertTrue(true);
356             } else if (i == 2) {
357                 // integration hits the pole
358                 Assert.assertTrue(integrated.isNaN());
359             } else if (i < 6) {
360                 // integration and Carlson agree and are most probably right
361                 // Wolfram Alpha gives a different result which seems to be wrong
362                 Assert.assertEquals(0.0, carlson.subtract(integrated).norm(), 4.1e-7);
363             } else if (i < 16) {
364                 // integration, Carlson and Wolfram Alpha all agree and are most probably right
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                 // integration and Carlson agree and are most probably right
369                 // Wolfram Alpha gives a different result which seems to be wrong
370                 Assert.assertEquals(0.0, carlson.subtract(integrated).norm(), 2.0e-6);
371             } else if (i < 35) {
372                 // integration, Carlson and Wolfram Alpha all give different results
373                 Assert.assertTrue(true);
374             } else {
375                 // integration and Wolfram Alpha agree and are most probably right
376                 // Carlson gives a different result which seems to be wrong
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 }