1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.propagation.semianalytical.dsst.forces;
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.hamcrest.MatcherAssert;
26 import org.hamcrest.Matchers;
27 import org.hipparchus.geometry.euclidean.threed.Rotation;
28 import org.hipparchus.geometry.euclidean.threed.RotationOrder;
29 import org.hipparchus.geometry.euclidean.threed.Vector3D;
30 import org.junit.jupiter.api.Assertions;
31 import org.junit.jupiter.api.BeforeEach;
32 import org.junit.jupiter.api.Test;
33 import org.orekit.Utils;
34 import org.orekit.attitudes.Attitude;
35 import org.orekit.attitudes.AttitudeProvider;
36 import org.orekit.attitudes.FrameAlignedProvider;
37 import org.orekit.attitudes.LofOffset;
38 import org.orekit.bodies.CelestialBody;
39 import org.orekit.bodies.CelestialBodyFactory;
40 import org.orekit.bodies.OneAxisEllipsoid;
41 import org.orekit.errors.OrekitException;
42 import org.orekit.forces.BoxAndSolarArraySpacecraft;
43 import org.orekit.forces.gravity.potential.GravityFieldFactory;
44 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
45 import org.orekit.frames.Frame;
46 import org.orekit.frames.FramesFactory;
47 import org.orekit.frames.LOFType;
48 import org.orekit.models.earth.atmosphere.Atmosphere;
49 import org.orekit.models.earth.atmosphere.HarrisPriester;
50 import org.orekit.models.earth.atmosphere.SimpleExponentialAtmosphere;
51 import org.orekit.orbits.EquinoctialOrbit;
52 import org.orekit.orbits.Orbit;
53 import org.orekit.orbits.PositionAngleType;
54 import org.orekit.propagation.PropagationType;
55 import org.orekit.propagation.SpacecraftState;
56 import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
57 import org.orekit.time.AbsoluteDate;
58 import org.orekit.time.DateComponents;
59 import org.orekit.time.TimeComponents;
60 import org.orekit.time.TimeScalesFactory;
61 import org.orekit.utils.Constants;
62 import org.orekit.utils.IERSConventions;
63 import org.orekit.utils.TimeStampedAngularCoordinates;
64
65 class DSSTAtmosphericDragTest {
66
67 @Test
68 void testGetMeanElementRate() throws IllegalArgumentException, OrekitException {
69
70
71 final UnnormalizedSphericalHarmonicsProvider provider =
72 GravityFieldFactory.getUnnormalizedProvider(2, 0);
73
74 final Frame earthFrame = FramesFactory.getEME2000();
75 final AbsoluteDate initDate = new AbsoluteDate(2003, 07, 01, 0, 0, 0, TimeScalesFactory.getUTC());
76 final double mu = 3.986004415E14;
77
78
79
80
81
82
83 final Orbit orbit = new EquinoctialOrbit(7204535.84810944,
84 -0.001119677138261611,
85 5.333650671984143E-4,
86 0.847841707880348,
87 0.7998014061193262,
88 3.897842092486239,
89 PositionAngleType.TRUE,
90 earthFrame,
91 initDate,
92 mu);
93
94
95 final OneAxisEllipsoid earth = new OneAxisEllipsoid(provider.getAe(),
96 Constants.WGS84_EARTH_FLATTENING,
97 CelestialBodyFactory.getEarth().getBodyOrientedFrame());
98 final Atmosphere atm = new HarrisPriester(CelestialBodyFactory.getSun(), earth, 6);
99 final double cd = 2.0;
100 final double area = 25.0;
101 DSSTForceModel drag = new DSSTAtmosphericDrag(atm, cd, area, mu);
102
103
104 Rotation rotation = new Rotation(1.0, 0.0, 0.0, 0.0, false);
105 Vector3D rotationRate = new Vector3D(0., 0., 0.);
106 Vector3D rotationAcceleration = new Vector3D(0., 0., 0.);
107 TimeStampedAngularCoordinates orientation = new TimeStampedAngularCoordinates(initDate, rotation, rotationRate, rotationAcceleration);
108 final Attitude att = new Attitude(earthFrame, orientation);
109
110 final SpacecraftState state = new SpacecraftState(orbit, att).withMass(1000.0);
111 final AuxiliaryElements auxiliaryElements = new AuxiliaryElements(state.getOrbit(), 1);
112
113
114 final double[] parameters = drag.getParameters(orbit.getDate());
115
116 drag.initializeShortPeriodTerms(auxiliaryElements, PropagationType.MEAN, parameters);
117
118
119 AttitudeProvider attitudeProvider = new FrameAlignedProvider(rotation);
120 drag.registerAttitudeProvider(attitudeProvider );
121
122
123 final double[] elements = new double[7];
124 Arrays.fill(elements, 0.0);
125 final double[] daidt = drag.getMeanElementRate(state, auxiliaryElements, parameters);
126 for (int i = 0; i < daidt.length; i++) {
127 elements[i] = daidt[i];
128 }
129
130 MatcherAssert.assertThat(elements[0], Matchers.closeTo(-3.415320567871035E-5, 4.e-20));
131 MatcherAssert.assertThat(elements[1], Matchers.closeTo(6.276312897745139E-13, 3e-26));
132 MatcherAssert.assertThat(elements[2], Matchers.closeTo(-9.303357008691404E-13, 3e-26));
133 MatcherAssert.assertThat(elements[3], Matchers.closeTo(-7.052316604063199E-14, 1.e-28));
134 MatcherAssert.assertThat(elements[4], Matchers.closeTo(-6.793277250493389E-14, 3.e-28));
135 MatcherAssert.assertThat(elements[5], Matchers.closeTo(-1.3565284454826392E-15, 1.e-27));
136
137 }
138
139 @Test
140 void testShortPeriodTerms() throws IllegalArgumentException, OrekitException {
141
142 final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2003, 03, 21), new TimeComponents(1, 0, 0.), TimeScalesFactory.getUTC());
143
144 final Orbit orbit = new EquinoctialOrbit(7069219.9806427825,
145 -4.5941811292223825E-4,
146 1.309932339472599E-4,
147 -1.002996107003202,
148 0.570979900577994,
149 2.62038786211518,
150 PositionAngleType.TRUE,
151 FramesFactory.getEME2000(),
152 initDate,
153 3.986004415E14);
154
155 final SpacecraftState meanState = new SpacecraftState(orbit);
156
157 final CelestialBody sun = CelestialBodyFactory.getSun();
158 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
159 Constants.WGS84_EARTH_FLATTENING,
160 FramesFactory.getITRF(IERSConventions.IERS_2010,
161 true));
162
163 final BoxAndSolarArraySpacecraft boxAndWing = new BoxAndSolarArraySpacecraft(5.0, 2.0, 2.0,
164 sun,
165 50.0, Vector3D.PLUS_J,
166 2.0, 0.1,
167 0.2, 0.6);
168
169 final Atmosphere atmosphere = new HarrisPriester(CelestialBodyFactory.getSun(), earth, 6);
170 final AttitudeProvider attitudeProvider = new LofOffset(meanState.getFrame(),
171 LOFType.LVLH_CCSDS, RotationOrder.XYZ,
172 0.0, 0.0, 0.0);
173
174 final DSSTForceModel drag = new DSSTAtmosphericDrag(atmosphere, boxAndWing, meanState.getOrbit().getMu());
175
176
177 final AuxiliaryElements aux = new AuxiliaryElements(meanState.getOrbit(), 1);
178
179
180 final List<ShortPeriodTerms> shortPeriodTerms = new ArrayList<ShortPeriodTerms>();
181
182 drag.registerAttitudeProvider(attitudeProvider);
183 shortPeriodTerms.addAll(drag.initializeShortPeriodTerms(aux, PropagationType.OSCULATING, drag.getParameters(meanState.getDate())));
184 drag.updateShortPeriodTerms(drag.getParametersAllValues(), meanState);
185
186 double[] y = new double[6];
187 for (final ShortPeriodTerms spt : shortPeriodTerms) {
188 final double[] shortPeriodic = spt.value(meanState.getOrbit());
189 for (int i = 0; i < shortPeriodic.length; i++) {
190 y[i] += shortPeriodic[i];
191 }
192 }
193
194 MatcherAssert.assertThat(y[0], Matchers.closeTo(0.03966657234224339, 1.e-15));
195 MatcherAssert.assertThat(y[1], Matchers.closeTo(-1.5294381443907498E-8, 1.e-23));
196 MatcherAssert.assertThat(y[2], Matchers.closeTo(-2.3614929828139117E-8, 1.e-23));
197 MatcherAssert.assertThat(y[3], Matchers.closeTo(-5.901580338876397E-11, 1.e-25));
198 MatcherAssert.assertThat(y[4], Matchers.closeTo(1.0287639737403387E-11, 1.e-24));
199 MatcherAssert.assertThat(y[5], Matchers.closeTo(2.5384275235864555E-8, 1.e-23));
200 }
201
202 @Test
203 public void testSetAtmosphereUpperLimit() {
204
205
206 final Frame eme2000Frame = FramesFactory.getEME2000();
207 final AbsoluteDate initDate = new AbsoluteDate(2003, 07, 01, 0, 0, 0, TimeScalesFactory.getUTC());
208 final double mu = 3.986004415E14;
209 final Orbit orbit = new EquinoctialOrbit(8204535.84810944,
210 -0.001119677138261611,
211 5.333650671984143E-4,
212 0.847841707880348,
213 0.7998014061193262,
214 3.897842092486239,
215 PositionAngleType.TRUE,
216 eme2000Frame,
217 initDate,
218 mu);
219
220
221 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
222 Constants.WGS84_EARTH_FLATTENING,
223 CelestialBodyFactory.getEarth().getBodyOrientedFrame());
224 final Atmosphere atm = new SimpleExponentialAtmosphere(earth, 2.6e-10, 200000, 26000);
225 final double cd = 2.0;
226 final double area = 25.0;
227 DSSTAtmosphericDrag drag = new DSSTAtmosphericDrag(atm, cd, area, mu);
228
229
230 Rotation rotation = new Rotation(1.0, 0.0, 0.0, 0.0, false);
231 Vector3D rotationRate = new Vector3D(0., 0., 0.);
232 Vector3D rotationAcceleration = new Vector3D(0., 0., 0.);
233 TimeStampedAngularCoordinates orientation = new TimeStampedAngularCoordinates(initDate, rotation, rotationRate, rotationAcceleration);
234 final Attitude att = new Attitude(eme2000Frame, orientation);
235
236 final SpacecraftState state = new SpacecraftState(orbit, att).withMass(1000.0);
237 final AuxiliaryElements auxiliaryElements = new AuxiliaryElements(state.getOrbit(), 1);
238
239
240 final double[] parameters = drag.getParameters(orbit.getDate());
241
242 drag.initializeShortPeriodTerms(auxiliaryElements, PropagationType.MEAN, parameters);
243
244
245 AttitudeProvider attitudeProvider = new FrameAlignedProvider(rotation);
246 drag.registerAttitudeProvider(attitudeProvider);
247
248
249 final double atmosphericMaxConstant = 1000000.0;
250 Assertions.assertEquals(atmosphericMaxConstant + Constants.WGS84_EARTH_EQUATORIAL_RADIUS, drag.getRbar(), 1e-9);
251
252
253 final double[] daidt = drag.getMeanElementRate(state, auxiliaryElements, parameters);
254 Assertions.assertEquals(0.0, daidt[0]);
255 Assertions.assertEquals(0.0, daidt[1]);
256 Assertions.assertEquals(0.0, daidt[2]);
257 Assertions.assertEquals(0.0, daidt[3]);
258 Assertions.assertEquals(0.0, daidt[4]);
259 Assertions.assertEquals(0.0, daidt[5]);
260
261
262 final double expectedAtmosphereLimit = 2000000.0 + Constants.WGS84_EARTH_EQUATORIAL_RADIUS;
263 drag.setRbar(expectedAtmosphereLimit);
264 Assertions.assertEquals(expectedAtmosphereLimit, drag.getRbar());
265
266
267 final double[] daidtNew = drag.getMeanElementRate(state, auxiliaryElements, parameters);
268 Assertions.assertEquals(-3.7465296003917817E-28, daidtNew[0], 1.0e-33);
269 Assertions.assertEquals(7.316039091705292E-36, daidtNew[1], 1.0e-41);
270 Assertions.assertEquals(-2.195983299144844E-36, daidtNew[2], 1.0e-41);
271 Assertions.assertEquals(-9.80724158695418E-37, daidtNew[3], 1.0e-42);
272 Assertions.assertEquals(-9.059767879911556E-37, daidtNew[4], 1.0e-42);
273 Assertions.assertEquals(1.4486591760431082E-38, daidtNew[5], 1.0e-43);
274 }
275
276 @BeforeEach
277 public void setUp() throws IOException, ParseException {
278 Utils.setDataRoot("regular-data:potential/shm-format");
279 }
280 }