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.Collection;
24  import java.util.List;
25  
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.bodies.CelestialBodyFactory;
32  import org.orekit.frames.Frame;
33  import org.orekit.frames.FramesFactory;
34  import org.orekit.orbits.EquinoctialOrbit;
35  import org.orekit.orbits.Orbit;
36  import org.orekit.orbits.PositionAngle;
37  import org.orekit.propagation.PropagationType;
38  import org.orekit.propagation.SpacecraftState;
39  import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
40  import org.orekit.propagation.semianalytical.dsst.forces.DSSTThirdBody;
41  import org.orekit.propagation.semianalytical.dsst.forces.ShortPeriodTerms;
42  import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
43  import org.orekit.time.AbsoluteDate;
44  import org.orekit.time.DateComponents;
45  import org.orekit.time.TimeComponents;
46  import org.orekit.time.TimeScalesFactory;
47  
48  public class DSSTThirdBodyTest {
49      
50      private static final double eps  = 3.5e-25;
51  
52      @Test
53      public void testGetMeanElementRate() throws IllegalArgumentException {
54          
55          final Frame earthFrame = FramesFactory.getEME2000();
56          final AbsoluteDate initDate = new AbsoluteDate(2003, 07, 01, 0, 0, 00.000, TimeScalesFactory.getUTC());
57          
58          final double mu = 3.986004415E14;
59          // a    = 42163393.0 m
60          // ex =  -0.25925449177598586
61          // ey =  -0.06946703170551687
62          // hx =   0.15995912655021305
63          // hy =  -0.5969755874197339
64          // lM   = 15.47576793123677 rad
65          final Orbit orbit = new EquinoctialOrbit(4.2163393E7,
66                                                   -0.25925449177598586,
67                                                   -0.06946703170551687,
68                                                   0.15995912655021305,
69                                                   -0.5969755874197339,
70                                                   15.47576793123677,
71                                                   PositionAngle.TRUE,
72                                                   earthFrame,
73                                                   initDate,
74                                                   mu);
75          
76          final SpacecraftState state = new SpacecraftState(orbit, 1000.0);
77          
78          final AuxiliaryElements auxiliaryElements = new AuxiliaryElements(state.getOrbit(), 1);
79          
80          final DSSTForceModel moon = new DSSTThirdBody(CelestialBodyFactory.getMoon(), mu);
81  
82          // Force model parameters
83          final double[] parameters = moon.getParameters();
84  
85          // Initialize force model
86          moon.initializeShortPeriodTerms(auxiliaryElements, PropagationType.MEAN, parameters);
87  
88          final double[] elements = new double[7];
89          Arrays.fill(elements, 0.0);
90          
91          final double[] daidt = moon.getMeanElementRate(state, auxiliaryElements, parameters);
92          for (int i = 0; i < daidt.length; i++) {
93              elements[i] = daidt[i];
94          }
95          
96          Assert.assertEquals(0.0,                    elements[0], eps);
97          Assert.assertEquals(4.346622384804537E-10,  elements[1], eps);
98          Assert.assertEquals(7.293879548440941E-10,  elements[2], eps);
99          Assert.assertEquals(7.465699631747887E-11,  elements[3], eps);
100         Assert.assertEquals(3.9170221137233836E-10, elements[4], eps);
101         Assert.assertEquals(-3.178319341840074E-10, elements[5], eps);
102 
103     }
104 
105     @Test
106     public void testShortPeriodTerms() throws IllegalArgumentException {
107         final SpacecraftState meanState = getGEOState();
108 
109         final DSSTForceModel moon = new DSSTThirdBody(CelestialBodyFactory.getMoon(), meanState.getMu());
110 
111         final Collection<DSSTForceModel> forces = new ArrayList<DSSTForceModel>();
112         forces.add(moon);
113 
114         //Create the auxiliary object
115         final AuxiliaryElements aux = new AuxiliaryElements(meanState.getOrbit(), 1);
116 
117         // Set the force models
118         final List<ShortPeriodTerms> shortPeriodTerms = new ArrayList<ShortPeriodTerms>();
119 
120         for (final DSSTForceModel force : forces) {
121             force.registerAttitudeProvider(null);
122             shortPeriodTerms.addAll(force.initializeShortPeriodTerms(aux, PropagationType.OSCULATING, force.getParameters()));
123             force.updateShortPeriodTerms(force.getParameters(), meanState);
124         }
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(-413.20633326933154,    y[0], 1.e-14);
135         Assert.assertEquals(-1.8060137920197483E-5, y[1], 1.e-20);
136         Assert.assertEquals(-2.8416367511811057E-5, y[2], 1.e-20);
137         Assert.assertEquals(-2.791424363476855E-6,  y[3], 1.e-21);
138         Assert.assertEquals(1.8817187527805853E-6,  y[4], 1.e-21);
139         Assert.assertEquals(-3.423664701811889E-5,  y[5], 1.e-20);
140     }
141 
142     private SpacecraftState getGEOState() throws IllegalArgumentException {
143         // No shadow at this date
144         final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2003, 05, 21), new TimeComponents(1, 0, 0.),
145                                                        TimeScalesFactory.getUTC());
146         final Orbit orbit = new EquinoctialOrbit(42164000,
147                                                  10e-3,
148                                                  10e-3,
149                                                  FastMath.tan(0.001745329) * FastMath.cos(2 * FastMath.PI / 3),
150                                                  FastMath.tan(0.001745329) * FastMath.sin(2 * FastMath.PI / 3), 0.1,
151                                                  PositionAngle.TRUE,
152                                                  FramesFactory.getEME2000(),
153                                                  initDate,
154                                                  3.986004415E14);
155         return new SpacecraftState(orbit);
156     }
157 
158     @Before
159     public void setUp() throws IOException, ParseException {
160         Utils.setDataRoot("regular-data");
161     }
162 
163 }