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.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 }