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;
18  
19  import org.hipparchus.CalculusFieldElement;
20  import org.hipparchus.Field;
21  import org.hipparchus.util.CombinatoricsUtils;
22  import org.hipparchus.util.Binary64Field;
23  import org.hipparchus.util.FastMath;
24  import org.junit.jupiter.api.Assertions;
25  import org.junit.jupiter.api.BeforeEach;
26  import org.junit.jupiter.api.Test;
27  
28  import java.lang.reflect.InvocationTargetException;
29  import java.lang.reflect.Method;
30  
31  public class FieldGammaMnsFunctionTest {
32  
33      int      nMax;
34  
35      @Test
36      public void testIndex()
37          throws NoSuchMethodException, SecurityException,
38                 IllegalAccessException, IllegalArgumentException, InvocationTargetException {
39          Method indexM = FieldGammaMnsFunction.class.getDeclaredMethod("index",
40                                                                   Integer.TYPE, Integer.TYPE, Integer.TYPE);
41          indexM.setAccessible(true);
42          int i = 0;
43          for (int n = 0; n <= nMax; ++n) {
44              for (int m = 0; m <= n; ++m) {
45                  for (int s = -n; s <= n; ++s) {
46                      Assertions.assertEquals(i++, indexM.invoke(null, m, n, s));
47                  }
48              }
49          }
50      }
51  
52      @Test
53      public void testPrecomputedRatios()
54          throws NoSuchFieldException, SecurityException, IllegalArgumentException, IllegalAccessException {
55          doTestPrecomputedRatios(Binary64Field.getInstance());
56      }
57  
58      private <T extends CalculusFieldElement<T>> void doTestPrecomputedRatios(Field<T> field)
59          throws NoSuchFieldException, SecurityException, IllegalArgumentException, IllegalAccessException {
60          final T zero = field.getZero();
61          java.lang.reflect.Field precomputedF = FieldGammaMnsFunction.class.getDeclaredField("PRECOMPUTED_RATIOS");
62          precomputedF.setAccessible(true);
63          new FieldGammaMnsFunction<>(nMax, zero.add(0.5), +1, field);
64          double[] precomputed = (double[]) precomputedF.get(null);
65          int i = 0;
66          for (int n = 0; n <= nMax; ++n) {
67              for (int m = 0; m <= n; ++m) {
68                  for (int s = -n; s <= n; ++s) {
69                      // compare against naive implementation
70                      double r = naiveRatio(m, n, s);
71                      Assertions.assertEquals(r, precomputed[i++], 2.0e-14 * r);
72                  }
73              }
74          }
75      }
76  
77      @Test
78      public void testReallocate()
79          throws NoSuchFieldException, SecurityException, IllegalArgumentException, IllegalAccessException {
80          doTestReallocate(Binary64Field.getInstance());
81      }
82  
83      private <T extends CalculusFieldElement<T>> void doTestReallocate(Field<T> field)
84          throws NoSuchFieldException, SecurityException, IllegalArgumentException, IllegalAccessException {
85          final T zero = field.getZero();
86          java.lang.reflect.Field precomputedF = FieldGammaMnsFunction.class.getDeclaredField("PRECOMPUTED_RATIOS");
87          precomputedF.setAccessible(true);
88          precomputedF.set(null, new double[0]);
89          new FieldGammaMnsFunction<>(nMax, zero.add(0.5), +1, field);
90          double[] orginalPrecomputed = (double[]) precomputedF.get(null);
91          Assertions.assertEquals((nMax + 1) * (nMax + 2) * (4 * nMax + 3) / 6, orginalPrecomputed.length);
92          new FieldGammaMnsFunction<>(nMax + 3, zero.add(0.5), +1, field);
93          double[] reallocatedPrecomputed = (double[]) precomputedF.get(null);
94          Assertions.assertEquals((nMax + 4) * (nMax + 5) * (4 * nMax + 15) / 6, reallocatedPrecomputed.length);
95          for (int i = 0; i < orginalPrecomputed.length; ++i) {
96              Assertions.assertEquals(orginalPrecomputed[i], reallocatedPrecomputed[i],
97                                  1.0e-15 * orginalPrecomputed[i]);
98          }
99      }
100 
101 
102     @Test
103     public void testValue() {
104         doTestValue(Binary64Field.getInstance());
105     }
106 
107     private <T extends CalculusFieldElement<T>> void doTestValue(Field<T> field) {
108         final T zero = field.getZero();
109         for (int bigI : new int[] { -1, +1 }) {
110             for (T gamma = zero; gamma.getReal() <= 1; gamma = gamma.add(1.0 / 64.0)) {
111                 FieldGammaMnsFunction<T> gammaMNS = new FieldGammaMnsFunction<>(nMax, gamma, bigI, field);
112                 for (int n = 0; n <= nMax; ++n) {
113                     for (int m = 0; m <= n; ++m) {
114                         for (int s = -n; s <= n; ++s) {
115                             // compare against naive implementation
116                             final T v = naiveValue(bigI, gamma, m, n, s, field);
117                             final T g = gammaMNS.getValue(m, n, s);
118                             if (Double.isInfinite(v.getReal())) {
119                                 Assertions.assertTrue(Double.isInfinite(g.getReal()));
120                                 Assertions.assertTrue(v.multiply(g).getReal() > 0);
121                             } else {
122                                 Assertions.assertEquals(v.getReal(), g.getReal(), FastMath.abs(v).multiply(2.0e-14).getReal());
123                             }
124                         }
125                     }
126                 }
127             }
128         }
129     }
130 
131     private <T extends CalculusFieldElement<T>> T naiveValue(final int bigI, final T gamma,
132                                        final int m, final int n, final int s,
133                                        final Field<T> field) {
134         final T zero = field.getZero();
135         final T f = gamma.multiply(bigI).add(1.);
136         if (s <= -m) {
137             return FastMath.pow(zero.subtract(1.), zero.add(m - s)).multiply(FastMath.pow(zero.add(2.), zero.add(s))).multiply(FastMath.pow(f, -bigI * m));
138         } else if (s < m) {
139             return FastMath.pow(zero.subtract(1.), zero.add(m - s)).multiply(FastMath.pow(zero.add(2), zero.subtract(m))).multiply(naiveRatio(m, n, s)).multiply(FastMath.pow(f, bigI * s));
140         } else {
141             return  FastMath.pow(zero.add(2.), zero.subtract(s)).multiply(FastMath.pow(f, bigI * m));
142         }
143     }
144 
145     private double naiveRatio(final int m, final int n, final int s) {
146         return (CombinatoricsUtils.factorialDouble(n + m) * CombinatoricsUtils.factorialDouble(n - m)) /
147                (CombinatoricsUtils.factorialDouble(n + s) * CombinatoricsUtils.factorialDouble(n - s));
148     }
149 
150     @BeforeEach
151     public void setUp() {
152         nMax = 12;
153     }
154 }