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.util.FastMath;
21 import org.hipparchus.util.FieldSinCos;
22
23 /** Algorithm computing Jacobi theta functions.
24 * @param <T> the type of the field elements
25 * @since 2.0
26 */
27 public class FieldJacobiTheta<T extends CalculusFieldElement<T>> {
28
29 /** Maximum number of terms in the Fourier series. */
30 private static final int N_MAX = 100;
31
32 /** Nome. */
33 private final T q;
34
35 /** q². */
36 private final T qSquare;
37
38 /** ∜q. */
39 private final T qFourth;
40
41 /** Simple constructor.
42 * <p>
43 * The nome {@code q} can be computed using ratios of complete elliptic integrals
44 * ({@link org.hipparchus.special.elliptic.legendre.LegendreEllipticIntegral#nome(CalculusFieldElement)
45 * LegendreEllipticIntegral.nome(m)} which are themselves defined in term of parameter m,
46 * where m=k² and k is the elliptic modulus.
47 * </p>
48 * @param q nome
49 */
50 public FieldJacobiTheta(final T q) {
51 this.q = q;
52 this.qSquare = q.multiply(q);
53 this.qFourth = FastMath.sqrt(FastMath.sqrt(q));
54 }
55
56 /** Get the nome.
57 * @return nome
58 */
59 public T getQ() {
60 return q;
61 }
62
63 /** Evaluate the Jacobi theta functions.
64 * @param z argument of the functions
65 * @return container for the four Jacobi theta functions θ₁(z|τ), θ₂(z|τ), θ₃(z|τ), and θ₄(z|τ)
66 */
67 public FieldTheta<T> values(final T z) {
68
69 // the computation is based on Fourier series,
70 // see Digital Library of Mathematical Functions section 20.2
71 // https://dlmf.nist.gov/20.2
72 final T zero = q.getField().getZero();
73 final T one = q.getField().getOne();
74
75 // base angle for Fourier Series
76 final FieldSinCos<T> sc1 = FastMath.sinCos(z);
77
78 // recursion rules initialization
79 double sgn = 1.0;
80 T qNN = one;
81 T qTwoN = one;
82 T qNNp1 = one;
83 FieldSinCos<T> sc2n1 = sc1;
84 final double eps = FastMath.ulp(one).getReal();
85
86 // Fourier series
87 T sum1 = sc1.sin();
88 T sum2 = sc1.cos();
89 T sum3 = zero;
90 T sum4 = zero;
91 for (int n = 1; n < N_MAX; ++n) {
92
93 sgn = -sgn; // (-1)ⁿ⁻¹ ← (-1)ⁿ
94 qNN = qNN.multiply(qTwoN).multiply(q); // q⁽ⁿ⁻¹⁾⁽ⁿ⁻¹⁾ ← qⁿⁿ
95 qTwoN = qTwoN.multiply(qSquare); // q²⁽ⁿ⁻¹⁾ ← q²ⁿ
96 qNNp1 = qNNp1.multiply(qTwoN); // q⁽ⁿ⁻¹⁾ⁿ ← qⁿ⁽ⁿ⁺¹⁾
97
98 sc2n1 = FieldSinCos.sum(sc2n1, sc1); // {sin|cos}([2n-1] z) ← {sin|cos}(2n z)
99 sum3 = sum3.add(sc2n1.cos().multiply(qNN));
100 sum4 = sum4.add(sc2n1.cos().multiply(qNN.multiply(sgn)));
101
102 sc2n1 = FieldSinCos.sum(sc2n1, sc1); // {sin|cos}(2n z) ← {sin|cos}([2n+1] z)
103 sum1 = sum1.add(sc2n1.sin().multiply(qNNp1.multiply(sgn)));
104 sum2 = sum2.add(sc2n1.cos().multiply(qNNp1));
105
106 if (qNNp1.norm() <= eps) {
107 // we have reach convergence
108 break;
109 }
110
111 }
112
113 return new FieldTheta<>(sum1.multiply(qFourth.multiply(2)),
114 sum2.multiply(qFourth.multiply(2)),
115 sum3.multiply(2).add(1),
116 sum4.multiply(2).add(1));
117
118 }
119
120 }