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.jacobi;
18  
19  import org.hipparchus.analysis.UnivariateFunction;
20  import org.hipparchus.util.FastMath;
21  import org.junit.Assert;
22  import org.junit.Test;
23  
24  public class JacobiEllipticTest {
25  
26      @Test
27      public void testCircular() {
28          for (double m : new double[] { -1.0e-10, 0.0, 1.0e-10 }) {
29           final double eps = 3 * FastMath.max(1.0e-14, FastMath.abs(m));
30              final JacobiElliptic je = JacobiEllipticBuilder.build(m);
31              for (double t = -10; t < 10; t += 0.01) {
32                  final CopolarN n = je.valuesN(t);
33                  Assert.assertEquals(FastMath.sin(t), n.sn(), eps);
34                  Assert.assertEquals(FastMath.cos(t), n.cn(), eps);
35                  Assert.assertEquals(1.0,             n.dn(), eps);
36              }
37          }
38      }
39  
40      @Test
41      public void testHyperbolic() {
42          for (double m1 : new double[] { -1.0e-12, 0.0, 1.0e-12 }) {
43              final double eps = 3 * FastMath.max(1.0e-14, FastMath.abs(m1));
44              final JacobiElliptic je = JacobiEllipticBuilder.build(1.0 - m1);
45              for (double t = -3; t < 3; t += 0.01) {
46                  final CopolarN n = je.valuesN(t);
47                  Assert.assertEquals(FastMath.tanh(t),       n.sn(), eps);
48                  Assert.assertEquals(1.0 / FastMath.cosh(t), n.cn(), eps);
49                  Assert.assertEquals(1.0 / FastMath.cosh(t), n.dn(), eps);
50              }
51          }
52      }
53  
54      @Test
55      public void testNoConvergence() {
56          Assert.assertTrue(Double.isNaN(JacobiEllipticBuilder.build(Double.NaN).valuesS(0.0).cs()));
57      }
58  
59      @Test
60      public void testNegativeParameter() {
61          Assert.assertEquals(0.49781366219021166315, JacobiEllipticBuilder.build(-4.5).valuesN(8.3).sn(), 1.5e-10);
62          Assert.assertEquals(0.86728401215332559984, JacobiEllipticBuilder.build(-4.5).valuesN(8.3).cn(), 1.5e-10);
63          Assert.assertEquals(1.45436686918553524215, JacobiEllipticBuilder.build(-4.5).valuesN(8.3).dn(), 1.5e-10);
64      }
65  
66      @Test
67      public void testAbramowitzStegunExample1() {
68          // Abramowitz and Stegun give a result of -1667, but Wolfram Alpha gives the following value
69          Assert.assertEquals(-1392.11114434139393839735, JacobiEllipticBuilder.build(0.64).valuesC(1.99650).nc(), 6.0e-10);
70      }
71  
72      @Test
73      public void testAbramowitzStegunExample2() {
74          Assert.assertEquals(0.996253, JacobiEllipticBuilder.build(0.19).valuesN(0.20).dn(), 1.0e-6);
75      }
76  
77      @Test
78      public void testAbramowitzStegunExample3() {
79          Assert.assertEquals(0.984056, JacobiEllipticBuilder.build(0.81).valuesN(0.20).dn(), 1.0e-6);
80      }
81  
82      @Test
83      public void testAbramowitzStegunExample4() {
84          Assert.assertEquals(0.980278, JacobiEllipticBuilder.build(0.81).valuesN(0.20).cn(), 1.0e-6);
85      }
86  
87      @Test
88      public void testAbramowitzStegunExample5() {
89          Assert.assertEquals(0.60952, JacobiEllipticBuilder.build(0.36).valuesN(0.672).sn(), 1.0e-5);
90          Assert.assertEquals(1.1740, JacobiEllipticBuilder.build(0.36).valuesC(0.672).dc(), 1.0e-4);
91      }
92  
93      @Test
94      public void testAbramowitzStegunExample7() {
95          Assert.assertEquals(1.6918083, JacobiEllipticBuilder.build(0.09).valuesS(0.5360162).cs(), 1.0e-7);
96      }
97  
98      @Test
99      public void testAbramowitzStegunExample8() {
100         Assert.assertEquals(0.56458, JacobiEllipticBuilder.build(0.5).valuesN(0.61802).sn(), 1.0e-5);
101     }
102 
103     @Test
104     public void testAbramowitzStegunExample9() {
105         Assert.assertEquals(0.68402, JacobiEllipticBuilder.build(0.5).valuesC(0.61802).sc(), 1.0e-5);
106     }
107 
108     @Test
109     public void testAllFunctions() {
110         // reference was computed from Wolfram Alpha, using the square relations
111         // from Abramowitz and Stegun section 16.9 for the functions Wolfram Alpha
112         // did not understood (i.e. for the sake of validation we did *not* use the
113         // relations from section 16.3 which are used in the Hipparchus implementation)
114         final double u = 1.4;
115         final double m = 0.7;
116         final double[] reference = {
117               0.92516138673582827365, 0.37957398289798418747, 0.63312991237590996850,
118               0.41027866958131945027, 0.68434537093007175683, 1.08089249544689572795,
119               1.66800134071905681841, 2.63453251554589286796, 2.43736775548306830513,
120               1.57945467502452678756, 1.46125047743207819361, 0.59951990180590090343
121         };
122         final JacobiElliptic je = JacobiEllipticBuilder.build(m);
123         Assert.assertEquals(reference[ 0], je.valuesN(u).sn(), 4 * FastMath.ulp(reference[ 0]));
124         Assert.assertEquals(reference[ 1], je.valuesN(u).cn(), 4 * FastMath.ulp(reference[ 1]));
125         Assert.assertEquals(reference[ 2], je.valuesN(u).dn(), 4 * FastMath.ulp(reference[ 2]));
126         Assert.assertEquals(reference[ 3], je.valuesS(u).cs(), 4 * FastMath.ulp(reference[ 3]));
127         Assert.assertEquals(reference[ 4], je.valuesS(u).ds(), 4 * FastMath.ulp(reference[ 4]));
128         Assert.assertEquals(reference[ 5], je.valuesS(u).ns(), 4 * FastMath.ulp(reference[ 5]));
129         Assert.assertEquals(reference[ 6], je.valuesC(u).dc(), 4 * FastMath.ulp(reference[ 6]));
130         Assert.assertEquals(reference[ 7], je.valuesC(u).nc(), 4 * FastMath.ulp(reference[ 7]));
131         Assert.assertEquals(reference[ 8], je.valuesC(u).sc(), 4 * FastMath.ulp(reference[ 8]));
132         Assert.assertEquals(reference[ 9], je.valuesD(u).nd(), 4 * FastMath.ulp(reference[ 9]));
133         Assert.assertEquals(reference[10], je.valuesD(u).sd(), 4 * FastMath.ulp(reference[10]));
134         Assert.assertEquals(reference[11], je.valuesD(u).cd(), 4 * FastMath.ulp(reference[11]));
135     }
136 
137     @Test
138     public void testInverseCopolarN() {
139         final double m = 0.7;
140         final JacobiElliptic je = JacobiEllipticBuilder.build(m);
141         doTestInverse(-0.80,  0.80, 100, u -> je.valuesN(u).sn(), x -> je.arcsn(x), 1.0e-14);
142         doTestInverse(-1.00,  1.00, 100, u -> je.valuesN(u).cn(), x -> je.arccn(x), 1.0e-14);
143         doTestInverse( 0.55,  1.00, 100, u -> je.valuesN(u).dn(), x -> je.arcdn(x), 1.0e-14);
144     }
145 
146     @Test
147     public void testInverseCopolarS() {
148         final double m = 0.7;
149         final JacobiElliptic je = JacobiEllipticBuilder.build(m);
150         doTestInverse(-2.00,  2.00, 100, u -> je.valuesS(u).cs(), x -> je.arccs(x), 1.0e-14);
151         doTestInverse( 0.55,  2.00, 100, u -> je.valuesS(u).ds(), x -> je.arcds(x), 1.0e-14);
152         doTestInverse(-2.00, -0.55, 100, u -> je.valuesS(u).ds(), x -> je.arcds(x), 1.0e-14);
153         doTestInverse( 1.00,  2.00, 100, u -> je.valuesS(u).ns(), x -> je.arcns(x), 1.0e-11);
154         doTestInverse(-2.00, -1.00, 100, u -> je.valuesS(u).ns(), x -> je.arcns(x), 1.0e-11);
155     }
156 
157     @Test
158     public void testInverseCopolarC() {
159         final double m = 0.7;
160         final JacobiElliptic je = JacobiEllipticBuilder.build(m);
161         doTestInverse( 1.00,  2.00, 100, u -> je.valuesC(u).dc(), x -> je.arcdc(x), 1.0e-14);
162         doTestInverse(-2.00, -1.00, 100, u -> je.valuesC(u).dc(), x -> je.arcdc(x), 1.0e-14);
163         doTestInverse( 1.00,  2.00, 100, u -> je.valuesC(u).nc(), x -> je.arcnc(x), 1.0e-14);
164         doTestInverse(-2.00, -1.00, 100, u -> je.valuesC(u).nc(), x -> je.arcnc(x), 1.0e-14);
165         doTestInverse(-2.00,  2.00, 100, u -> je.valuesC(u).sc(), x -> je.arcsc(x), 1.0e-14);
166     }
167 
168     @Test
169     public void testInverseCopolarD() {
170         final double m = 0.7;
171         final JacobiElliptic je = JacobiEllipticBuilder.build(m);
172         doTestInverse( 1.00,  1.80, 100, u -> je.valuesD(u).nd(), x -> je.arcnd(x), 1.0e-14);
173         doTestInverse(-1.80,  1.80, 100, u -> je.valuesD(u).sd(), x -> je.arcsd(x), 1.0e-14);
174         doTestInverse(-1.00,  1.00, 100, u -> je.valuesD(u).cd(), x -> je.arccd(x), 1.0e-14);
175     }
176 
177     private void doTestInverse(final double xMin, final double xMax, final int n,
178                                final UnivariateFunction direct, final UnivariateFunction inverse,
179                                final double tolerance) {
180         for (int i = 0; i < n; ++i) {
181             final double x        = xMin + i * (xMax - xMin) / (n - 1);
182             final double xRebuilt = direct.value(inverse.value(x));
183             Assert.assertEquals(x, xRebuilt, tolerance);
184         }
185     }
186 
187 }