1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.hipparchus.special.elliptic.jacobi;
18
19 import org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.Field;
21 import org.hipparchus.analysis.differentiation.UnivariateDerivative1;
22 import org.hipparchus.dfp.Dfp;
23 import org.hipparchus.dfp.DfpField;
24 import org.hipparchus.util.Binary64Field;
25 import org.hipparchus.util.FastMath;
26 import org.junit.jupiter.api.Test;
27
28 import java.util.function.DoubleFunction;
29 import java.util.function.Function;
30
31 import static org.junit.jupiter.api.Assertions.assertEquals;
32 import static org.junit.jupiter.api.Assertions.assertTrue;
33
34 class FieldJacobiEllipticTest {
35
36 @Test
37 void testCircular() {
38 doTestCircular(Binary64Field.getInstance());
39 }
40
41 @Test
42 void testHyperbolic() {
43 doTestHyperbolic(Binary64Field.getInstance());
44 }
45
46 @Test
47 void testNoConvergence() {
48 doTestNoConvergence(Binary64Field.getInstance());
49 }
50
51 @Test
52 void testNegativeParameter() {
53 doTestNegativeParameter(Binary64Field.getInstance());
54 }
55
56 @Test
57 void testAbramowitzStegunExample1() {
58 doTestAbramowitzStegunExample1(Binary64Field.getInstance());
59 }
60
61 @Test
62 void testAbramowitzStegunExample2() {
63 doTestAbramowitzStegunExample2(Binary64Field.getInstance());
64 }
65
66 @Test
67 void testAbramowitzStegunExample3() {
68 doTestAbramowitzStegunExample3(Binary64Field.getInstance());
69 }
70
71 @Test
72 void testAbramowitzStegunExample4() {
73 doTestAbramowitzStegunExample4(Binary64Field.getInstance());
74 }
75
76 @Test
77 void testAbramowitzStegunExample5() {
78 doTestAbramowitzStegunExample5(Binary64Field.getInstance());
79 }
80
81 @Test
82 void testAbramowitzStegunExample7() {
83 doTestAbramowitzStegunExample7(Binary64Field.getInstance());
84 }
85
86 @Test
87 void testAbramowitzStegunExample8() {
88 doTestAbramowitzStegunExample8(Binary64Field.getInstance());
89 }
90
91 @Test
92 void testAbramowitzStegunExample9() {
93 doTestAbramowitzStegunExample9(Binary64Field.getInstance());
94 }
95
96 @Test
97 void testAllFunctions() {
98 doTestAllFunctions(Binary64Field.getInstance());
99 }
100
101 @Test
102 void testHighAccuracy() {
103 final DfpField field = new DfpField(60);
104 FieldJacobiElliptic<Dfp> je = JacobiEllipticBuilder.build(field.newDfp("0.75"));
105 Dfp sn = je.valuesN(field.newDfp("1.3")).sn();
106
107 Dfp reference = field.newDfp("0.8929235150418389265984488063926925504375953835259703680383");
108 assertTrue(sn.subtract(reference).abs().getReal() < 5.0e-58);
109 }
110
111 @Test
112 void testInverseCopolarN() {
113 doTestInverseCopolarN(Binary64Field.getInstance());
114 }
115
116 @Test
117 void testInverseCopolarS() {
118 doTestInverseCopolarS(Binary64Field.getInstance());
119 }
120
121 @Test
122 void testInverseCopolarC() {
123 doTestInverseCopolarC(Binary64Field.getInstance());
124 }
125
126 @Test
127 void testInverseCopolarD() {
128 doTestInverseCopolarD(Binary64Field.getInstance());
129 }
130
131 @Test
132 void testDerivatives() {
133
134 final double m = 0.75;
135 final double m1 = 1 - m;
136 final double u = 1.3;
137
138 JacobiElliptic jeD = JacobiEllipticBuilder.build(m);
139 CopolarN valuesND = jeD.valuesN(u);
140 CopolarC valuesCD = jeD.valuesC(u);
141 CopolarS valuesSD = jeD.valuesS(u);
142 CopolarD valuesDD = jeD.valuesD(u);
143
144 FieldJacobiElliptic<UnivariateDerivative1> jeU = JacobiEllipticBuilder.build(new UnivariateDerivative1(m, 0.0));
145 FieldCopolarN<UnivariateDerivative1> valuesNU = jeU.valuesN(new UnivariateDerivative1(u, 1.0));
146 FieldCopolarC<UnivariateDerivative1> valuesCU = jeU.valuesC(new UnivariateDerivative1(u, 1.0));
147 FieldCopolarS<UnivariateDerivative1> valuesSU = jeU.valuesS(new UnivariateDerivative1(u, 1.0));
148 FieldCopolarD<UnivariateDerivative1> valuesDU = jeU.valuesD(new UnivariateDerivative1(u, 1.0));
149
150
151 assertEquals( valuesND.cn() * valuesND.dn(), valuesNU.sn().getFirstDerivative(), 2.0e-15);
152 assertEquals(-1 * valuesND.sn() * valuesND.dn(), valuesNU.cn().getFirstDerivative(), 2.0e-15);
153 assertEquals(-m * valuesND.sn() * valuesND.cn(), valuesNU.dn().getFirstDerivative(), 2.0e-15);
154
155 assertEquals(-m1 * valuesDD.sd() * valuesDD.nd(), valuesDU.cd().getFirstDerivative(), 2.0e-15);
156 assertEquals( valuesDD.cd() * valuesDD.nd(), valuesDU.sd().getFirstDerivative(), 2.0e-15);
157 assertEquals( m * valuesDD.sd() * valuesDD.cd(), valuesDU.nd().getFirstDerivative(), 2.0e-15);
158
159 assertEquals( m1 * valuesCD.sc() * valuesCD.nc(), valuesCU.dc().getFirstDerivative(), 2.0e-15);
160 assertEquals( valuesCD.sc() * valuesCD.dc(), valuesCU.nc().getFirstDerivative(), 2.0e-15);
161 assertEquals( valuesCD.dc() * valuesCD.nc(), valuesCU.sc().getFirstDerivative(), 2.0e-15);
162
163 assertEquals(-1 * valuesSD.ds() * valuesSD.cs(), valuesSU.ns().getFirstDerivative(), 2.0e-15);
164 assertEquals(-1 * valuesSD.cs() * valuesSD.ns(), valuesSU.ds().getFirstDerivative(), 2.0e-15);
165 assertEquals(-1 * valuesSD.ns() * valuesSD.ds(), valuesSU.cs().getFirstDerivative(), 2.0e-15);
166
167 }
168
169 private <T extends CalculusFieldElement<T>> void doTestCircular(final Field<T> field) {
170 for (double m : new double[] { -1.0e-10, 0.0, 1.0e-10 }) {
171 final double eps = 3 * FastMath.max(1.0e-14, FastMath.abs(m));
172 final FieldJacobiElliptic<T> je = build(field, m);
173 for (double t = -10; t < 10; t += 0.01) {
174 final FieldCopolarN<T> n = je.valuesN(t);
175 assertEquals(FastMath.sin(t), n.sn().getReal(), eps);
176 assertEquals(FastMath.cos(t), n.cn().getReal(), eps);
177 assertEquals(1.0, n.dn().getReal(), eps);
178 }
179 }
180 }
181
182 private <T extends CalculusFieldElement<T>> void doTestHyperbolic(final Field<T> field) {
183 for (double m1 : new double[] { -1.0e-12, 0.0, 1.0e-12 }) {
184 final double eps = 3 * FastMath.max(1.0e-14, FastMath.abs(m1));
185 final FieldJacobiElliptic<T> je = build(field, 1.0 - m1);
186 for (double t = -3; t < 3; t += 0.01) {
187 final FieldCopolarN<T> n = je.valuesN(t);
188 assertEquals(FastMath.tanh(t), n.sn().getReal(), eps);
189 assertEquals(1.0 / FastMath.cosh(t), n.cn().getReal(), eps);
190 assertEquals(1.0 / FastMath.cosh(t), n.dn().getReal(), eps);
191 }
192 }
193 }
194
195 private <T extends CalculusFieldElement<T>> void doTestNoConvergence(final Field<T> field) {
196 assertTrue(build(field, Double.NaN).valuesS(0.0).cs().isNaN());
197 }
198
199 private <T extends CalculusFieldElement<T>> void doTestNegativeParameter(final Field<T> field) {
200 assertEquals(0.49781366219021166315, build(field, -4.5).valuesN(8.3).sn().getReal(), 1.5e-10);
201 assertEquals(0.86728401215332559984, build(field, -4.5).valuesN(8.3).cn().getReal(), 1.5e-10);
202 assertEquals(1.45436686918553524215, build(field, -4.5).valuesN(8.3).dn().getReal(), 1.5e-10);
203 }
204
205 private <T extends CalculusFieldElement<T>> void doTestAbramowitzStegunExample1(final Field<T> field) {
206
207 assertEquals(-1392.11114434139393839735, build(field, 0.64).valuesC(1.99650).nc().getReal(), 2.8e-10);
208 }
209
210 private <T extends CalculusFieldElement<T>> void doTestAbramowitzStegunExample2(final Field<T> field) {
211 assertEquals(0.996253, build(field, 0.19).valuesN(0.20).dn().getReal(), 1.0e-6);
212 }
213
214 private <T extends CalculusFieldElement<T>> void doTestAbramowitzStegunExample3(final Field<T> field) {
215 assertEquals(0.984056, build(field, 0.81).valuesN(0.20).dn().getReal(), 1.0e-6);
216 }
217
218 private <T extends CalculusFieldElement<T>> void doTestAbramowitzStegunExample4(final Field<T> field) {
219 assertEquals(0.980278, build(field, 0.81).valuesN(0.20).cn().getReal(), 1.0e-6);
220 }
221
222 private <T extends CalculusFieldElement<T>> void doTestAbramowitzStegunExample5(final Field<T> field) {
223 assertEquals(0.60952, build(field, 0.36).valuesN(0.672).sn().getReal(), 1.0e-5);
224 assertEquals(1.1740, build(field, 0.36).valuesC(0.672).dc().getReal(), 1.0e-4);
225 }
226
227 private <T extends CalculusFieldElement<T>> void doTestAbramowitzStegunExample7(final Field<T> field) {
228 assertEquals(1.6918083, build(field, 0.09).valuesS(0.5360162).cs().getReal(), 1.0e-7);
229 }
230
231 private <T extends CalculusFieldElement<T>> void doTestAbramowitzStegunExample8(final Field<T> field) {
232 assertEquals(0.56458, build(field, 0.5).valuesN(0.61802).sn().getReal(), 1.0e-5);
233 }
234
235 private <T extends CalculusFieldElement<T>> void doTestAbramowitzStegunExample9(final Field<T> field) {
236 assertEquals(0.68402, build(field, 0.5).valuesC(0.61802).sc().getReal(), 1.0e-5);
237 }
238
239 private <T extends CalculusFieldElement<T>> void doTestAllFunctions(final Field<T> field) {
240
241
242
243
244 final double u = 1.4;
245 final double m = 0.7;
246 final double[] reference = {
247 0.92516138673582827365, 0.37957398289798418747, 0.63312991237590996850,
248 0.41027866958131945027, 0.68434537093007175683, 1.08089249544689572795,
249 1.66800134071905681841, 2.63453251554589286796, 2.43736775548306830513,
250 1.57945467502452678756, 1.46125047743207819361, 0.59951990180590090343
251 };
252 final FieldJacobiElliptic<T> je = build(field, m);
253 assertEquals(reference[ 0], je.valuesN(u).sn().getReal(), 4 * FastMath.ulp(reference[ 0]));
254 assertEquals(reference[ 1], je.valuesN(u).cn().getReal(), 4 * FastMath.ulp(reference[ 1]));
255 assertEquals(reference[ 2], je.valuesN(u).dn().getReal(), 4 * FastMath.ulp(reference[ 2]));
256 assertEquals(reference[ 3], je.valuesS(u).cs().getReal(), 4 * FastMath.ulp(reference[ 3]));
257 assertEquals(reference[ 4], je.valuesS(u).ds().getReal(), 4 * FastMath.ulp(reference[ 4]));
258 assertEquals(reference[ 5], je.valuesS(u).ns().getReal(), 4 * FastMath.ulp(reference[ 5]));
259 assertEquals(reference[ 6], je.valuesC(u).dc().getReal(), 4 * FastMath.ulp(reference[ 6]));
260 assertEquals(reference[ 7], je.valuesC(u).nc().getReal(), 4 * FastMath.ulp(reference[ 7]));
261 assertEquals(reference[ 8], je.valuesC(u).sc().getReal(), 4 * FastMath.ulp(reference[ 8]));
262 assertEquals(reference[ 9], je.valuesD(u).nd().getReal(), 4 * FastMath.ulp(reference[ 9]));
263 assertEquals(reference[10], je.valuesD(u).sd().getReal(), 4 * FastMath.ulp(reference[10]));
264 assertEquals(reference[11], je.valuesD(u).cd().getReal(), 4 * FastMath.ulp(reference[11]));
265 }
266
267 private <T extends CalculusFieldElement<T>> FieldJacobiElliptic<T> build(final Field<T> field, final double m) {
268 return JacobiEllipticBuilder.build(field.getZero().newInstance(m));
269 }
270
271 private <T extends CalculusFieldElement<T>> void doTestInverseCopolarN(final Field<T> field) {
272 final double m = 0.7;
273 final FieldJacobiElliptic<T> je = build(field, m);
274 doTestInverse(-0.80, 0.80, 100, field, u -> je.valuesN(u).sn(), x -> je.arcsn(x), x -> je.arcsn(x), 1.0e-14);
275 doTestInverse(-1.00, 1.00, 100, field, u -> je.valuesN(u).cn(), x -> je.arccn(x), x -> je.arccn(x), 1.0e-14);
276 doTestInverse( 0.55, 1.00, 100, field, u -> je.valuesN(u).dn(), x -> je.arcdn(x), x -> je.arcdn(x), 1.0e-14);
277 }
278
279 private <T extends CalculusFieldElement<T>> void doTestInverseCopolarS(final Field<T> field) {
280 final double m = 0.7;
281 final FieldJacobiElliptic<T> je = build(field, m);
282 doTestInverse(-2.00, 2.00, 100, field, u -> je.valuesS(u).cs(), x -> je.arccs(x), x -> je.arccs(x), 1.0e-14);
283 doTestInverse( 0.55, 2.00, 100, field, u -> je.valuesS(u).ds(), x -> je.arcds(x), x -> je.arcds(x), 1.0e-14);
284 doTestInverse(-2.00, -0.55, 100, field, u -> je.valuesS(u).ds(), x -> je.arcds(x), x -> je.arcds(x), 1.0e-14);
285 doTestInverse( 1.00, 2.00, 100, field, u -> je.valuesS(u).ns(), x -> je.arcns(x), x -> je.arcns(x), 1.0e-11);
286 doTestInverse(-2.00, -1.00, 100, field, u -> je.valuesS(u).ns(), x -> je.arcns(x), x -> je.arcns(x), 1.0e-11);
287 }
288
289 private <T extends CalculusFieldElement<T>> void doTestInverseCopolarC(final Field<T> field) {
290 final double m = 0.7;
291 final FieldJacobiElliptic<T> je = build(field, m);
292 doTestInverse( 1.00, 2.00, 100, field, u -> je.valuesC(u).dc(), x -> je.arcdc(x), x -> je.arcdc(x), 1.0e-14);
293 doTestInverse(-2.00, -1.00, 100, field, u -> je.valuesC(u).dc(), x -> je.arcdc(x), x -> je.arcdc(x), 1.0e-14);
294 doTestInverse( 1.00, 2.00, 100, field, u -> je.valuesC(u).nc(), x -> je.arcnc(x), x -> je.arcnc(x), 1.0e-14);
295 doTestInverse(-2.00, -1.00, 100, field, u -> je.valuesC(u).nc(), x -> je.arcnc(x), x -> je.arcnc(x), 1.0e-14);
296 doTestInverse(-2.00, 2.00, 100, field, u -> je.valuesC(u).sc(), x -> je.arcsc(x), x -> je.arcsc(x), 1.0e-14);
297 }
298
299 private <T extends CalculusFieldElement<T>> void doTestInverseCopolarD(final Field<T> field) {
300 final double m = 0.7;
301 final FieldJacobiElliptic<T> je = build(field, m);
302 doTestInverse( 1.00, 1.80, 100, field, u -> je.valuesD(u).nd(), x -> je.arcnd(x), x -> je.arcnd(x), 1.0e-14);
303 doTestInverse(-1.80, 1.80, 100, field, u -> je.valuesD(u).sd(), x -> je.arcsd(x), x -> je.arcsd(x), 1.0e-14);
304 doTestInverse(-1.00, 1.00, 100, field, u -> je.valuesD(u).cd(), x -> je.arccd(x), x -> je.arccd(x), 1.0e-14);
305 }
306
307 private <T extends CalculusFieldElement<T>> void doTestInverse(final double xMin, final double xMax, final int n,
308 final Field<T> field,
309 final Function<T, T> direct,
310 final Function<T, T> inverseField,
311 final DoubleFunction<T> inverseDouble,
312 final double tolerance) {
313 for (int i = 0; i < n; ++i) {
314 final T x = field.getZero().newInstance(xMin + i * (xMax - xMin) / (n - 1));
315 final T xFieldRebuilt = direct.apply(inverseField.apply(x));
316 assertEquals(x.getReal(), xFieldRebuilt.getReal(), tolerance);
317 final T xDoubleRebuilt = direct.apply(inverseDouble.apply(x.getReal()));
318 assertEquals(x.getReal(), xDoubleRebuilt.getReal(), tolerance);
319 }
320 }
321
322 }