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.CalculusFieldElement;
20 import org.hipparchus.special.elliptic.carlson.CarlsonEllipticIntegral;
21 import org.hipparchus.special.elliptic.legendre.LegendreEllipticIntegral;
22 import org.hipparchus.util.FastMath;
23
24 /** Computation of Jacobi elliptic functions.
25 * The Jacobi elliptic functions are related to elliptic integrals.
26 * @param <T> the type of the field elements
27 * @since 2.0
28 */
29 public abstract class FieldJacobiElliptic<T extends CalculusFieldElement<T>> {
30
31 /** Parameter of the function. */
32 private final T m;
33
34 /** Simple constructor.
35 * @param m parameter of the function
36 */
37 protected FieldJacobiElliptic(final T m) {
38 this.m = m;
39 }
40
41 /** Get the parameter of the function.
42 * @return parameter of the function
43 */
44 public T getM() {
45 return m;
46 }
47
48 /** Evaluate the three principal Jacobi elliptic functions with pole at point n in Glaisher’s Notation.
49 * @param u argument of the functions
50 * @return copolar trio containing the three principal Jacobi
51 * elliptic functions {@code sn(u|m)}, {@code cn(u|m)}, and {@code dn(u|m)}.
52 */
53 public abstract FieldCopolarN<T> valuesN(T u);
54
55 /** Evaluate the three principal Jacobi elliptic functions with pole at point n in Glaisher’s Notation.
56 * @param u argument of the functions
57 * @return copolar trio containing the three principal Jacobi
58 * elliptic functions {@code sn(u|m)}, {@code cn(u|m)}, and {@code dn(u|m)}.
59 */
60 public FieldCopolarN<T> valuesN(final double u) {
61 return valuesN(m.newInstance(u));
62 }
63
64 /** Evaluate the three subsidiary Jacobi elliptic functions with pole at point s in Glaisher’s Notation.
65 * @param u argument of the functions
66 * @return copolar trio containing the three subsidiary Jacobi
67 * elliptic functions {@code cs(u|m)}, {@code ds(u|m)} and {@code ns(u|m)}.
68 */
69 public FieldCopolarS<T> valuesS(final T u) {
70 return new FieldCopolarS<>(valuesN(u));
71 }
72
73 /** Evaluate the three subsidiary Jacobi elliptic functions with pole at point s in Glaisher’s Notation.
74 * @param u argument of the functions
75 * @return copolar trio containing the three subsidiary Jacobi
76 * elliptic functions {@code cs(u|m)}, {@code ds(u|m)} and {@code ns(u|m)}.
77 */
78 public FieldCopolarS<T> valuesS(final double u) {
79 return new FieldCopolarS<>(valuesN(u));
80 }
81
82 /** Evaluate the three subsidiary Jacobi elliptic functions with pole at point c in Glaisher’s Notation.
83 * @param u argument of the functions
84 * @return copolar trio containing the three subsidiary Jacobi
85 * elliptic functions {@code dc(u|m)}, {@code nc(u|m)}, and {@code sc(u|m)}.
86 */
87 public FieldCopolarC<T> valuesC(final T u) {
88 return new FieldCopolarC<>(valuesN(u));
89 }
90
91 /** Evaluate the three subsidiary Jacobi elliptic functions with pole at point c in Glaisher’s Notation.
92 * @param u argument of the functions
93 * @return copolar trio containing the three subsidiary Jacobi
94 * elliptic functions {@code dc(u|m)}, {@code nc(u|m)}, and {@code sc(u|m)}.
95 */
96 public FieldCopolarC<T> valuesC(final double u) {
97 return new FieldCopolarC<>(valuesN(u));
98 }
99
100 /** Evaluate the three subsidiary Jacobi elliptic functions with pole at point d in Glaisher’s Notation.
101 * @param u argument of the functions
102 * @return copolar trio containing the three subsidiary Jacobi
103 * elliptic functions {@code nd(u|m)}, {@code sd(u|m)}, and {@code cd(u|m)}.
104 */
105 public FieldCopolarD<T> valuesD(final T u) {
106 return new FieldCopolarD<>(valuesN(u));
107 }
108
109 /** Evaluate the three subsidiary Jacobi elliptic functions with pole at point d in Glaisher’s Notation.
110 * @param u argument of the functions
111 * @return copolar trio containing the three subsidiary Jacobi
112 * elliptic functions {@code nd(u|m)}, {@code sd(u|m)}, and {@code cd(u|m)}.
113 */
114 public FieldCopolarD<T> valuesD(final double u) {
115 return new FieldCopolarD<>(valuesN(u));
116 }
117
118 /** Evaluate inverse of Jacobi elliptic function sn.
119 * @param x value of Jacobi elliptic function {@code sn(u|m)}
120 * @return u such that {@code x=sn(u|m)}
121 * @since 2.1
122 */
123 public T arcsn(final T x) {
124 // p = n, q = c, r = d, see DLMF 19.25.29 for evaluating Δ(q, p) and Δ(r, p)
125 return arcsp(x, x.getField().getOne().negate(), getM().negate());
126 }
127
128 /** Evaluate inverse of Jacobi elliptic function sn.
129 * @param x value of Jacobi elliptic function {@code sn(u|m)}
130 * @return u such that {@code x=sn(u|m)}
131 * @since 2.1
132 */
133 public T arcsn(final double x) {
134 return arcsn(getM().getField().getZero().newInstance(x));
135 }
136
137 /** Evaluate inverse of Jacobi elliptic function cn.
138 * @param x value of Jacobi elliptic function {@code cn(u|m)}
139 * @return u such that {@code x=cn(u|m)}
140 * @since 2.1
141 */
142 public T arccn(final T x) {
143 // p = c, q = n, r = d, see DLMF 19.25.29 for evaluating Δ(q, p) and Δ(r, q)
144 return arcpqNoDivision(x, x.getField().getOne(), getM().negate());
145 }
146
147 /** Evaluate inverse of Jacobi elliptic function cn.
148 * @param x value of Jacobi elliptic function {@code cn(u|m)}
149 * @return u such that {@code x=cn(u|m)}
150 * @since 2.1
151 */
152 public T arccn(final double x) {
153 return arccn(getM().getField().getZero().newInstance(x));
154 }
155
156 /** Evaluate inverse of Jacobi elliptic function dn.
157 * @param x value of Jacobi elliptic function {@code dn(u|m)}
158 * @return u such that {@code x=dn(u|m)}
159 * @since 2.1
160 */
161 public T arcdn(final T x) {
162 // p = d, q = n, r = c, see DLMF 19.25.29 for evaluating Δ(q, p) and Δ(r, q)
163 return arcpqNoDivision(x, getM(), x.getField().getOne().negate());
164 }
165
166 /** Evaluate inverse of Jacobi elliptic function dn.
167 * @param x value of Jacobi elliptic function {@code dn(u|m)}
168 * @return u such that {@code x=dn(u|m)}
169 * @since 2.1
170 */
171 public T arcdn(final double x) {
172 return arcdn(getM().getField().getZero().newInstance(x));
173 }
174
175 /** Evaluate inverse of Jacobi elliptic function cs.
176 * @param x value of Jacobi elliptic function {@code cs(u|m)}
177 * @return u such that {@code x=cs(u|m)}
178 * @since 2.1
179 */
180 public T arccs(final T x) {
181 // p = c, q = n, r = d, see DLMF 19.25.29 for evaluating Δ(q, p) and Δ(r, p)
182 return arcps(x, x.getField().getOne(), getM().subtract(1).negate());
183 }
184
185 /** Evaluate inverse of Jacobi elliptic function cs.
186 * @param x value of Jacobi elliptic function {@code cs(u|m)}
187 * @return u such that {@code x=cs(u|m)}
188 * @since 2.1
189 */
190 public T arccs(final double x) {
191 return arccs(getM().getField().getZero().newInstance(x));
192 }
193
194 /** Evaluate inverse of Jacobi elliptic function ds.
195 * @param x value of Jacobi elliptic function {@code ds(u|m)}
196 * @return u such that {@code x=ds(u|m)}
197 * @since 2.1
198 */
199 public T arcds(final T x) {
200 // p = d, q = c, r = n, see DLMF 19.25.29 for evaluating Δ(q, p) and Δ(r, p)
201 return arcps(x, getM().subtract(1), getM());
202 }
203
204 /** Evaluate inverse of Jacobi elliptic function ds.
205 * @param x value of Jacobi elliptic function {@code ds(u|m)}
206 * @return u such that {@code x=ds(u|m)}
207 * @since 2.1
208 */
209 public T arcds(final double x) {
210 return arcds(getM().getField().getZero().newInstance(x));
211 }
212
213 /** Evaluate inverse of Jacobi elliptic function ns.
214 * @param x value of Jacobi elliptic function {@code ns(u|m)}
215 * @return u such that {@code x=ns(u|m)}
216 * @since 2.1
217 */
218 public T arcns(final T x) {
219 // p = n, q = c, r = d, see DLMF 19.25.29 for evaluating Δ(q, p) and Δ(r, p)
220 return arcps(x, x.getField().getOne().negate(), getM().negate());
221 }
222
223 /** Evaluate inverse of Jacobi elliptic function ns.
224 * @param x value of Jacobi elliptic function {@code ns(u|m)}
225 * @return u such that {@code x=ns(u|m)}
226 * @since 2.1
227 */
228 public T arcns(final double x) {
229 return arcns(getM().getField().getZero().newInstance(x));
230 }
231
232 /** Evaluate inverse of Jacobi elliptic function dc.
233 * @param x value of Jacobi elliptic function {@code dc(u|m)}
234 * @return u such that {@code x=dc(u|m)}
235 * @since 2.1
236 */
237 public T arcdc(final T x) {
238 // p = d, q = c, r = n, see DLMF 19.25.29 for evaluating Δ(q, p) and Δ(r, q)
239 return arcpq(x, getM().subtract(1), x.getField().getOne());
240 }
241
242 /** Evaluate inverse of Jacobi elliptic function dc.
243 * @param x value of Jacobi elliptic function {@code dc(u|m)}
244 * @return u such that {@code x=dc(u|m)}
245 * @since 2.1
246 */
247 public T arcdc(final double x) {
248 return arcdc(getM().getField().getZero().newInstance(x));
249 }
250
251 /** Evaluate inverse of Jacobi elliptic function nc.
252 * @param x value of Jacobi elliptic function {@code nc(u|m)}
253 * @return u such that {@code x=nc(u|m)}
254 * @since 2.1
255 */
256 public T arcnc(final T x) {
257 // p = n, q = c, r = d, see DLMF 19.25.29 for evaluating Δ(q, p) and Δ(r, q)
258 return arcpq(x, x.getField().getOne().negate(), getM().subtract(1).negate());
259 }
260
261 /** Evaluate inverse of Jacobi elliptic function nc.
262 * @param x value of Jacobi elliptic function {@code nc(u|m)}
263 * @return u such that {@code x=nc(u|m)}
264 * @since 2.1
265 */
266 public T arcnc(final double x) {
267 return arcnc(getM().getField().getZero().newInstance(x));
268 }
269
270 /** Evaluate inverse of Jacobi elliptic function sc.
271 * @param x value of Jacobi elliptic function {@code sc(u|m)}
272 * @return u such that {@code x=sc(u|m)}
273 * @since 2.1
274 */
275 public T arcsc(final T x) {
276 // p = c, q = n, r = d, see DLMF 19.25.29 for evaluating Δ(q, p) and Δ(r, p)
277 return arcsp(x, x.getField().getOne(), getM().subtract(1).negate());
278 }
279
280 /** Evaluate inverse of Jacobi elliptic function sc.
281 * @param x value of Jacobi elliptic function {@code sc(u|m)}
282 * @return u such that {@code x=sc(u|m)}
283 * @since 2.1
284 */
285 public T arcsc(final double x) {
286 return arcsc(getM().getField().getZero().newInstance(x));
287 }
288
289 /** Evaluate inverse of Jacobi elliptic function nd.
290 * @param x value of Jacobi elliptic function {@code nd(u|m)}
291 * @return u such that {@code x=nd(u|m)}
292 * @since 2.1
293 */
294 public T arcnd(final T x) {
295 // p = n, q = d, r = c, see DLMF 19.25.29 for evaluating Δ(q, p) and Δ(r, q)
296 return arcpq(x, getM().negate(), getM().subtract(1));
297 }
298
299 /** Evaluate inverse of Jacobi elliptic function nd.
300 * @param x value of Jacobi elliptic function {@code nd(u|m)}
301 * @return u such that {@code x=nd(u|m)}
302 * @since 2.1
303 */
304 public T arcnd(final double x) {
305 return arcnd(getM().getField().getZero().newInstance(x));
306 }
307
308 /** Evaluate inverse of Jacobi elliptic function sd.
309 * @param x value of Jacobi elliptic function {@code sd(u|m)}
310 * @return u such that {@code x=sd(u|m)}
311 * @since 2.1
312 */
313 public T arcsd(final T x) {
314 // p = d, q = n, r = c, see DLMF 19.25.29 for evaluating Δ(q, p) and Δ(r, p)
315 return arcsp(x, getM(), getM().subtract(1));
316 }
317
318 /** Evaluate inverse of Jacobi elliptic function sd.
319 * @param x value of Jacobi elliptic function {@code sd(u|m)}
320 * @return u such that {@code x=sd(u|m)}
321 * @since 2.1
322 */
323 public T arcsd(final double x) {
324 return arcsd(getM().getField().getZero().newInstance(x));
325 }
326
327 /** Evaluate inverse of Jacobi elliptic function cd.
328 * @param x value of Jacobi elliptic function {@code cd(u|m)}
329 * @return u such that {@code x=cd(u|m)}
330 * @since 2.1
331 */
332 public T arccd(final T x) {
333 // p = c, q = d, r = n, see DLMF 19.25.29 for evaluating Δ(q, p) and Δ(r, q)
334 return arcpq(x, getM().subtract(1).negate(), getM());
335 }
336
337 /** Evaluate inverse of Jacobi elliptic function cd.
338 * @param x value of Jacobi elliptic function {@code cd(u|m)}
339 * @return u such that {@code x=cd(u|m)}
340 * @since 2.1
341 */
342 public T arccd(final double x) {
343 return arccd(getM().getField().getZero().newInstance(x));
344 }
345
346 /** Evaluate inverse of Jacobi elliptic function ps.
347 * <p>
348 * Here p, q, r are any permutation of the letters c, d, n.
349 * </p>
350 * @param x value of Jacobi elliptic function {@code ps(u|m)}
351 * @param deltaQP Δ(q, p) = qs²(u|m) - ps²(u|m) (equation 19.5.28 of DLMF)
352 * @param deltaRP Δ(r, p) = rs²(u|m) - ps²(u|m) (equation 19.5.28 of DLMF)
353 * @return u such that {@code x=ps(u|m)}
354 * @since 2.1
355 */
356 private T arcps(final T x, final T deltaQP, final T deltaRP) {
357 // see equation 19.25.32 in Digital Library of Mathematical Functions
358 // https://dlmf.nist.gov/19.25.E32
359 final T x2 = x.square();
360 final T rf = CarlsonEllipticIntegral.rF(x2, x2.add(deltaQP), x2.add(deltaRP));
361 return FastMath.copySign(1.0, rf.getReal()) * FastMath.copySign(1.0, x.getReal()) < 0 ?
362 rf.negate() : rf;
363 }
364
365 /** Evaluate inverse of Jacobi elliptic function sp.
366 * <p>
367 * Here p, q, r are any permutation of the letters c, d, n.
368 * </p>
369 * @param x value of Jacobi elliptic function {@code sp(u|m)}
370 * @param deltaQP Δ(q, p) = qs²(u|m) - ps²(u|m) (equation 19.5.28 of DLMF)
371 * @param deltaRP Δ(r, p) = rs²(u|m) - ps²(u|m) (equation 19.5.28 of DLMF)
372 * @return u such that {@code x=sp(u|m)}
373 * @since 2.1
374 */
375 private T arcsp(final T x, final T deltaQP, final T deltaRP) {
376 // see equation 19.25.33 in Digital Library of Mathematical Functions
377 // https://dlmf.nist.gov/19.25.E33
378 final T x2 = x.square();
379 return x.multiply(CarlsonEllipticIntegral.rF(x.getField().getOne(),
380 deltaQP.multiply(x2).add(1),
381 deltaRP.multiply(x2).add(1)));
382 }
383
384 /** Evaluate inverse of Jacobi elliptic function pq.
385 * <p>
386 * Here p, q, r are any permutation of the letters c, d, n.
387 * </p>
388 * @param x value of Jacobi elliptic function {@code pq(u|m)}
389 * @param deltaQP Δ(q, p) = qs²(u|m) - ps²(u|m) (equation 19.5.28 of DLMF)
390 * @param deltaRQ Δ(r, q) = rs²(u|m) - qs²(u|m) (equation 19.5.28 of DLMF)
391 * @return u such that {@code x=pq(u|m)}
392 * @since 2.1
393 */
394 private T arcpq(final T x, final T deltaQP, final T deltaRQ) {
395 // see equation 19.25.34 in Digital Library of Mathematical Functions
396 // https://dlmf.nist.gov/19.25.E34
397 final T x2 = x.square();
398 final T w = x2.subtract(1).negate().divide(deltaQP);
399 final T rf = CarlsonEllipticIntegral.rF(x2, x.getField().getOne(), deltaRQ.multiply(w).add(1));
400 final T positive = w.sqrt().multiply(rf);
401 return x.getReal() < 0 ? LegendreEllipticIntegral.bigK(getM()).multiply(2).subtract(positive) : positive;
402 }
403
404 /** Evaluate inverse of Jacobi elliptic function pq.
405 * <p>
406 * Here p, q, r are any permutation of the letters c, d, n.
407 * </p>
408 * <p>
409 * This computed the same thing as {@link #arcpq(CalculusFieldElement, CalculusFieldElement, CalculusFieldElement)}
410 * but uses the homogeneity property Rf(x, y, z) = Rf(ax, ay, az) / √a to get rid of the division
411 * by deltaRQ. This division induces problems in the complex case as it may lose the sign
412 * of zero for values exactly along the real or imaginary axis, hence perturbing branch cuts.
413 * </p>
414 * @param x value of Jacobi elliptic function {@code pq(u|m)}
415 * @param deltaQP Δ(q, p) = qs²(u|m) - ps²(u|m) (equation 19.5.28 of DLMF)
416 * @param deltaRQ Δ(r, q) = rs²(u|m) - qs²(u|m) (equation 19.5.28 of DLMF)
417 * @return u such that {@code x=pq(u|m)}
418 * @since 2.1
419 */
420 private T arcpqNoDivision(final T x, final T deltaQP, final T deltaRQ) {
421 // see equation 19.25.34 in Digital Library of Mathematical Functions
422 // https://dlmf.nist.gov/19.25.E34
423 final T x2 = x.square();
424 final T wDeltaQP = x2.subtract(1).negate();
425 final T rf = CarlsonEllipticIntegral.rF(x2.multiply(deltaQP), deltaQP, deltaRQ.multiply(wDeltaQP).add(deltaQP));
426 final T positive = wDeltaQP.sqrt().multiply(rf);
427 return FastMath.copySign(1.0, x.getReal()) < 0 ?
428 LegendreEllipticIntegral.bigK(getM()).multiply(2).subtract(positive) :
429 positive;
430 }
431
432 }