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.forces;
18  
19  import org.hipparchus.util.FastMath;
20  import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
21  import org.orekit.propagation.semianalytical.dsst.utilities.UpperBounds;
22  
23  /**
24   * This class is a container for the common parameters used in
25   * {@link DSSTThirdBody}.
26   * <p>
27   * It performs parameters initialization at each integration step for the third
28   * body attraction perturbation. These parameters are initialize as soon as
29   * possible. In fact, they are initialized once with short period terms and
30   * don't evolve during propagation.
31   * </p>
32   * @author Bryan Cazabonne
33   * @since 11.3.3
34   */
35  public class DSSTThirdBodyStaticContext extends ForceModelContext {
36  
37      /** Max power for a/R3 in the serie expansion. */
38      private int maxAR3Pow;
39  
40      /** Max power for e in the serie expansion. */
41      private int maxEccPow;
42  
43      /** Max frequency of F. */
44      private int maxFreqF;
45  
46      /**
47       * Constructor.
48       *
49       * @param aux auxiliary elements
50       * @param x DSST Chi element
51       * @param r3 distance from center of mass of the central body to the 3rd body
52       * @param parameters force model parameters
53       */
54      public DSSTThirdBodyStaticContext(final AuxiliaryElements aux,
55                                        final double x, final double r3,
56                                        final double[] parameters) {
57          super(aux);
58  
59          // Factorials computation
60          final int dim = 2 * DSSTThirdBody.MAX_POWER;
61          final double[] fact = new double[dim];
62          fact[0] = 1.;
63          for (int i = 1; i < dim; i++) {
64              fact[i] = i * fact[i - 1];
65          }
66  
67          // Truncation tolerance.
68          final double aor = aux.getSma() / r3;
69          final double tol = (aor > .3 || aor > .15 && aux.getEcc() > .25) ? DSSTThirdBody.BIG_TRUNCATION_TOLERANCE : DSSTThirdBody.SMALL_TRUNCATION_TOLERANCE;
70  
71          // Utilities for truncation
72          // Set a lower bound for eccentricity
73          final double eo2 = FastMath.max(0.0025, 0.5 * aux.getEcc());
74          final double x2o2 = 0.5 * x * x;
75          final double[] eccPwr = new double[DSSTThirdBody.MAX_POWER];
76          final double[] chiPwr = new double[DSSTThirdBody.MAX_POWER];
77          eccPwr[0] = 1.;
78          chiPwr[0] = x;
79          for (int i = 1; i < DSSTThirdBody.MAX_POWER; i++) {
80              eccPwr[i] = eccPwr[i - 1] * eo2;
81              chiPwr[i] = chiPwr[i - 1] * x2o2;
82          }
83  
84          // Auxiliary quantities.
85          final double ao2rxx = aor / (2. * x * x);
86          double xmuarn = ao2rxx * ao2rxx * parameters[0] / (x * r3);
87          double term = 0.;
88  
89          // Compute max power for a/R3 and e.
90          maxAR3Pow = 2;
91          maxEccPow = 0;
92          int n = 2;
93          int m = 2;
94          int nsmd2 = 0;
95  
96          do {
97              // Upper bound for Tnm.
98              term =
99                  xmuarn *
100                    (fact[n + m] / (fact[nsmd2] * fact[nsmd2 + m])) *
101                    (fact[n + m + 1] / (fact[m] * fact[n + 1])) *
102                    (fact[n - m + 1] / fact[n + 1]) * eccPwr[m] *
103                    UpperBounds.getDnl(x * x, chiPwr[m], n + 2, m);
104 
105             if (term < tol) {
106                 if (m == 0) {
107                     break;
108                 } else if (m < 2) {
109                     xmuarn *= ao2rxx;
110                     m = 0;
111                     n++;
112                     nsmd2++;
113                 } else {
114                     m -= 2;
115                     nsmd2++;
116                 }
117             } else {
118                 maxAR3Pow = n;
119                 maxEccPow = FastMath.max(m, maxEccPow);
120                 xmuarn *= ao2rxx;
121                 m++;
122                 n++;
123             }
124         } while (n < DSSTThirdBody.MAX_POWER);
125 
126         maxEccPow = FastMath.min(maxAR3Pow, maxEccPow);
127         maxFreqF = maxAR3Pow + 1;
128 
129     }
130 
131     /**
132      * Get the value of max power for a/R3 in the serie expansion.
133      *
134      * @return maxAR3Pow
135      */
136     public int getMaxAR3Pow() {
137         return maxAR3Pow;
138     }
139 
140     /**
141      * Get the value of max power for e in the serie expansion.
142      *
143      * @return maxEccPow
144      */
145     public int getMaxEccPow() {
146         return maxEccPow;
147     }
148 
149     /**
150      * Get the value of max frequency of F.
151      *
152      * @return maxFreqF
153      */
154     public int getMaxFreqF() {
155         return maxFreqF;
156     }
157 
158 }