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