UpperBounds.java

  1. /* Copyright 2002-2025 CS GROUP
  2.  * Licensed to CS GROUP (CS) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * CS 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.  *   http://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.orekit.propagation.semianalytical.dsst.utilities;

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


  23. /** Utility class to compute upper bounds for truncation algorithms.
  24.  *
  25.  *  @author Pascal Parraud
  26.  */
  27. public class UpperBounds {

  29.     /** Private constructor as the class is a utility class. */
  30.     private UpperBounds() {
  31.     }

  33.     /** Get the upper bound value D<sub>n</sub><sup>l</sup>(&Chi;).
  34.     *
  35.     * @param xx value of &Chi;²
  36.     * @param xpl value of &Chi; * (&Chi;² / 2)<sup>l</sup>
  37.     * @param n index n (power of a/R)
  38.     * @param l index l (power of eccentricity)
  39.     * @return the upper bound D<sub>n</sub><sup>l</sup>(&Chi;)
  40.     */
  41.     public static double getDnl(final double xx, final double xpl, final int n, final int l) {
  42.         final int lp2 = l + 2;
  43.         if (n > lp2) {
  44.             final int ll = l * l;
  45.             double dM = xpl;
  46.             double dL = dM;
  47.             double dB = (l + 1) * xx * xpl;
  48.             for (int j = l + 3; j <= n; j++) {
  49.                 final int jm1 = j - 1;
  50.                 dL = dM;
  51.                 dM = dB;
  52.                 dB = jm1 * xx * ((2 * j - 3) * dM - (j - 2) * dL) / (jm1 * jm1 - ll);
  53.             }
  54.             return dB;
  55.         } else if (n == lp2) {
  56.             return  (l + 1) * xx * xpl;
  57.         } else {
  58.             return xpl;
  59.         }
  60.     }

  62.     /** Get the upper bound value D<sub>n</sub><sup>l</sup>(&Chi;).
  63.      *
  64.      * @param xx value of &Chi;²
  65.      * @param xpl value of &Chi; * (&Chi;² / 2)<sup>l</sup>
  66.      * @param n index n (power of a/R)
  67.      * @param l index l (power of eccentricity)
  68.      * @param <T> the type of the field elements
  69.      * @return the upper bound D<sub>n</sub><sup>l</sup>(&Chi;)
  70.      */
  71.     public static <T extends CalculusFieldElement<T>> T getDnl(final T xx, final T xpl, final int n, final int l) {
  72.         final int lp2 = l + 2;
  73.         if (n > lp2) {
  74.             final int ll = l * l;
  75.             T dM = xpl;
  76.             T dL = dM;
  77.             T dB = xx.multiply(xpl).multiply(l + 1);
  78.             for (int j = l + 3; j <= n; j++) {
  79.                 final int jm1 = j - 1;
  80.                 dL = dM;
  81.                 dM = dB;
  82.                 dB = xx.multiply(jm1).multiply(dM.multiply(2 * j - 3).subtract(dL.multiply(j - 2))).divide(jm1 * jm1 - ll);
  83.             }
  84.             return dB;
  85.         } else if (n == lp2) {
  86.             return  xx.multiply(xpl).multiply(l + 1);
  87.         } else {
  88.             return xpl;
  89.         }
  90.     }

  92.     /** Get the upper bound value R<sup>ε</sup><sub>n,m,l</sub>(γ).
  93.      *
  94.      * @param gamma value of γ
  95.      * @param n index n
  96.      * @param l index l
  97.      * @param m index m
  98.      * @param eps ε value (+1/-1)
  99.      * @param irf retrograde factor I (+1/-1)
  100.      * @return the upper bound R<sup>ε</sup><sub>n,m,l</sub>(γ)
  101.      */
  102.     public static double getRnml(final double gamma,
  103.                                  final int n, final int l, final int m,
  104.                                  final int eps, final int irf) {
  105.         // Initialization
  106.         final int mei = m * eps * irf;
  107.         final double sinisq = 1. - gamma * gamma;
  108.         // Set a lower bound for inclination
  109.         final double sininc = FastMath.max(0.03, FastMath.sqrt(sinisq));
  110.         final double onepig = 1. + gamma * irf;
  111.         final double sinincPowLmMEI = FastMath.pow(sininc, l - mei);
  112.         final double onepigPowLmMEI = FastMath.pow(onepig, mei);

  114.         // Bound for index 0
  115.         double rBound = sinincPowLmMEI * onepigPowLmMEI;

  117.         // If index > 0
  118.         if (n > l) {
  119.             final int lp1 = l + 1;

  121.             double dpnml  = lp1 * eps;
  122.             double pnml   = dpnml * gamma - m;

  124.             // If index > 1
  125.             if (n > l + 1) {
  126.                 final int ll  = l * l;
  127.                 final int ml  = m * l;
  128.                 final int mm  = m * m;

  130.                 double pn1ml  = 1.;
  131.                 double dpn1ml = 0.;
  132.                 double pn2ml  = 1.;
  133.                 double dpn2ml = 0.;
  134.                 for (int in = l + 2; in <= n; in++) {
  135.                     final int nm1   = in - 1;
  136.                     final int tnm1  = in + nm1;
  137.                     final int nnlnl = nm1 * (in * in - ll);
  138.                     final int nnmnm = in * (nm1 * nm1 - mm);
  139.                     final int c2nne = tnm1 * in * nm1 * eps;
  140.                     final int c2nml = tnm1 * ml;
  141.                     final double coef = c2nne * gamma - c2nml;

  143.                     pn2ml  = pn1ml;
  144.                     dpn2ml = dpn1ml;
  145.                     pn1ml  = pnml;
  146.                     dpn1ml = dpnml;
  147.                     pnml   = (coef * pn1ml  - nnmnm * pn2ml) / nnlnl;
  148.                     dpnml  = (coef * dpn1ml - nnmnm * dpn2ml + c2nne * pn1ml) / nnlnl;
  149.                 }
  150.             }
  151.             // Bound for index > 0
  152.             rBound *= FastMath.sqrt(pnml * pnml + dpnml * dpnml * sinisq / ((n - l) * (n + lp1)));
  153.         }

  155.         return rBound;
  156.     }

  158.     /** Get the upper bound value R<sup>ε</sup><sub>n,m,l</sub>(γ).
  159.     *
  160.     * @param gamma value of γ
  161.     * @param n index n
  162.     * @param l index l
  163.     * @param m index m
  164.     * @param eps ε value (+1/-1)
  165.     * @param irf retrograde factor I (+1/-1)
  166.     * @param <T> the type of the field elements
  167.     * @return the upper bound R<sup>ε</sup><sub>n,m,l</sub>(γ)
  168.     */
  169.     public static <T extends CalculusFieldElement<T>> T getRnml(final T gamma,
  170.                                 final int n, final int l, final int m,
  171.                                 final int eps, final int irf) {
  172.         final T zero = gamma.getField().getZero();
  173.         // Initialization
  174.         final int mei = m * eps * irf;
  175.         final T sinisq = gamma.square().negate().add(1.);
  176.         // Set a lower bound for inclination
  177.         final T sininc = FastMath.max(zero.newInstance(0.03), FastMath.sqrt(sinisq));
  178.         final T onepig = gamma.multiply(irf).add(1.);
  179.         final T sinincPowLmMEI = FastMath.pow(sininc, l - mei);
  180.         final T onepigPowLmMEI = FastMath.pow(onepig, mei);

  182.         // Bound for index 0
  183.         T rBound = sinincPowLmMEI.multiply(onepigPowLmMEI);

  185.         // If index > 0
  186.         if (n > l) {
  187.             final int lp1 = l + 1;

  189.             T dpnml  = zero.newInstance(lp1 * eps);
  190.             T pnml   = gamma.multiply(dpnml).subtract(m);

  192.             // If index > 1
  193.             if (n > l + 1) {
  194.                 final int ll  = l * l;
  195.                 final int ml  = m * l;
  196.                 final int mm  = m * m;

  198.                 T pn1ml  = zero.newInstance(1.);
  199.                 T dpn1ml = zero;
  200.                 T pn2ml  = zero.newInstance(1.);
  201.                 T dpn2ml = zero;
  202.                 for (int in = l + 2; in <= n; in++) {
  203.                     final int nm1   = in - 1;
  204.                     final int tnm1  = in + nm1;
  205.                     final int nnlnl = nm1 * (in * in - ll);
  206.                     final int nnmnm = in * (nm1 * nm1 - mm);
  207.                     final int c2nne = tnm1 * in * nm1 * eps;
  208.                     final int c2nml = tnm1 * ml;
  209.                     final T coef = gamma.multiply(c2nne).subtract(c2nml);

  211.                     pn2ml  = pn1ml;
  212.                     dpn2ml = dpn1ml;
  213.                     pn1ml  = pnml;
  214.                     dpn1ml = dpnml;
  215.                     pnml   = (coef.multiply(pn1ml).subtract(pn2ml.multiply(nnmnm))).divide(nnlnl);
  216.                     dpnml  = (coef.multiply(dpn1ml).subtract(dpn2ml.multiply(nnmnm)).add(pn1ml.multiply(c2nne))).divide(nnlnl);
  217.                 }
  218.             }
  219.             // Bound for index > 0
  220.             rBound = rBound.multiply(FastMath.sqrt(pnml.multiply(pnml).add(dpnml.multiply(dpnml).multiply(sinisq).divide((n - l) * (n + lp1)))));
  221.         }

  223.         return rBound;
  224.     }

  226. }
Created with JaCoCo 0.8.12.202403310830