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.text.ParseException;
21  import java.util.ArrayList;
22  import java.util.Arrays;
23  import java.util.List;
24  
25  import org.hipparchus.exception.LocalizedCoreFormats;
26  import org.hipparchus.util.FastMath;
27  import org.junit.Assert;
28  import org.junit.Before;
29  import org.junit.Test;
30  import org.orekit.Utils;
31  import org.orekit.errors.OrekitException;
32  import org.orekit.forces.gravity.potential.GRGSFormatReader;
33  import org.orekit.forces.gravity.potential.GravityFieldFactory;
34  import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
35  import org.orekit.frames.Frame;
36  import org.orekit.frames.FramesFactory;
37  import org.orekit.orbits.EquinoctialOrbit;
38  import org.orekit.orbits.Orbit;
39  import org.orekit.orbits.PositionAngle;
40  import org.orekit.propagation.PropagationType;
41  import org.orekit.propagation.SpacecraftState;
42  import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
43  import org.orekit.propagation.semianalytical.dsst.forces.DSSTZonal;
44  import org.orekit.propagation.semianalytical.dsst.forces.ShortPeriodTerms;
45  import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
46  import org.orekit.time.AbsoluteDate;
47  import org.orekit.time.DateComponents;
48  import org.orekit.time.TimeComponents;
49  import org.orekit.time.TimeScalesFactory;
50  
51  public class DSSTZonalTest {
52  
53      @Test
54      public void testGetMeanElementRate() {
55          
56          // Central Body geopotential 4x4
57          final UnnormalizedSphericalHarmonicsProvider provider =
58                  GravityFieldFactory.getUnnormalizedProvider(4, 4);
59          
60          final Frame earthFrame = FramesFactory.getEME2000();
61          final AbsoluteDate initDate = new AbsoluteDate(2007, 04, 16, 0, 46, 42.400, TimeScalesFactory.getUTC());
62          
63          // a  = 26559890 m
64          // ex = 2.719455286199036E-4
65          // ey = 0.0041543085910249414
66          // hx = -0.3412974060023717
67          // hy = 0.3960084733107685
68          // lM = 8.566537840341699 rad
69          final Orbit orbit = new EquinoctialOrbit(2.655989E7,
70                                                   2.719455286199036E-4,
71                                                   0.0041543085910249414,
72                                                   -0.3412974060023717,
73                                                   0.3960084733107685,
74                                                   8.566537840341699,
75                                                   PositionAngle.TRUE,
76                                                   earthFrame,
77                                                   initDate,
78                                                   3.986004415E14);
79          
80          final SpacecraftState state = new SpacecraftState(orbit, 1000.0);
81          
82          final DSSTForceModel zonal = new DSSTZonal(provider, 4, 3, 9);
83  
84          // Force model parameters
85          final double[] parameters = zonal.getParameters();
86  
87          final AuxiliaryElements auxiliaryElements = new AuxiliaryElements(state.getOrbit(), 1);
88          
89          // Initialize force model
90          zonal.initializeShortPeriodTerms(auxiliaryElements, PropagationType.MEAN, parameters);
91          
92          final double[] elements = new double[7];
93          Arrays.fill(elements, 0.0);
94          
95          final double[] daidt = zonal.getMeanElementRate(state, auxiliaryElements, parameters);
96          for (int i = 0; i < daidt.length; i++) {
97              elements[i] = daidt[i];
98          }
99          
100         Assert.assertEquals(0.0,                     elements[0], 1.e-25);
101         Assert.assertEquals(1.3909396722346468E-11,  elements[1], 3.e-26);
102         Assert.assertEquals(-2.0275977261372793E-13, elements[2], 3.e-27);
103         Assert.assertEquals(3.087141512018238E-9,    elements[3], 1.e-24);
104         Assert.assertEquals(2.6606317310148797E-9,   elements[4], 4.e-24);
105         Assert.assertEquals(-3.659904725206694E-9,   elements[5], 1.e-24);
106         
107     }
108     
109     @Test
110     public void testShortPeriodTerms() {
111         final SpacecraftState meanState = getGEOState();
112         
113         final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(2, 0);
114         final DSSTForceModel zonal    = new DSSTZonal(provider, 2, 1, 5);
115         
116         //Create the auxiliary object
117         final AuxiliaryElements aux = new AuxiliaryElements(meanState.getOrbit(), 1);
118 
119         // Set the force models
120         final List<ShortPeriodTerms> shortPeriodTerms = new ArrayList<ShortPeriodTerms>();
121 
122         zonal.registerAttitudeProvider(null);
123         shortPeriodTerms.addAll(zonal.initializeShortPeriodTerms(aux, PropagationType.OSCULATING, zonal.getParameters()));
124         zonal.updateShortPeriodTerms(zonal.getParameters(), meanState);
125 
126         double[] y = new double[6];
127         for (final ShortPeriodTerms spt : shortPeriodTerms) {
128             final double[] shortPeriodic = spt.value(meanState.getOrbit());
129             for (int i = 0; i < shortPeriodic.length; i++) {
130                 y[i] += shortPeriodic[i];
131             }
132         }
133         
134         Assert.assertEquals(35.005618980090276,     y[0], 1.e-15);
135         Assert.assertEquals(3.75891551882889E-5,    y[1], 1.e-20);
136         Assert.assertEquals(3.929119925563796E-6,   y[2], 1.e-21);
137         Assert.assertEquals(-1.1781951949124315E-8, y[3], 1.e-24);
138         Assert.assertEquals(-3.2134924513679615E-8, y[4], 1.e-24);
139         Assert.assertEquals(-1.1607392915997098E-6, y[5], 1.e-21);
140     }
141 
142     @Test
143     public void testIssue625() {
144 
145         Utils.setDataRoot("regular-data:potential/grgs-format");
146         GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
147 
148         final Frame earthFrame = FramesFactory.getEME2000();
149         final AbsoluteDate initDate = new AbsoluteDate(2007, 04, 16, 0, 46, 42.400, TimeScalesFactory.getUTC());
150         
151         // a  = 2.655989E6 m
152         // ex = 2.719455286199036E-4
153         // ey = 0.0041543085910249414
154         // hx = -0.3412974060023717
155         // hy = 0.3960084733107685
156         // lM = 8.566537840341699 rad
157         final Orbit orbit = new EquinoctialOrbit(2.655989E6,
158                                                  2.719455286199036E-4,
159                                                  0.0041543085910249414,
160                                                  -0.3412974060023717,
161                                                  0.3960084733107685,
162                                                  8.566537840341699,
163                                                  PositionAngle.TRUE,
164                                                  earthFrame,
165                                                  initDate,
166                                                  3.986004415E14);
167         
168         final SpacecraftState state = new SpacecraftState(orbit, 1000.0);
169 
170         final AuxiliaryElements auxiliaryElements = new AuxiliaryElements(state.getOrbit(), 1);
171 
172         // Central Body geopotential 32x32
173         final UnnormalizedSphericalHarmonicsProvider provider =
174                 GravityFieldFactory.getUnnormalizedProvider(32, 32);
175 
176         // Zonal force model
177         final DSSTZonal zonal = new DSSTZonal(provider, 32, 4, 65);
178         zonal.initializeShortPeriodTerms(auxiliaryElements, PropagationType.MEAN, zonal.getParameters());
179 
180         // Zonal force model with default constructor
181         final DSSTZonal zonalDefault = new DSSTZonal(provider);
182         zonalDefault.initializeShortPeriodTerms(auxiliaryElements, PropagationType.MEAN, zonalDefault.getParameters());
183 
184         // Compute mean element rate for the zonal force model
185         final double[] elements = zonal.getMeanElementRate(state, auxiliaryElements, zonal.getParameters());
186 
187         // Compute mean element rate for the "default" zonal force model
188         final double[] elementsDefault = zonalDefault.getMeanElementRate(state, auxiliaryElements, zonalDefault.getParameters());
189 
190         // Verify
191         for (int i = 0; i < 6; i++) {
192             Assert.assertEquals(elements[i], elementsDefault[i], Double.MIN_VALUE);
193         }
194 
195     }
196 
197     @Test
198     public void testOutOfRangeException() {
199         try {
200             @SuppressWarnings("unused")
201             final DSSTZonal zonal = new DSSTZonal(GravityFieldFactory.getUnnormalizedProvider(1, 0));
202             Assert.fail("An exception should have been thrown");
203         } catch (OrekitException oe) {
204             Assert.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
205         }
206     }
207 
208     private SpacecraftState getGEOState() {
209         // No shadow at this date
210         final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2003, 05, 21), new TimeComponents(1, 0, 0.),
211                                                        TimeScalesFactory.getUTC());
212         final Orbit orbit = new EquinoctialOrbit(42164000,
213                                                  10e-3,
214                                                  10e-3,
215                                                  FastMath.tan(0.001745329) * FastMath.cos(2 * FastMath.PI / 3),
216                                                  FastMath.tan(0.001745329) * FastMath.sin(2 * FastMath.PI / 3), 0.1,
217                                                  PositionAngle.TRUE,
218                                                  FramesFactory.getEME2000(),
219                                                  initDate,
220                                                  3.986004415E14);
221         return new SpacecraftState(orbit);
222     }
223     
224     @Before
225     public void setUp() throws IOException, ParseException {
226         Utils.setDataRoot("regular-data:potential/shm-format");
227     }
228 
229 }