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