RdRealDuplication.java

  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;

  19. import org.hipparchus.util.FastMath;
  20. import org.hipparchus.util.MathArrays;

  22. /** Duplication algorithm for Carlson R<sub>D</sub> elliptic integral.
  23.  * @since 2.0
  24.  */
  25. class RdRealDuplication extends RealDuplication {

  27.     /** Constant term in R<sub>J</sub> and R<sub>D</sub> polynomials. */
  28.     static final double CONSTANT = 4084080;

  30.     /** Coefficient of E₂ in R<sub>J</sub> and R<sub>D</sub> polynomials. */
  31.     static final double E2 = -875160;

  33.     /** Coefficient of E₃ in R<sub>J</sub> and R<sub>D</sub> polynomials. */
  34.     static final double E3 = 680680;

  36.     /** Coefficient of E₂² in R<sub>J</sub> and R<sub>D</sub> polynomials. */
  37.     static final double E2_E2 = 417690;

  39.     /** Coefficient of E₄ in R<sub>J</sub> and R<sub>D</sub> polynomials. */
  40.     static final double E4 = -556920;

  42.     /** Coefficient of E₂E₃ in R<sub>J</sub> and R<sub>D</sub> polynomials. */
  43.     static final double E2_E3 = -706860;

  45.     /** Coefficient of E₅ in R<sub>J</sub> and R<sub>D</sub> polynomials. */
  46.     static final double E5 = 471240;

  48.     /** Coefficient of E₂³ in R<sub>J</sub> and R<sub>D</sub> polynomials. */
  49.     static final double E2_E2_E2 = -255255;

  51.     /** Coefficient of E₃² in R<sub>J</sub> and R<sub>D</sub> polynomials. */
  52.     static final double E3_E3 = 306306;

  54.     /** Coefficient of E₂E₄ in R<sub>J</sub> and R<sub>D</sub> polynomials. */
  55.     static final double E2_E4 = 612612;

  57.     /** Coefficient of E₂²E₃ in R<sub>J</sub> and R<sub>D</sub> polynomials. */
  58.     static final double E2_E2_E3 = 675675;

  60.     /** Coefficient of E₃E₄+E₂E₅ in R<sub>J</sub> and R<sub>D</sub> polynomials. */
  61.     static final double E3_E4_P_E2_E5 = -540540;

  63.     /** Denominator in R<sub>J</sub> and R<sub>D</sub> polynomials. */
  64.     static final double DENOMINATOR = 4084080;

  66.     /** Partial sum. */
  67.     private double sum;

  69.     /** Simple constructor.
  70.      * @param x first symmetric variable of the integral
  71.      * @param y second symmetric variable of the integral
  72.      * @param z third symmetric variable of the integral
  73.      */
  74.     RdRealDuplication(final double x, final double y, final double z) {
  75.         super(x, y, z);
  76.         sum = 0;
  77.     }

  79.     /** {@inheritDoc} */
  80.     @Override
  81.     protected void initialMeanPoint(final double[] va) {
  82.         va[3] = (va[0] + va[1] + va[2] * 3.0) / 5.0;
  83.     }

  85.     /** {@inheritDoc} */
  86.     @Override
  87.     protected double convergenceCriterion(final double r, final double max) {
  88.         return max / (FastMath.sqrt(FastMath.sqrt(FastMath.sqrt(r * 0.25))));
  89.     }

  91.     /** {@inheritDoc} */
  92.     @Override
  93.     protected void update(final int m, final double[] vaM, final double[] sqrtM, final  double fourM) {

  95.         // equation 2.29 in Carlson[1995]
  96.         final double lambdaA = sqrtM[0] * sqrtM[1];
  97.         final double lambdaB = sqrtM[0] * sqrtM[2];
  98.         final double lambdaC = sqrtM[1] * sqrtM[2];

  100.         // running sum in equation 2.34 in Carlson[1995]
  101.         final double lambda = lambdaA + lambdaB + lambdaC;
  102.         sum += 1.0 / ((vaM[2] + lambda) * sqrtM[2] * fourM);

  104.         // equations 2.29 and 2.30 in Carlson[1995]
  105.         vaM[0] = MathArrays.linearCombination(0.25, vaM[0], 0.25, lambdaA, 0.25, lambdaB, 0.25, lambdaC); // xₘ
  106.         vaM[1] = MathArrays.linearCombination(0.25, vaM[1], 0.25, lambdaA, 0.25, lambdaB, 0.25, lambdaC); // yₘ
  107.         vaM[2] = MathArrays.linearCombination(0.25, vaM[2], 0.25, lambdaA, 0.25, lambdaB, 0.25, lambdaC); // zₘ
  108.         vaM[3] = MathArrays.linearCombination(0.25, vaM[3], 0.25, lambdaA, 0.25, lambdaB, 0.25, lambdaC); // aₘ

  110.     }

  112.     /** {@inheritDoc} */
  113.     @Override
  114.     protected double evaluate(final double[] va0, final double aM, final  double fourM) {

  116.         // compute symmetric differences
  117.         final double inv   = 1.0 / (aM * fourM);
  118.         final double bigX  = (va0[3] - va0[0]) * inv;
  119.         final double bigY  = (va0[3] - va0[1]) * inv;
  120.         final double bigZ  = (bigX + bigY) / -3;
  121.         final double bigXY = bigX * bigY;
  122.         final double bigZ2 = bigZ * bigZ;

  124.         // compute elementary symmetric functions (we already know e1 = 0 by construction)
  125.         final double e2  = bigXY - bigZ2 * 6;
  126.         final double e3  = (bigXY * 3 - bigZ2 * 8) * bigZ;
  127.         final double e4  = (bigXY - bigZ2) * 3 * bigZ2;
  128.         final double e5  = bigXY * bigZ2 * bigZ;

  130.         final double e2e2   =   e2 * e2;
  131.         final double e2e3   =   e2 * e3;
  132.         final double e2e4   =   e2 * e4;
  133.         final double e2e5   =   e2 * e5;
  134.         final double e3e3   =   e3 * e3;
  135.         final double e3e4   =   e3 * e4;
  136.         final double e2e2e2 = e2e2 * e2;
  137.         final double e2e2e3 = e2e2 * e3;

  139.         // evaluate integral using equation 19.36.1 in DLMF
  140.         // (which add more terms than equation 2.7 in Carlson[1995])
  141.         final double poly = ((e3e4 + e2e5) * E3_E4_P_E2_E5 +
  142.                               e2e2e3 * E2_E2_E3 +
  143.                               e2e4 * E2_E4 +
  144.                               e3e3 * E3_E3 +
  145.                               e2e2e2 * E2_E2_E2 +
  146.                               e5 * E5 +
  147.                               e2e3 * E2_E3 +
  148.                               e4 * E4 +
  149.                               e2e2 * E2_E2 +
  150.                               e3 * E3 +
  151.                               e2 * E2 +
  152.                               CONSTANT) /
  153.                              DENOMINATOR;
  154.         final double polyTerm = poly / (aM * FastMath.sqrt(aM) * fourM);

  156.         return polyTerm + sum * 3;

  158.     }

  160. }
Created with JaCoCo 0.8.12.202403310830