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.propagation.semianalytical.dsst;
18  
19  import java.io.IOException;
20  import java.lang.reflect.InvocationTargetException;
21  import java.lang.reflect.Method;
22  import java.text.ParseException;
23  import java.util.ArrayList;
24  import java.util.Arrays;
25  import java.util.List;
26  
27  import org.hipparchus.exception.LocalizedCoreFormats;
28  import org.hipparchus.util.FastMath;
29  import org.junit.Assert;
30  import org.junit.Before;
31  import org.junit.Test;
32  import org.orekit.Utils;
33  import org.orekit.bodies.CelestialBodyFactory;
34  import org.orekit.bodies.OneAxisEllipsoid;
35  import org.orekit.errors.OrekitException;
36  import org.orekit.forces.gravity.potential.GravityFieldFactory;
37  import org.orekit.forces.gravity.potential.ICGEMFormatReader;
38  import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
39  import org.orekit.frames.Frame;
40  import org.orekit.frames.FramesFactory;
41  import org.orekit.orbits.EquinoctialOrbit;
42  import org.orekit.orbits.KeplerianOrbit;
43  import org.orekit.orbits.Orbit;
44  import org.orekit.orbits.PositionAngle;
45  import org.orekit.propagation.PropagationType;
46  import org.orekit.propagation.SpacecraftState;
47  import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
48  import org.orekit.propagation.semianalytical.dsst.forces.DSSTTesseral;
49  import org.orekit.propagation.semianalytical.dsst.forces.ShortPeriodTerms;
50  import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
51  import org.orekit.time.AbsoluteDate;
52  import org.orekit.time.TimeScalesFactory;
53  import org.orekit.utils.Constants;
54  
55  public class DSSTTesseralTest {
56  
57      @Test
58      public void testGetMeanElementRate() {
59          
60          // Central Body geopotential 4x4
61          final UnnormalizedSphericalHarmonicsProvider provider =
62                  GravityFieldFactory.getUnnormalizedProvider(4, 4);
63          
64          final Frame frame = FramesFactory.getEME2000();
65          final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
66          final AbsoluteDate initDate = new AbsoluteDate(2007, 4, 16, 0, 46, 42.400, TimeScalesFactory.getUTC());
67          
68          // a  = 26559890 m
69          // ex = 2.719455286199036E-4
70          // ey = 0.0041543085910249414
71          // hx = -0.3412974060023717
72          // hy = 0.3960084733107685
73          // lM = 8.566537840341699 rad
74          final Orbit orbit = new EquinoctialOrbit(2.655989E7,
75                                                   2.719455286199036E-4,
76                                                   0.0041543085910249414,
77                                                   -0.3412974060023717,
78                                                   0.3960084733107685,
79                                                   8.566537840341699,
80                                                   PositionAngle.TRUE,
81                                                   frame,
82                                                   initDate,
83                                                   provider.getMu());
84          
85          final SpacecraftState state = new SpacecraftState(orbit, 1000.0);
86          
87          final AuxiliaryElements auxiliaryElements = new AuxiliaryElements(state.getOrbit(), 1);
88          
89          final DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
90                                                           Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider,
91                                                           4, 4, 4, 8, 4, 4, 2);
92          // Force model parameters
93          final double[] parameters = tesseral.getParameters();
94  
95          // Initialize force model
96          tesseral.initializeShortPeriodTerms(auxiliaryElements, PropagationType.MEAN, parameters);
97  
98          final double[] elements = new double[7];
99          Arrays.fill(elements, 0.0);
100         
101         final double[] daidt = tesseral.getMeanElementRate(state, auxiliaryElements, parameters);
102         for (int i = 0; i < daidt.length; i++) {
103             elements[i] = daidt[i];
104         }
105         
106         Assert.assertEquals(7.120011500375922E-5,   elements[0], 1.e-20);
107         Assert.assertEquals(-1.109767646425212E-11, elements[1], 1.e-26);
108         Assert.assertEquals(2.3036711391089307E-11, elements[2], 1.e-26);
109         Assert.assertEquals(2.499304852807308E-12,  elements[3], 1.e-27);
110         Assert.assertEquals(1.3899097178558372E-13, elements[4], 1.e-28);
111         Assert.assertEquals(5.795522421338584E-12,  elements[5], 1.e-27);
112         
113     }
114     
115     @Test
116     public void testShortPeriodTerms() throws IllegalArgumentException {
117         
118         Utils.setDataRoot("regular-data:potential/icgem-format");
119         GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("^eigen-6s-truncated$", false));
120         UnnormalizedSphericalHarmonicsProvider nshp = GravityFieldFactory.getUnnormalizedProvider(8, 8);
121         Orbit orbit = new KeplerianOrbit(13378000, 0.05, 0, 0, FastMath.PI, 0, PositionAngle.MEAN,
122                                          FramesFactory.getTOD(false),
123                                          new AbsoluteDate(2003, 5, 6, TimeScalesFactory.getUTC()),
124                                          nshp.getMu());
125         
126         OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
127                                                       Constants.WGS84_EARTH_FLATTENING,
128                                                       FramesFactory.getGTOD(false));
129         
130         // Force model
131         final DSSTForceModel force = new DSSTTesseral(earth.getBodyFrame(),
132                                                       Constants.WGS84_EARTH_ANGULAR_VELOCITY,
133                                                       nshp, 8, 8, 4, 12, 8, 8, 4);
134         
135         // Initial state
136         final SpacecraftState meanState = new SpacecraftState(orbit, 45.0);
137         
138         //Create the auxiliary object
139         final AuxiliaryElements aux = new AuxiliaryElements(orbit, 1);
140        
141         final List<ShortPeriodTerms> shortPeriodTerms = new ArrayList<ShortPeriodTerms>();
142 
143         force.registerAttitudeProvider(null);
144         shortPeriodTerms.addAll(force.initializeShortPeriodTerms(aux, PropagationType.OSCULATING, force.getParameters()));
145         force.updateShortPeriodTerms(force.getParameters(), meanState);
146         
147         double[] y = new double[6];
148         for (final ShortPeriodTerms spt : shortPeriodTerms) {
149             final double[] shortPeriodic = spt.value(meanState.getOrbit());
150             for (int i = 0; i < shortPeriodic.length; i++) {
151                 y[i] += shortPeriodic[i];
152             }
153         }
154         
155         Assert.assertEquals(5.192409957353236,      y[0], 1.e-15);
156         Assert.assertEquals(9.660364749662076E-7,   y[1], 1.e-22);
157         Assert.assertEquals(1.542008987162059E-6,   y[2], 1.e-21);
158         Assert.assertEquals(-4.9944146013126755E-8, y[3], 1.e-22);
159         Assert.assertEquals(-4.500974242661177E-8,  y[4], 1.e-22);
160         Assert.assertEquals(-2.785213556107612E-7,  y[5], 1.e-21);
161     }
162 
163     @Test
164     public void testIssue625() {
165 
166         // Central Body geopotential 4x4
167         final UnnormalizedSphericalHarmonicsProvider provider =
168                 GravityFieldFactory.getUnnormalizedProvider(4, 4);
169         
170         final Frame frame = FramesFactory.getEME2000();
171         final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
172         final AbsoluteDate initDate = new AbsoluteDate(2007, 4, 16, 0, 46, 42.400, TimeScalesFactory.getUTC());
173         
174         // a  = 26559890 m
175         // ex = 2.719455286199036E-4
176         // ey = 0.0041543085910249414
177         // hx = -0.3412974060023717
178         // hy = 0.3960084733107685
179         // lM = 8.566537840341699 rad
180         final Orbit orbit = new EquinoctialOrbit(2.655989E7,
181                                                  2.719455286199036E-4,
182                                                  0.0041543085910249414,
183                                                  -0.3412974060023717,
184                                                  0.3960084733107685,
185                                                  8.566537840341699,
186                                                  PositionAngle.TRUE,
187                                                  frame,
188                                                  initDate,
189                                                  provider.getMu());
190         
191         final SpacecraftState state = new SpacecraftState(orbit, 1000.0);
192         
193         final AuxiliaryElements auxiliaryElements = new AuxiliaryElements(state.getOrbit(), 1);
194 
195         // Tesseral force model
196         final DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
197                                                          Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider,
198                                                          4, 4, 4, 8, 4, 4, 2);
199         tesseral.initializeShortPeriodTerms(auxiliaryElements, PropagationType.MEAN, tesseral.getParameters());
200 
201         // Tesseral force model with default constructor
202         final DSSTForceModel tesseralDefault = new DSSTTesseral(earthFrame,
203                                                                 Constants.WGS84_EARTH_ANGULAR_VELOCITY,
204                                                                 provider);
205         tesseralDefault.initializeShortPeriodTerms(auxiliaryElements, PropagationType.MEAN, tesseralDefault.getParameters());
206 
207         // Compute mean element rate for the tesseral force model
208         final double[] elements = tesseral.getMeanElementRate(state, auxiliaryElements, tesseral.getParameters());
209 
210         // Compute mean element rate for the "default" tesseral force model
211         final double[] elementsDefault = tesseralDefault.getMeanElementRate(state, auxiliaryElements, tesseralDefault.getParameters());
212 
213         // Verify
214         for (int i = 0; i < 6; i++) {
215             Assert.assertEquals(elements[i], elementsDefault[i], Double.MIN_VALUE);
216         }
217 
218     }
219 
220     @Test
221     public void testIssue736() {
222 
223         // Central Body geopotential 4x4
224         final UnnormalizedSphericalHarmonicsProvider provider =
225                 GravityFieldFactory.getUnnormalizedProvider(4, 4);
226 
227         // Frames and epoch
228         final Frame frame = FramesFactory.getEME2000();
229         final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
230         final AbsoluteDate initDate = new AbsoluteDate(2007, 4, 16, 0, 46, 42.400, TimeScalesFactory.getUTC());
231 
232         // Initial orbit
233         final Orbit orbit = new EquinoctialOrbit(2.655989E7, 2.719455286199036E-4, 0.0041543085910249414,
234                                                  -0.3412974060023717, 0.3960084733107685,
235                                                  8.566537840341699, PositionAngle.TRUE,
236                                                  frame, initDate, provider.getMu());
237 
238         // Force model
239         final DSSTForceModel tesseral = new DSSTTesseral(earthFrame, Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider);
240         final double[] parameters = tesseral.getParameters();
241 
242         // Initialize force model
243         tesseral.initializeShortPeriodTerms(new AuxiliaryElements(orbit, 1), PropagationType.MEAN, parameters);
244 
245         // Eccentricity shift
246         final Orbit shfitedOrbit = new EquinoctialOrbit(2.655989E7, 0.02, 0.0041543085910249414,
247                                                  -0.3412974060023717, 0.3960084733107685,
248                                                  8.566537840341699, PositionAngle.TRUE,
249                                                  frame, initDate, provider.getMu());
250 
251         final double[] elements = tesseral.getMeanElementRate(new SpacecraftState(shfitedOrbit), new AuxiliaryElements(shfitedOrbit, 1), parameters);
252 
253         // The purpose of this test is not to verify a specific value.
254         // Its purpose is to verify that a NullPointerException does not
255         // occur when calculating initial values of Hansen Coefficients
256         for (int i = 0; i < elements.length; i++) {
257             Assert.assertTrue(elements[i] != 0);
258         }
259 
260     }
261 
262     @Test
263     public void testOutOfRangeException() {
264         // Central Body geopotential 1x0
265         final UnnormalizedSphericalHarmonicsProvider provider =
266                 GravityFieldFactory.getUnnormalizedProvider(1, 0);
267         // Earth
268         final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
269                                                             Constants.WGS84_EARTH_FLATTENING,
270                                                             FramesFactory.getGTOD(false));
271         try {
272             @SuppressWarnings("unused")
273             final DSSTForceModel tesseral = new DSSTTesseral(earth.getBodyFrame(),
274                                                           Constants.WGS84_EARTH_ANGULAR_VELOCITY,
275                                                           provider);
276             Assert.fail("An exception should have been thrown");
277         } catch (OrekitException oe) {
278             Assert.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
279         }
280     }
281 
282     @Test
283     public void testGetMaxEccPow()
284         throws NoSuchMethodException, SecurityException, IllegalAccessException, IllegalArgumentException, InvocationTargetException {
285         final UnnormalizedSphericalHarmonicsProvider provider =
286                         GravityFieldFactory.getUnnormalizedProvider(4, 4);;
287         final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
288         final DSSTTesseral force = new DSSTTesseral(earthFrame, Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider);
289         Method getMaxEccPow = DSSTTesseral.class.getDeclaredMethod("getMaxEccPow", Double.TYPE);
290         getMaxEccPow.setAccessible(true);
291         Assert.assertEquals(3,  getMaxEccPow.invoke(force, 0.0));
292         Assert.assertEquals(4,  getMaxEccPow.invoke(force, 0.01));
293         Assert.assertEquals(7,  getMaxEccPow.invoke(force, 0.08));
294         Assert.assertEquals(10, getMaxEccPow.invoke(force, 0.15));
295         Assert.assertEquals(12, getMaxEccPow.invoke(force, 0.25));
296         Assert.assertEquals(15, getMaxEccPow.invoke(force, 0.35));
297         Assert.assertEquals(20, getMaxEccPow.invoke(force, 1.0));
298     }
299 
300     @Before
301     public void setUp() throws IOException, ParseException {
302         Utils.setDataRoot("regular-data:potential/shm-format");
303     }
304 }