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.complex.Complex;
20 import org.hipparchus.util.FastMath;
21 import org.hipparchus.util.FieldSinCos;
22 import org.hipparchus.util.Precision;
23
24 /** Algorithm computing Jacobi theta functions.
25 * @since 2.0
26 */
27 public class JacobiTheta {
28
29 /** Maximum number of terms in the Fourier series. */
30 private static final int N_MAX = 100;
31
32 /** Nome. */
33 private final double q;
34
35 /** q². */
36 private final double qSquare;
37
38 /** ∜q. */
39 private final double 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(double)
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 JacobiTheta(final double q) {
51 this.q = q;
52 this.qSquare = q * q;
53 this.qFourth = FastMath.sqrt(FastMath.sqrt(q));
54 }
55
56 /** Get the nome.
57 * @return nome
58 */
59 public double 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 Theta values(final Complex 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
73 // base angle for Fourier Series
74 final FieldSinCos<Complex> sc1 = FastMath.sinCos(z);
75
76 // recursion rules initialization
77 double sgn = 1.0;
78 double qNN = 1.0;
79 double qTwoN = 1.0;
80 double qNNp1 = 1.0;
81 FieldSinCos<Complex> sc2n1 = sc1;
82
83 // Fourier series
84 Complex sum1 = sc1.sin();
85 Complex sum2 = sc1.cos();
86 Complex sum3 = Complex.ZERO;
87 Complex sum4 = Complex.ZERO;
88 for (int n = 1; n < N_MAX; ++n) {
89
90 sgn = -sgn; // (-1)ⁿ⁻¹ ← (-1)ⁿ
91 qNN = qNN * qTwoN * q; // q⁽ⁿ⁻¹⁾⁽ⁿ⁻¹⁾ ← qⁿⁿ
92 qTwoN = qTwoN * qSquare; // q²⁽ⁿ⁻¹⁾ ← q²ⁿ
93 qNNp1 = qNNp1 * qTwoN; // q⁽ⁿ⁻¹⁾ⁿ ← qⁿ⁽ⁿ⁺¹⁾
94
95 sc2n1 = FieldSinCos.sum(sc2n1, sc1); // {sin|cos}([2n-1] z) ← {sin|cos}(2n z)
96 sum3 = sum3.add(sc2n1.cos().multiply(qNN));
97 sum4 = sum4.add(sc2n1.cos().multiply(sgn * qNN));
98
99 sc2n1 = FieldSinCos.sum(sc2n1, sc1); // {sin|cos}(2n z) ← {sin|cos}([2n+1] z)
100 sum1 = sum1.add(sc2n1.sin().multiply(sgn * qNNp1));
101 sum2 = sum2.add(sc2n1.cos().multiply(qNNp1));
102
103 if (FastMath.abs(qNNp1) <= Precision.EPSILON) {
104 // we have reach convergence
105 break;
106 }
107
108 }
109
110 return new Theta(sum1.multiply(2 * qFourth),
111 sum2.multiply(2 * qFourth),
112 sum3.multiply(2).add(1),
113 sum4.multiply(2).add(1));
114
115 }
116
117 }