1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.propagation.events;
18
19 import java.util.ArrayList;
20 import java.util.Comparator;
21
22 import org.hipparchus.geometry.euclidean.twod.PolygonsSet;
23 import org.hipparchus.geometry.euclidean.twod.Vector2D;
24 import org.hipparchus.geometry.partitioning.Region.Location;
25 import org.hipparchus.ode.events.Action;
26 import org.hipparchus.util.FastMath;
27 import org.junit.jupiter.api.AfterEach;
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.bodies.GeodeticPoint;
33 import org.orekit.bodies.OneAxisEllipsoid;
34 import org.orekit.data.DataContext;
35 import org.orekit.frames.Frame;
36 import org.orekit.frames.FramesFactory;
37 import org.orekit.models.earth.GeoMagneticField;
38 import org.orekit.models.earth.GeoMagneticFieldFactory;
39 import org.orekit.models.earth.GeoMagneticFieldFactory.FieldModel;
40 import org.orekit.orbits.CircularOrbit;
41 import org.orekit.orbits.Orbit;
42 import org.orekit.orbits.PositionAngleType;
43 import org.orekit.propagation.Propagator;
44 import org.orekit.propagation.SpacecraftState;
45 import org.orekit.propagation.analytical.KeplerianPropagator;
46 import org.orekit.propagation.events.handlers.EventHandler;
47 import org.orekit.time.AbsoluteDate;
48 import org.orekit.time.TimeScale;
49 import org.orekit.time.TimeScalesFactory;
50 import org.orekit.utils.Constants;
51 import org.orekit.utils.IERSConventions;
52 import org.orekit.utils.units.UnitsConverter;
53
54
55
56
57
58 public class MagneticFieldDetectorTest {
59
60 Propagator propagator;
61 Frame eme2000;
62 Frame itrf;
63 TimeScale utc;
64 AbsoluteDate initialDate;
65 OneAxisEllipsoid earth;
66 GeoMagneticField wmm;
67 GeoMagneticField igrf;
68
69 double saaValidationThreshold = UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(500);
70 double saaValidationWidth = UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(200);
71
72 @BeforeEach
73 public void setup() {
74
75 Utils.setDataRoot("regular-data:earth");
76
77
78 eme2000 = FramesFactory.getEME2000();
79 itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, false);
80 utc = TimeScalesFactory.getUTC();
81 initialDate = new AbsoluteDate("2019-01-01T12:00:00.000", utc);
82 earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, itrf);
83
84
85 wmm = GeoMagneticFieldFactory.getField(FieldModel.WMM, 2019);
86 igrf = GeoMagneticFieldFactory.getField(FieldModel.IGRF, 2019);
87
88 }
89
90 @AfterEach
91 public void tearDown() {
92
93 DataContext.getDefault().getDataProvidersManager().clearProviders();
94 DataContext.getDefault().getDataProvidersManager().clearLoadedDataNames();
95 }
96
97
98
99
100 private void initializePropagator() {
101 double a = Constants.WGS84_EARTH_EQUATORIAL_RADIUS + 600000;
102 double i = FastMath.toRadians(80);
103 Orbit initialOrbit = new CircularOrbit(a, 0, 0, i, 0, 0, PositionAngleType.TRUE, eme2000, initialDate, Constants.WGS84_EARTH_MU);
104 propagator = new KeplerianPropagator(initialOrbit);
105 }
106
107
108
109
110 @Test
111 public void magneticFieldDetectorWMMSeaLevelTest() {
112 initializePropagator();
113 double threshold = UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(45000);
114
115 CustomEventHandler handler = new CustomEventHandler();
116 MagneticFieldDetector detector = new MagneticFieldDetector(threshold, FieldModel.WMM, earth, true).withHandler(handler);
117 propagator.addEventDetector(detector);
118 propagator.propagate(initialDate, initialDate.shiftedBy(864000));
119
120 checkEvents(handler.getEvents(), threshold, wmm, true);
121 }
122
123
124
125
126 @Test
127 public void magneticFieldDetectorIGRFSeaLevelTest() {
128 initializePropagator();
129 double threshold = UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(45000);
130
131 CustomEventHandler handler = new CustomEventHandler();
132 MagneticFieldDetector detector = new MagneticFieldDetector(threshold, FieldModel.IGRF, earth, true).withHandler(handler);
133 propagator.addEventDetector(detector);
134 propagator.propagate(initialDate, initialDate.shiftedBy(864000));
135
136 checkEvents(handler.getEvents(), threshold, igrf, true);
137 }
138
139
140
141
142 @Test
143 public void magneticFieldDetectorWMMTrueAltitudeTest() {
144 initializePropagator();
145 double threshold = UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(45000);
146
147 CustomEventHandler handler = new CustomEventHandler();
148 MagneticFieldDetector detector = new MagneticFieldDetector(threshold, FieldModel.WMM, earth).withHandler(handler);
149 propagator.addEventDetector(detector);
150 propagator.propagate(initialDate, initialDate.shiftedBy(864000));
151
152 checkEvents(handler.getEvents(), threshold, wmm, false);
153 }
154
155
156
157
158 @Test
159 public void magneticFieldDetectorIGRFTrueAltitudeTest() {
160 initializePropagator();
161 double threshold = UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(45000);
162
163 CustomEventHandler handler = new CustomEventHandler();
164 MagneticFieldDetector detector = new MagneticFieldDetector(threshold, FieldModel.IGRF, earth).withHandler(handler);
165 propagator.addEventDetector(detector);
166 propagator.propagate(initialDate, initialDate.shiftedBy(864000));
167
168 checkEvents(handler.getEvents(), threshold, igrf, false);
169 }
170
171
172
173
174 @Test
175 public void saaDetectorWMMSeaLevelTest() {
176 initializePropagator();
177
178 double altitude = 0;
179 double threshold = UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(32000);
180
181 PolygonsSet saaIn = generateGeomagneticMap(wmm, altitude, threshold - saaValidationThreshold, saaValidationWidth);
182 PolygonsSet saaOut = generateGeomagneticMap(wmm, altitude, threshold + saaValidationThreshold, saaValidationWidth);
183
184 CustomEventHandler handler = new CustomEventHandler();
185 MagneticFieldDetector detector = new MagneticFieldDetector(threshold, FieldModel.WMM, earth, true).withHandler(handler);
186 propagator.addEventDetector(detector);
187 propagator.propagate(initialDate, initialDate.shiftedBy(864000));
188
189 checkEvents(handler.getEvents(), threshold, wmm, true);
190 checkSaa(handler.getEvents(), saaIn, saaOut);
191 }
192
193
194
195
196 @Test
197 public void saaDetectorIGRFSeaLevelTest() {
198 initializePropagator();
199
200 double altitude = 0;
201 double threshold = UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(32000);
202
203 PolygonsSet saaIn = generateGeomagneticMap(igrf, altitude, threshold - saaValidationThreshold, saaValidationWidth);
204 PolygonsSet saaOut = generateGeomagneticMap(igrf, altitude, threshold + saaValidationThreshold, saaValidationWidth);
205
206 CustomEventHandler handler = new CustomEventHandler();
207 MagneticFieldDetector detector = new MagneticFieldDetector(threshold, FieldModel.IGRF, earth, true).withHandler(handler);
208 propagator.addEventDetector(detector);
209 propagator.propagate(initialDate, initialDate.shiftedBy(864000));
210
211 checkEvents(handler.getEvents(), threshold, igrf, true);
212 checkSaa(handler.getEvents(), saaIn, saaOut);
213 }
214
215
216
217
218 @Test
219 public void saaDetectorWMMTrueAltitudeTest() {
220 initializePropagator();
221
222 double altitude = 600000;
223 double threshold = UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(23000);
224
225 PolygonsSet saaIn = generateGeomagneticMap(wmm, altitude, threshold - saaValidationThreshold, saaValidationWidth);
226 PolygonsSet saaOut = generateGeomagneticMap(wmm, altitude, threshold + saaValidationThreshold, saaValidationWidth);
227
228 CustomEventHandler handler = new CustomEventHandler();
229 MagneticFieldDetector detector = new MagneticFieldDetector(threshold, FieldModel.WMM, earth).withHandler(handler);
230 propagator.addEventDetector(detector);
231 propagator.propagate(initialDate, initialDate.shiftedBy(864000));
232
233 checkEvents(handler.getEvents(), threshold, wmm, false);
234 checkSaa(handler.getEvents(), saaIn, saaOut);
235 }
236
237
238
239
240 @Test
241 public void saaDetectorIGRFTrueAltitudeTest() {
242 initializePropagator();
243
244 double altitude = 600000;
245 double threshold = UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(23000);
246
247 PolygonsSet saaIn = generateGeomagneticMap(igrf, altitude, threshold - saaValidationThreshold, saaValidationWidth);
248 PolygonsSet saaOut = generateGeomagneticMap(igrf, altitude, threshold + saaValidationThreshold, saaValidationWidth);
249
250 CustomEventHandler handler = new CustomEventHandler();
251 MagneticFieldDetector detector = new MagneticFieldDetector(threshold, FieldModel.IGRF, earth).withHandler(handler);
252 propagator.addEventDetector(detector);
253 propagator.propagate(initialDate, initialDate.shiftedBy(864000));
254
255 checkEvents(handler.getEvents(), threshold, igrf, false);
256 checkSaa(handler.getEvents(), saaIn, saaOut);
257 }
258
259
260
261
262
263
264
265
266
267
268 private PolygonsSet generateGeomagneticMap(GeoMagneticField field, double altitude, double threshold, double width) {
269
270
271 ArrayList<Double[]> points = new ArrayList<Double[]>();
272
273 for(int latitude = -89; latitude < 90; latitude++) {
274 for(int longitude = -179; longitude < 180; longitude++) {
275
276
277 double value = field.calculateField(FastMath.toRadians(latitude), FastMath.toRadians(longitude), altitude).getTotalIntensity();
278
279
280 if (value - threshold > -0.5*width && value - threshold < 0.5*width) {
281 Double[] point = {(double)latitude, (double)longitude};
282 points.add(point);
283 }
284 }
285 }
286
287
288 double[][] arrayPoints = new double[points.size()][2];
289
290 for(int i = 0; i < points.size(); i++) {
291 arrayPoints[i][0] = points.get(i)[0];
292 arrayPoints[i][1] = points.get(i)[1];
293 }
294
295
296 return buildZone(arrayPoints);
297 }
298
299
300
301
302
303
304 private PolygonsSet buildZone(double[][] points) {
305
306 final Vector2D[] vertices = new Vector2D[points.length];
307 for (int i = 0; i < points.length; ++i) {
308 vertices[i] = new Vector2D(FastMath.toRadians(points[i][1]),
309 FastMath.toRadians(points[i][0]));
310 }
311
312
313 final Vector2D[] sortedVertices = sortVertices(vertices);
314
315
316 return new PolygonsSet(1.0e-10, sortedVertices);
317 }
318
319
320
321
322
323
324
325
326
327 private Vector2D[] sortVertices(Vector2D[] vertices) {
328
329
330
331 double midLat = 0;
332 double midLon = 0;
333
334 double minLon = vertices[0].getX();
335 double maxLon = vertices[0].getX();
336 double minLat = vertices[0].getY();
337 double maxLat = vertices[0].getY();
338 for(int i = 0; i < vertices.length; i++) {
339 if(vertices[i].getX() < minLon) minLon = vertices[i].getX();
340 if(vertices[i].getX() > maxLon) maxLon = vertices[i].getX();
341 if(vertices[i].getY() < minLat) minLat = vertices[i].getY();
342 if(vertices[i].getY() > maxLat) maxLat = vertices[i].getY();
343 }
344 midLat = (minLat + maxLat)/2;
345 midLon = (minLon + maxLon)/2;
346
347 Vector2D mid = new Vector2D(midLon, midLat);
348
349
350 ArrayList<Vector2D> angularVerticesList = new ArrayList<Vector2D>();
351 Vector2D ref = new Vector2D(1, 0);
352 Vector2D angularVertice;
353 Vector2D centeredVertice;
354 for(int i = 0; i < vertices.length; i++) {
355 centeredVertice = vertices[i].subtract(mid);
356 double norm = centeredVertice.getNorm();
357 double angle = Vector2D.angle(ref, centeredVertice);
358 if(centeredVertice.getY() < 0) {
359 angle = -angle;
360 }
361 angularVertice = new Vector2D(norm, angle);
362 angularVerticesList.add(angularVertice);
363 }
364
365
366 Comparator<Vector2D> comparator = new Comparator<Vector2D>() {
367
368 @Override
369 public int compare(Vector2D o1, Vector2D o2) {
370 return Double.compare(o1.getY(),o2.getY());
371 }
372 };
373
374 angularVerticesList.sort(comparator);
375
376
377 Vector2D[] sortedVertices = new Vector2D[vertices.length];
378 Vector2D originalVertice;
379 for(int i = 0; i < vertices.length; i++) {
380 angularVertice = angularVerticesList.get(i);
381 double longitude = angularVertice.getX()*Math.cos(angularVertice.getY());
382 double latitude = angularVertice.getX()*Math.sin(angularVertice.getY());
383 centeredVertice = new Vector2D(longitude, latitude);
384 originalVertice = centeredVertice.add(mid);
385 sortedVertices[i] = originalVertice;
386 }
387
388 return sortedVertices;
389 }
390
391
392
393
394 private class CustomEventHandler implements EventHandler {
395
396 ArrayList<SpacecraftState> events = new ArrayList<SpacecraftState>();
397
398 @Override
399 public Action eventOccurred(SpacecraftState s, EventDetector detector, boolean increasing) {
400
401 events.add(s);
402
403 return Action.CONTINUE;
404 }
405
406 public ArrayList<SpacecraftState> getEvents(){
407 return events;
408 }
409 }
410
411
412
413
414
415
416
417
418 private void checkEvents(ArrayList<SpacecraftState> events, double threshold, GeoMagneticField field, boolean sea) {
419
420 for (SpacecraftState s : events) {
421
422 GeodeticPoint geo = earth.transform(s.getPosition(), s.getFrame(), s.getDate());
423 double altitude = geo.getAltitude();
424 if (sea) {
425 altitude = 0;
426 }
427 double meas = field.calculateField(geo.getLatitude(), geo.getLongitude(), altitude).getTotalIntensity();
428 Assertions.assertEquals(threshold, meas, 1e-12);
429 }
430 }
431
432
433
434
435
436
437
438 private void checkSaa(ArrayList<SpacecraftState> events, PolygonsSet saaIn, PolygonsSet saaOut) {
439
440 for (SpacecraftState s : events) {
441
442 GeodeticPoint geo = earth.transform(s.getPosition(itrf), itrf, s.getDate());
443 Vector2D point = new Vector2D(geo.getLongitude(), geo.getLatitude());
444
445
446 Assertions.assertTrue(saaIn.checkPoint(point).equals(Location.OUTSIDE));
447
448 Assertions.assertTrue(saaOut.checkPoint(point).equals(Location.INSIDE));
449
450 }
451 }
452 }