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.carlson;
18  
19  import org.hipparchus.util.FastMath;
20  import org.hipparchus.util.MathArrays;
21  
22  /** Duplication algorithm for Carlson R<sub>J</sub> elliptic integral.
23   * @since 2.0
24   */
25  class RjRealDuplication extends RealDuplication {
26  
27      /** Delta product. */
28      private double delta;
29  
30      /** sₘ iteration parameter. */
31      private double sM;
32  
33      /** Simple constructor.
34       * @param x first symmetric variable of the integral
35       * @param y second symmetric variable of the integral
36       * @param z third symmetric variable of the integral
37       * @param p fourth <em>not</em> symmetric variable of the integral
38       * @param delta precomputed value of (p-x)(p-y)(p-z)
39       */
40      RjRealDuplication(final double x, final double y, final double z, final double p, final double delta) {
41          super(x, y, z, p);
42          this.delta = delta;
43      }
44  
45      /** {@inheritDoc} */
46      @Override
47      protected void initialMeanPoint(final double[] va) {
48          va[4] = (va[0] + va[1] + va[2] + va[3] * 2) / 5.0;
49      }
50  
51      /** {@inheritDoc} */
52      @Override
53      protected double convergenceCriterion(final double r, final double max) {
54          return max / (FastMath.sqrt(FastMath.sqrt(FastMath.sqrt(r * 0.25))));
55      }
56  
57      /** {@inheritDoc} */
58      @Override
59      protected void update(final int m, final double[] vaM, final double[] sqrtM, final  double fourM) {
60          final double dM = (sqrtM[3] + sqrtM[0]) * (sqrtM[3] + sqrtM[1]) * (sqrtM[3] + sqrtM[2]);
61          if (m == 0) {
62              sM = dM * 0.5;
63          } else {
64              // equation A.3 in Carlson[2000]
65              final double rM = sM * (FastMath.sqrt(delta / (sM * sM * fourM) + 1.0) + 1.0);
66              sM = (dM * rM - delta / (fourM * fourM)) / ((dM + rM / fourM) * 2);
67          }
68  
69          // equation 2.19 in Carlson[1995]
70          final double lambdaA = sqrtM[0] * sqrtM[1];
71          final double lambdaB = sqrtM[0] * sqrtM[2];
72          final double lambdaC = sqrtM[1] * sqrtM[2];
73  
74          // equations 2.19 and 2.20 in Carlson[1995]
75          vaM[0] = MathArrays.linearCombination(0.25, vaM[0], 0.25, lambdaA, 0.25, lambdaB, 0.25, lambdaC); // xₘ
76          vaM[1] = MathArrays.linearCombination(0.25, vaM[1], 0.25, lambdaA, 0.25, lambdaB, 0.25, lambdaC); // yₘ
77          vaM[2] = MathArrays.linearCombination(0.25, vaM[2], 0.25, lambdaA, 0.25, lambdaB, 0.25, lambdaC); // zₘ
78          vaM[3] = MathArrays.linearCombination(0.25, vaM[3], 0.25, lambdaA, 0.25, lambdaB, 0.25, lambdaC); // pₘ
79          vaM[4] = MathArrays.linearCombination(0.25, vaM[4], 0.25, lambdaA, 0.25, lambdaB, 0.25, lambdaC); // aₘ
80  
81      }
82  
83      /** {@inheritDoc} */
84      @Override
85      protected double evaluate(final double[] va0, final double aM, final  double fourM) {
86  
87          // compute symmetric differences
88          final double inv    = 1.0 / (aM * fourM);
89          final double bigX   = (va0[4] - va0[0]) * inv;
90          final double bigY   = (va0[4] - va0[1]) * inv;
91          final double bigZ   = (va0[4] - va0[2]) * inv;
92          final double bigP   = (bigX + bigY + bigZ) * -0.5;
93          final double bigP2  = bigP * bigP;
94  
95          // compute elementary symmetric functions (we already know e1 = 0 by construction)
96          final double xyz    = bigX * bigY * bigZ;
97          final double e2     = bigX * (bigY + bigZ) + bigY * bigZ - bigP * bigP * 3;
98          final double e3     = xyz + bigP * 2 * (e2 + bigP2 * 2);
99          final double e4     = (xyz * 2 + bigP * (e2 + bigP2 * 3)) * bigP;
100         final double e5     = xyz * bigP2;
101 
102         final double e2e2   = e2   * e2;
103         final double e2e3   = e2   * e3;
104         final double e2e4   = e2   * e4;
105         final double e2e5   = e2   * e5;
106         final double e3e3   = e3   * e3;
107         final double e3e4   = e3   * e4;
108         final double e2e2e2 = e2e2 * e2;
109         final double e2e2e3 = e2e2 * e3;
110 
111         // evaluate integral using equation 19.36.1 in DLMF
112         // (which add more terms than equation 2.7 in Carlson[1995])
113         final double poly = ((e3e4 + e2e5) * RdRealDuplication.E3_E4_P_E2_E5 +
114                               e2e2e3       * RdRealDuplication.E2_E2_E3 +
115                               e2e4         * RdRealDuplication.E2_E4 +
116                               e3e3         * RdRealDuplication.E3_E3 +
117                               e2e2e2       * RdRealDuplication.E2_E2_E2 +
118                               e5           * RdRealDuplication.E5 +
119                               e2e3         * RdRealDuplication.E2_E3 +
120                               e4           * RdRealDuplication.E4 +
121                               e2e2         * RdRealDuplication.E2_E2 +
122                               e3           * RdRealDuplication.E3 +
123                               e2           * RdRealDuplication.E2 +
124                               RdRealDuplication.CONSTANT) /
125                              RdRealDuplication.DENOMINATOR;
126         final double polyTerm = poly / (aM * FastMath.sqrt(aM) * fourM);
127 
128         // compute a single R_C term
129         final double rcTerm = new RcRealDuplication(1.0, delta / (sM * sM * fourM) + 1.0).integral() * 3 / sM;
130 
131         return polyTerm + rcTerm;
132 
133     }
134 
135 }