1   /* Copyright 2002-2022 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.forces.gravity.potential;
18  
19  
20  import java.lang.reflect.Field;
21  
22  import org.hipparchus.util.FastMath;
23  import org.hipparchus.util.MathUtils;
24  import org.junit.Assert;
25  import org.junit.Test;
26  import org.orekit.Utils;
27  import org.orekit.data.DataContext;
28  import org.orekit.forces.gravity.potential.RawSphericalHarmonicsProvider.RawSphericalHarmonics;
29  import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics;
30  import org.orekit.time.AbsoluteDate;
31  import org.orekit.time.TimeScale;
32  import org.orekit.utils.TimeSpanMap;
33  
34  @Deprecated
35  public class PulsatingSphericalHarmonicsTest {
36  
37      @Deprecated
38      @Test
39      public void testDeprecated() {
40          Utils.setDataRoot("potential");
41          TimeScale tt = DataContext.getDefault().getTimeScales().getTT();
42          GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("EIGEN-6S4-v2-truncated", true, tt));
43          UnnormalizedSphericalHarmonicsProvider ush      = GravityFieldFactory.getUnnormalizedProvider(1, 1);
44          UnnormalizedSphericalHarmonics         parsed   = ush.onDate(AbsoluteDate.J2000_EPOCH);
45          PulsatingSphericalHarmonics            psh      = toPulsatingSphericalHarmonics(ush, AbsoluteDate.J2000_EPOCH);
46          RawSphericalHarmonics                  rebuilt  = psh.onDate(AbsoluteDate.J2000_EPOCH);
47          Assert.assertEquals(AbsoluteDate.J2000_EPOCH, rebuilt.getDate());
48          Assert.assertEquals(new AbsoluteDate(1950, 1, 1, 0, 0, 0.0, tt), psh.getReferenceDate());
49          Assert.assertEquals(ush.getMu(), psh.getMu(), 1.0e-20);
50          Assert.assertEquals(ush.getAe(), psh.getAe(), 1.0e-20);
51          Assert.assertEquals(ush.getMaxDegree(),   psh.getMaxDegree());
52          Assert.assertEquals(ush.getMaxOrder(),    psh.getMaxOrder());
53          Assert.assertEquals(ush.getTideSystem(),  psh.getTideSystem());
54          for (int n = 0; n <= ush.getMaxDegree(); ++n) {
55              for (int m = 0; m <= FastMath.min(n, ush.getMaxOrder()); ++m) {
56                  Assert.assertEquals(parsed.getUnnormalizedCnm(n, m), rebuilt.getRawCnm(n, m), 1.0e-15);
57                  Assert.assertEquals(parsed.getUnnormalizedSnm(n, m), rebuilt.getRawSnm(n, m), 1.0e-15);
58              }
59          }
60      }
61  
62      private PulsatingSphericalHarmonics toPulsatingSphericalHarmonics(final UnnormalizedSphericalHarmonicsProvider ush,
63                                                                        final AbsoluteDate date) {
64          PiecewiseSphericalHarmonics     psh        = extractField(ush,            "rawProvider");
65          ConstantSphericalHarmonics      constant   = extractField(psh,            "constant");
66          AbsoluteDate[]                  references = extractField(psh,            "references");
67          double[]                        pulsations = extractField(psh,            "pulsations");
68          TimeSpanMap<PiecewisePart>      pieces     = extractField(psh,            "pieces");
69          TimeSpanMap.Span<PiecewisePart> span       = pieces.getSpan(date);
70          Flattener                       flattener  = extractField(span.getData(), "flattener");
71          TimeDependentHarmonic[]         components = extractField(span.getData(), "components");
72  
73          // extract constant part
74          double[] rawC = ((double[]) extractField(constant, "rawC")).clone();
75          double[] rawS = ((double[]) extractField(constant, "rawS")).clone();
76  
77          // extract secular part (patching constant part as required)
78          int trendReferenceIndex = -1;
79          double[] cTrend = new double[flattener.arraySize()];
80          double[] sTrend = new double[flattener.arraySize()];
81          for (int n = 0; n <= flattener.getDegree(); ++n) {
82              for (int m = 0; m <= FastMath.min(n, flattener.getOrder()); ++m) {
83                  final int index = flattener.index(n, m);
84                  if (components[index] != null) {
85                      trendReferenceIndex = ((Integer) extractField(components[index], "trendReferenceIndex")).intValue();
86                      rawC[index]  += ((Double) extractField(components[index], "cBase")).doubleValue();
87                      rawS[index]  += ((Double) extractField(components[index], "sBase")).doubleValue();
88                      cTrend[index] = ((Double) extractField(components[index], "cTrend")).doubleValue();
89                      sTrend[index] = ((Double) extractField(components[index], "sTrend")).doubleValue();
90                  }
91              }
92          }
93          RawSphericalHarmonicsProvider raw =
94                          new ConstantSphericalHarmonics(constant.getAe(), constant.getMu(), constant.getTideSystem(),
95                                                         flattener, rawC, rawS);
96          raw = new SecularTrendSphericalHarmonics(raw, references[trendReferenceIndex],
97                                                   flattener, cTrend, sTrend);
98  
99          // extract harmonic parts
100         for (int i = 0; i < pulsations.length; ++i) {
101             double[][] cosC = new double[flattener.getDegree() + 1][];
102             double[][] sinC = new double[flattener.getDegree() + 1][];
103             double[][] cosS = new double[flattener.getDegree() + 1][];
104             double[][] sinS = new double[flattener.getDegree() + 1][];
105             for (int n = 0; n <= flattener.getDegree(); ++n) {
106                 cosC[n] = new double[FastMath.min(n, flattener.getOrder()) + 1];
107                 sinC[n] = new double[FastMath.min(n, flattener.getOrder()) + 1];
108                 cosS[n] = new double[FastMath.min(n, flattener.getOrder()) + 1];
109                 sinS[n] = new double[FastMath.min(n, flattener.getOrder()) + 1];
110                 for (int m = 0; m <= FastMath.min(n, flattener.getOrder()); ++m) {
111                     final int index = flattener.index(n, m);
112                     if (components[index] != null) {
113                         cosC[n][m] = ((double[]) extractField(components[index], "cosC"))[i];
114                         sinC[n][m] = ((double[]) extractField(components[index], "sinC"))[i];
115                         cosS[n][m] = ((double[]) extractField(components[index], "cosS"))[i];
116                         sinS[n][m] = ((double[]) extractField(components[index], "sinS"))[i];
117                     }
118                 }
119             }
120             raw = new PulsatingSphericalHarmonics(raw, MathUtils.TWO_PI / pulsations[i],
121                                                   cosC, sinC, cosS, sinS);
122         }
123 
124         return (PulsatingSphericalHarmonics) raw;
125 
126     }
127 
128     @SuppressWarnings("unchecked")
129     private <T> T extractField(final Object o, final String fieldName) {
130         try {
131             Field field = o.getClass().getDeclaredField(fieldName);
132             field.setAccessible(true);
133             return (T) field.get(o);
134         } catch (NoSuchFieldException | IllegalArgumentException | IllegalAccessException e) {
135             Assert.fail(e.getLocalizedMessage());
136             return null;
137         }
138     }
139 
140 }