1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.forces.drag.atmosphere;
18
19 import org.hipparchus.analysis.UnivariateFunction;
20 import org.hipparchus.analysis.differentiation.DSFactory;
21 import org.hipparchus.analysis.differentiation.DerivativeStructure;
22 import org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator;
23 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
24 import org.hipparchus.geometry.euclidean.threed.Vector3D;
25 import org.junit.After;
26 import org.junit.Assert;
27 import org.junit.Before;
28 import org.junit.Test;
29 import org.orekit.Utils;
30 import org.orekit.bodies.CelestialBodyFactory;
31 import org.orekit.bodies.GeodeticPoint;
32 import org.orekit.bodies.OneAxisEllipsoid;
33 import org.orekit.errors.OrekitException;
34 import org.orekit.forces.drag.atmosphere.HarrisPriester;
35 import org.orekit.frames.Frame;
36 import org.orekit.frames.FramesFactory;
37 import org.orekit.time.AbsoluteDate;
38 import org.orekit.time.DateComponents;
39 import org.orekit.time.FieldAbsoluteDate;
40 import org.orekit.time.TimeComponents;
41 import org.orekit.time.TimeScale;
42 import org.orekit.time.TimeScalesFactory;
43 import org.orekit.utils.PVCoordinatesProvider;
44
45 public class HarrisPriesterTest {
46
47
48 private PVCoordinatesProvider sun;
49
50
51 private OneAxisEllipsoid earth;
52
53
54 private Frame earthFrame;
55
56
57 private TimeScale utc;
58
59
60 private AbsoluteDate date;
61
62 @Test
63 public void testStandard() {
64
65 final HarrisPriester hp = new HarrisPriester(sun, earth);
66
67
68 final GeodeticPoint point = new GeodeticPoint(0, 0, 500000.);
69 final Vector3D pos = earth.transform(point);
70
71
72 final double rho = hp.getDensity(date, pos, earthFrame);
73
74 Assert.assertEquals(3.9237E-13, rho, 1.0e-17);
75
76 }
77
78 @Test
79 public void testParameterN() {
80
81 final HarrisPriester hp = new HarrisPriester(sun, earth);
82
83
84 final GeodeticPoint point = new GeodeticPoint(0, 0, 500000.);
85 final Vector3D pos = earth.transform(point);
86
87
88 final double rho4 = hp.getDensity(date, pos, earthFrame);
89
90
91 final HarrisPriester hp2 = new HarrisPriester(sun, earth, 2);
92
93
94 final double rho2 = hp2.getDensity(date, pos, earthFrame);
95
96 final HarrisPriester hp6 = new HarrisPriester(sun, earth, 6);
97
98
99 final double rho6 = hp6.getDensity(date, pos, earthFrame);
100
101 final double c2Psi2 = 0.02163787;
102
103 Assert.assertEquals(c2Psi2, (rho6-rho2)/(rho4-rho2) - 1., 1.e-8);
104
105 }
106
107 @Test
108 public void testMaxAlt() {
109
110 final HarrisPriester hp = new HarrisPriester(sun, earth);
111
112
113 final GeodeticPoint point = new GeodeticPoint(0, 0, 1500000.);
114 final Vector3D pos = earth.transform(point);
115
116
117 final double rho = hp.getDensity(date, pos, earthFrame);
118
119 Assert.assertEquals(0.0, rho, 0.0);
120 }
121
122 @Test
123 public void testUserTab() {
124
125 final double[][] userTab = {
126 {100000., 4.974e+02, 4.974e+02},
127 {110000., 7.800e+01, 7.800e+01},
128 {120000., 2.490e+01, 2.400e+01},
129 {130000., 8.377e+00, 8.710e+00},
130 {140000., 3.899e+00, 4.059e+00},
131 {150000., 2.122e+00, 2.215e+00},
132 {160000., 1.263e+00, 1.344e+00},
133 {170000., 8.008e-01, 8.758e-01},
134 {180000., 5.283e-01, 6.010e-01},
135 {190000., 3.618e-01, 4.297e-01},
136 {200000., 2.557e-01, 3.162e-01},
137 {210000., 1.839e-01, 2.396e-01},
138 {220000., 1.341e-01, 1.853e-01},
139 {230000., 9.949e-02, 1.455e-01},
140 {240000., 7.488e-02, 1.157e-01},
141 {250000., 5.709e-02, 9.308e-02},
142 {260000., 4.403e-02, 7.555e-02},
143 {270000., 3.430e-02, 6.182e-02},
144 {280000., 2.697e-02, 5.095e-02},
145 {290000., 2.139e-02, 4.226e-02},
146 {300000., 1.708e-02, 3.526e-02},
147 {320000., 1.099e-02, 2.511e-02},
148 {340000., 7.214e-03, 1.819e-02},
149 {360000., 4.824e-03, 1.337e-02},
150 {380000., 3.274e-03, 9.955e-03},
151 {400000., 2.249e-03, 7.492e-03},
152 {420000., 1.558e-03, 5.684e-03},
153 {440000., 1.091e-03, 4.355e-03},
154 {460000., 7.701e-04, 3.362e-03},
155 {480000., 5.474e-04, 2.612e-03},
156 {500000., 3.916e-04, 2.042e-03},
157 {520000., 2.819e-04, 1.605e-03},
158 {540000., 2.042e-04, 1.267e-03},
159 {560000., 1.488e-04, 1.005e-03},
160 {580000., 1.092e-04, 7.997e-04},
161 {600000., 8.070e-05, 6.390e-04},
162 {620000., 6.012e-05, 5.123e-04},
163 {640000., 4.519e-05, 4.121e-04},
164 {660000., 3.430e-05, 3.325e-04},
165 {680000., 2.632e-05, 2.691e-04},
166 {700000., 2.043e-05, 2.185e-04},
167 {720000., 1.607e-05, 1.779e-04},
168 {740000., 1.281e-05, 1.452e-04},
169 {760000., 1.036e-05, 1.190e-04},
170 {780000., 8.496e-06, 9.776e-05},
171 {800000., 7.069e-06, 8.059e-05},
172 {850000., 4.800e-06, 5.500e-05},
173 {900000., 3.300e-06, 3.700e-05},
174 {950000., 2.450e-06, 2.400e-05},
175 {1000000., 1.900e-06, 1.700e-05},
176 {1100000., 1.180e-06, 8.700e-06},
177 {1200000., 7.500e-07, 4.800e-06},
178 {1300000., 5.300e-07, 3.200e-06},
179 {1400000., 4.100e-07, 2.000e-06},
180 {1500000., 2.900e-07, 1.350e-06},
181 {1600000., 2.000e-07, 9.500e-07},
182 {1700000., 1.600e-07, 7.700e-07},
183 {1800000., 1.200e-07, 6.300e-07},
184 {1900000., 9.600e-08, 5.200e-07},
185 {2000000., 7.300e-08, 4.400e-07}
186 };
187
188
189 final GeodeticPoint point = new GeodeticPoint(0, 0, 1500000.);
190 final Vector3D pos = earth.transform(point);
191
192 final HarrisPriester hp = new HarrisPriester(sun, earth, userTab);
193
194
195 final double rho = hp.getDensity(date, pos, earthFrame);
196
197 Assert.assertEquals(2.9049E-7, rho, 1.0e-11);
198
199 final HarrisPriester hp6 = new HarrisPriester(sun, earth, userTab, 6);
200 final double rho6 = hp6.getDensity(date, pos, earthFrame);
201
202 final HarrisPriester hp2 = new HarrisPriester(sun, earth, userTab, 2);
203 final double rho2 = hp2.getDensity(date, pos, earthFrame);
204
205 final double c2Psi2 = 0.02163787;
206
207 Assert.assertEquals(c2Psi2, (rho6-rho2)/(rho-rho2) - 1., 1.0e-8);
208
209 }
210
211 @Test(expected=OrekitException.class)
212 public void testOutOfRange() {
213
214 final HarrisPriester hp = new HarrisPriester(sun, earth);
215
216
217 final GeodeticPoint point = new GeodeticPoint(0, 0, 50000.);
218 final Vector3D pos = earth.transform(point);
219
220
221 hp.getDensity(date, pos, earthFrame);
222 }
223
224 @Test
225 public void testVelocityDerivative() {
226 final Frame eme2000 = FramesFactory.getEME2000();
227 final HarrisPriester hp = new HarrisPriester(sun, earth);
228 final Vector3D pos = earth.getBodyFrame().getTransformTo(eme2000, date).
229 transformPosition(earth.transform(new GeodeticPoint(-1.7, 4.2, 987654.321)));
230 double dP = 100.0;
231 double dVxdX = gradientComponent(hp, pos, Vector3D.PLUS_I, dP, eme2000, v -> v.getX());
232 double dVxdY = gradientComponent(hp, pos, Vector3D.PLUS_J, dP, eme2000, v -> v.getX());
233 double dVxdZ = gradientComponent(hp, pos, Vector3D.PLUS_K, dP, eme2000, v -> v.getX());
234 double dVydX = gradientComponent(hp, pos, Vector3D.PLUS_I, dP, eme2000, v -> v.getY());
235 double dVydY = gradientComponent(hp, pos, Vector3D.PLUS_J, dP, eme2000, v -> v.getY());
236 double dVydZ = gradientComponent(hp, pos, Vector3D.PLUS_K, dP, eme2000, v -> v.getY());
237 double dVzdX = gradientComponent(hp, pos, Vector3D.PLUS_I, dP, eme2000, v -> v.getZ());
238 double dVzdY = gradientComponent(hp, pos, Vector3D.PLUS_J, dP, eme2000, v -> v.getZ());
239 double dVzdZ = gradientComponent(hp, pos, Vector3D.PLUS_K, dP, eme2000, v -> v.getZ());
240
241 DSFactory factory = new DSFactory(3, 1);
242 FieldVector3D<DerivativeStructure> dsPos = new FieldVector3D<>(factory.variable(0, pos.getX()),
243 factory.variable(1, pos.getY()),
244 factory.variable(2, pos.getZ()));
245 FieldVector3D<DerivativeStructure> dsVel = hp.getVelocity(new FieldAbsoluteDate<>(factory.getDerivativeField(),
246 date),
247 dsPos, eme2000);
248 Assert.assertEquals(dVxdX, dsVel.getX().getPartialDerivative(1, 0, 0), 1.0e-16);
249 Assert.assertEquals(dVxdY, dsVel.getX().getPartialDerivative(0, 1, 0), 1.0e-16);
250 Assert.assertEquals(dVxdZ, dsVel.getX().getPartialDerivative(0, 0, 1), 1.0e-16);
251 Assert.assertEquals(dVydX, dsVel.getY().getPartialDerivative(1, 0, 0), 1.0e-16);
252 Assert.assertEquals(dVydY, dsVel.getY().getPartialDerivative(0, 1, 0), 1.0e-16);
253 Assert.assertEquals(dVydZ, dsVel.getY().getPartialDerivative(0, 0, 1), 1.0e-16);
254 Assert.assertEquals(dVzdX, dsVel.getZ().getPartialDerivative(1, 0, 0), 1.0e-16);
255 Assert.assertEquals(dVzdY, dsVel.getZ().getPartialDerivative(0, 1, 0), 1.0e-16);
256 Assert.assertEquals(dVzdZ, dsVel.getZ().getPartialDerivative(0, 0, 1), 1.0e-16);
257
258 }
259
260 private double gradientComponent(final Atmosphere atm, final Vector3D position, final Vector3D direction,
261 final double dP, final Frame frame, final ComponentGetter getter) {
262 FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(5, dP);
263 UnivariateFunction f = delta -> {
264 try {
265 return getter.get(atm.getVelocity(date, new Vector3D(1, position, delta, direction), frame));
266 } catch (OrekitException oe) {
267 return Double.NaN;
268 }
269 };
270 return differentiator.differentiate(f).value(new DSFactory(1, 1).variable(0, 0.0)).getPartialDerivative(1);
271 }
272
273 private interface ComponentGetter {
274 double get(Vector3D v);
275 }
276
277 @Before
278 public void setUp() {
279 Utils.setDataRoot("regular-data");
280 sun = CelestialBodyFactory.getSun();
281
282 earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
283 earth = new OneAxisEllipsoid(6378136.460, 1.0 / 298.257222101, earthFrame);
284
285
286 utc = TimeScalesFactory.getUTC();
287 date = new AbsoluteDate(new DateComponents(2003, 03, 21), new TimeComponents(1, 0, 0.), utc);
288 }
289
290 @After
291 public void tearDown() {
292 utc = null;
293 }
294
295 }