1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.ssa.collision.shorttermencounter.probability.twod;
18
19 import org.hipparchus.Field;
20 import org.hipparchus.analysis.differentiation.DSFactory;
21 import org.hipparchus.analysis.differentiation.DerivativeStructure;
22 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
23 import org.hipparchus.geometry.euclidean.threed.Vector3D;
24 import org.hipparchus.geometry.euclidean.twod.FieldVector2D;
25 import org.hipparchus.geometry.euclidean.twod.Vector2D;
26 import org.hipparchus.linear.BlockFieldMatrix;
27 import org.hipparchus.linear.BlockRealMatrix;
28 import org.hipparchus.linear.FieldMatrix;
29 import org.hipparchus.linear.RealMatrix;
30 import org.hipparchus.stat.descriptive.DescriptiveStatistics;
31 import org.hipparchus.util.Binary64;
32 import org.hipparchus.util.Binary64Field;
33 import org.junit.jupiter.api.Assertions;
34 import org.junit.jupiter.api.BeforeAll;
35 import org.junit.jupiter.api.Disabled;
36 import org.junit.jupiter.api.DisplayName;
37 import org.junit.jupiter.api.Test;
38 import org.mockito.Mockito;
39 import org.orekit.Utils;
40 import org.orekit.data.DataSource;
41 import org.orekit.files.ccsds.ndm.ParserBuilder;
42 import org.orekit.files.ccsds.ndm.cdm.Cdm;
43 import org.orekit.frames.Frame;
44 import org.orekit.frames.FramesFactory;
45 import org.orekit.frames.LOFType;
46 import org.orekit.orbits.CartesianOrbit;
47 import org.orekit.orbits.FieldCartesianOrbit;
48 import org.orekit.orbits.FieldOrbit;
49 import org.orekit.orbits.Orbit;
50 import org.orekit.propagation.FieldStateCovariance;
51 import org.orekit.propagation.StateCovariance;
52 import org.orekit.ssa.collision.shorttermencounter.probability.twod.armellinutils.ArmellinDataLoader;
53 import org.orekit.ssa.collision.shorttermencounter.probability.twod.armellinutils.ArmellinDataRow;
54 import org.orekit.ssa.collision.shorttermencounter.probability.twod.armellinutils.ArmellinStatistics;
55 import org.orekit.ssa.metrics.FieldProbabilityOfCollision;
56 import org.orekit.ssa.metrics.ProbabilityOfCollision;
57 import org.orekit.time.AbsoluteDate;
58 import org.orekit.time.FieldAbsoluteDate;
59 import org.orekit.utils.Constants;
60 import org.orekit.utils.FieldPVCoordinates;
61 import org.orekit.utils.PVCoordinates;
62
63 import java.io.IOException;
64 import java.util.List;
65
66 class Alfano2005Test {
67
68
69
70
71 private final ShortTermEncounter2DPOCMethod method = new Alfano2005();
72
73 @BeforeAll
74 static void initializeOrekitData() {
75 Utils.setDataRoot("regular-data");
76 }
77
78 @Test
79 @DisplayName("Chan test case 01")
80 void ChanTestCase01() {
81
82 final double xm = 0;
83 final double ym = 10;
84 final double sigmaX = 25;
85 final double sigmaY = 50;
86 final double radius = 5;
87
88
89 final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
90
91
92 Assertions.assertEquals(9.742e-3, result.getValue(), 1e-6);
93 }
94
95 @Test
96 @DisplayName("Chan test case 02")
97 void ChanTestCase02() {
98
99 final double xm = 10;
100 final double ym = 0;
101 final double sigmaX = 25;
102 final double sigmaY = 50;
103 final double radius = 5;
104
105
106 final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
107
108
109 Assertions.assertEquals(9.181e-3, result.getValue(), 1e-6);
110 }
111
112 @Test
113 @DisplayName("Chan test case 03")
114 void ChanTestCase03() {
115
116 final double xm = 0;
117 final double ym = 10;
118 final double sigmaX = 25;
119 final double sigmaY = 75;
120 final double radius = 5;
121
122
123 final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
124
125
126 Assertions.assertEquals(6.571e-3, result.getValue(), 1e-6);
127 }
128
129 @Test
130 @DisplayName("Chan test case 04")
131 void ChanTestCase04() {
132
133 final double xm = 10;
134 final double ym = 0;
135 final double sigmaX = 25;
136 final double sigmaY = 75;
137 final double radius = 5;
138
139
140 final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
141
142
143 Assertions.assertEquals(6.125e-3, result.getValue(), 1e-6);
144 }
145
146 @Test
147 @DisplayName("Chan test case 05")
148 void ChanTestCase05() {
149
150 final double xm = 0;
151 final double ym = 1000;
152 final double sigmaX = 1000;
153 final double sigmaY = 3000;
154 final double radius = 10;
155
156
157 final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
158
159
160 Assertions.assertEquals(1.577e-5, result.getValue(), 1e-8);
161 }
162
163 @Test
164 @DisplayName("Chan test case 06")
165 void ChanTestCase06() {
166
167 final double xm = 1000;
168 final double ym = 0;
169 final double sigmaX = 1000;
170 final double sigmaY = 3000;
171 final double radius = 10;
172
173
174 final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
175
176
177 Assertions.assertEquals(1.011e-5, result.getValue(), 1e-8);
178 }
179
180 @Test
181 @DisplayName("Chan test case 07")
182 void ChanTestCase07() {
183
184 final double xm = 0;
185 final double ym = 10000;
186 final double sigmaX = 1000;
187 final double sigmaY = 3000;
188 final double radius = 10;
189
190
191 final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
192
193
194 Assertions.assertEquals(6.443e-8, result.getValue(), 1e-11);
195 }
196
197 @Test
198 @DisplayName("Chan test case 08")
199 void ChanTestCase08() {
200
201 final double xm = 10000;
202 final double ym = 0;
203 final double sigmaX = 1000;
204 final double sigmaY = 3000;
205 final double radius = 10;
206
207
208 final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
209
210
211 Assertions.assertEquals(3.219e-27, result.getValue(), 1e-30);
212 }
213
214 @Test
215 @DisplayName("Chan test case 09")
216 void ChanTestCase09() {
217
218 final double xm = 0;
219 final double ym = 10000;
220 final double sigmaX = 1000;
221 final double sigmaY = 10000;
222 final double radius = 10;
223
224
225 final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
226
227
228 Assertions.assertEquals(3.033e-6, result.getValue(), 1e-9);
229 }
230
231 @Test
232 @DisplayName("Chan test case 10")
233 void ChanTestCase10() {
234
235 final double xm = 10000;
236 final double ym = 0;
237 final double sigmaX = 1000;
238 final double sigmaY = 10000;
239 final double radius = 10;
240
241
242 final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
243
244
245 Assertions.assertEquals(9.656e-28, result.getValue(), 1e-31);
246 }
247
248 @Test
249 @DisplayName("Chan test case 11")
250 void ChanTestCase11() {
251
252 final double xm = 0;
253 final double ym = 5000;
254 final double sigmaX = 1000;
255 final double sigmaY = 3000;
256 final double radius = 50;
257
258
259 final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
260
261
262 Assertions.assertEquals(1.039e-4, result.getValue(), 1e-7);
263 }
264
265 @Test
266 @DisplayName("Chan test case 12")
267 void ChanTestCase12() {
268
269 final double xm = 5000;
270 final double ym = 0;
271 final double sigmaX = 1000;
272 final double sigmaY = 3000;
273 final double radius = 50;
274
275
276 final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
277
278
279 Assertions.assertEquals(1.564e-9, result.getValue(), 1e-12);
280 }
281
282 @Test
283 @DisplayName("CSM test case 1")
284 void CsmTestCase1() {
285
286 final double xm = 84.875546;
287 final double ym = 60.583685;
288 final double sigmaX = 57.918666;
289 final double sigmaY = 152.8814468;
290 final double radius = 10.3;
291
292
293 final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
294
295
296 Assertions.assertEquals(1.9002e-3, result.getValue(), 1e-7);
297 }
298
299 @Test
300 @DisplayName("CSM test case 2")
301 void CsmTestCase2() {
302
303 final double xm = -81.618369;
304 final double ym = 115.055899;
305 final double sigmaX = 15.988242;
306 final double sigmaY = 5756.840725;
307 final double radius = 1.3;
308
309
310 final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
311
312
313 Assertions.assertEquals(2.0553e-11, result.getValue(), 1e-15);
314 }
315
316 @Test
317 @DisplayName("CSM test case 3")
318 void CsmTestCase3() {
319
320 final double xm = 102.177247;
321 final double ym = 693.405893;
322 final double sigmaX = 94.230921;
323 final double sigmaY = 643.409272;
324 final double radius = 5.3;
325
326
327 final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
328
329
330 Assertions.assertEquals(7.2004e-5, result.getValue(), 1e-9);
331 }
332
333 @Test
334 @DisplayName("CDM test case 1")
335 void CdmTestCase1() {
336
337 final double xm = -752.672701;
338 final double ym = 644.939441;
339 final double sigmaX = 445.859950;
340 final double sigmaY = 6095.858688;
341 final double radius = 3.5;
342
343
344 final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
345
346
347 Assertions.assertEquals(5.3904e-7, result.getValue(), 1e-11);
348 }
349
350 @Test
351 @DisplayName("CDM test case 2")
352 void CdmTestCase2() {
353
354 final double xm = -692.362272;
355 final double ym = 4475.456261;
356 final double sigmaX = 193.454603;
357 final double sigmaY = 562.027293;
358 final double radius = 13.2;
359
360
361 final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
362
363
364 Assertions.assertEquals(2.1652e-20, result.getValue(), 1e-24);
365 }
366
367 @Test
368 @DisplayName("Alfano test case 3")
369 void AlfanoTestCase3() {
370
371 final double xm = -3.8872073;
372 final double ym = 0.1591646;
373 final double sigmaX = 1.4101830;
374 final double sigmaY = 114.2585190;
375 final double radius = 15;
376
377
378 final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
379
380
381 Assertions.assertEquals(1.0038e-1, result.getValue(), 1e-5);
382 }
383
384 @Test
385 @DisplayName("Alfano test case 5")
386 void AlfanoTestCase5() {
387
388 final double xm = -1.2217895;
389 final double ym = 2.1230067;
390 final double sigmaX = 0.0373279;
391 final double sigmaY = 177.8109003;
392 final double radius = 10;
393
394
395 final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
396
397
398 Assertions.assertEquals(4.4510e-2, result.getValue(), 1e-6);
399 }
400
401 @Test
402 @DisplayName("Test that we can compute the Chan test case 01 from collision plane instance")
403 void testComputeChanTestCase01ProbabilityOfCollisionFromCollisionInstance() {
404
405 final ShortTermEncounter2DDefinition collision = Mockito.mock(ShortTermEncounter2DDefinition.class);
406
407
408 final double[][] covarianceMatrixData = { { 25 * 25, 0 }, { 0, 50 * 50 } };
409
410 Mockito.when(collision.computeProjectedAndDiagonalizedCombinedPositionalCovarianceMatrix())
411 .thenReturn(new BlockRealMatrix(covarianceMatrixData));
412
413 Mockito.when(collision.computeOtherPositionInRotatedCollisionPlane(Mockito.anyDouble()))
414 .thenReturn(new Vector2D(0, 10));
415
416 Mockito.when(collision.getCombinedRadius()).thenReturn(5.);
417
418 final ProbabilityOfCollision result = method.compute(collision);
419
420
421 Assertions.assertEquals(9.742e-3, result.getValue(), 1e-6);
422 }
423
424 @Test
425 @DisplayName("Test that we can compute expected probability of collision from primary and secondary collision object instance")
426 void testComputeExpectedProbabilityOfCollisionFromTimeAndCollisionObject() {
427
428
429
430 final AbsoluteDate timeOfClosestApproach = new AbsoluteDate();
431 final double mu = Constants.IERS2010_EARTH_MU;
432
433
434 final Frame primaryInertialFrame = FramesFactory.getEME2000();
435 final Orbit primary = new CartesianOrbit(
436 new PVCoordinates(new Vector3D(6778000, 0, 0), new Vector3D(0, 7668.631425, 0)), primaryInertialFrame,
437 timeOfClosestApproach, mu);
438 final RealMatrix primaryCovarianceMatrixInPrimaryRTN = new BlockRealMatrix(
439 new double[][] { { 100, 100, 100, 100, 100, 100 },
440 { 100, 100, 100, 100, 100, 100 },
441 { 100, 100, 200, 100, 100, 100 },
442 { 100, 100, 100, 100, 100, 100 },
443 { 100, 100, 100, 100, 100, 100 },
444 { 100, 100, 100, 100, 100, 100 } });
445 final StateCovariance primaryCovariance = new StateCovariance(primaryCovarianceMatrixInPrimaryRTN,
446 timeOfClosestApproach, LOFType.QSW_INERTIAL);
447 final double primaryRadius = 8;
448
449
450 final Frame secondaryInertialFrame = FramesFactory.getEME2000();
451 final Orbit secondary = new CartesianOrbit(
452 new PVCoordinates(new Vector3D(6778000 + 1, 0, 0), new Vector3D(0, 0, 7668.631425)),
453 secondaryInertialFrame, timeOfClosestApproach, mu);
454 final RealMatrix secondaryCovarianceMatrixInSecondaryRTN = new BlockRealMatrix(
455 new double[][] { { 100, 100, 100, 100, 100, 100 },
456 { 100, 100, 100, 100, 100, 100 },
457 { 100, 100, 200, 100, 100, 100 },
458 { 100, 100, 100, 100, 100, 100 },
459 { 100, 100, 100, 100, 100, 100 },
460 { 100, 100, 100, 100, 100, 100 } });
461 final StateCovariance secondaryCovariance = new StateCovariance(secondaryCovarianceMatrixInSecondaryRTN,
462 timeOfClosestApproach, LOFType.QSW_INERTIAL);
463 final double secondaryRadius = 2;
464
465
466 final ProbabilityOfCollision result = method.compute(primary, primaryCovariance, primaryRadius, secondary,
467 secondaryCovariance, secondaryRadius);
468
469
470 Assertions.assertEquals(0.21464810889751232, result.getValue(), 1e-17);
471 }
472
473 @Disabled("Statistics on alfano methods with armellin data : Values found in the data are" +
474 " said to be computed using the same Alfano method as implemented here but found results are order of" +
475 " magnitudes different from what is expected (even when computed with other methods). It is suspected that" +
476 " these values were computed using a maximum probability method similar to Alfriend1999 ")
477 @Test
478 @DisplayName("Test this method on Armellin's data and compare statistics about the probability of collision computed with the Alfano method")
479 void testCompareStatisticsAboutAlfanoProbabilityOfCollisionWithArmellinData() throws IOException {
480
481 final List<ArmellinDataRow> armellinDataRowList = ArmellinDataLoader.load();
482
483
484 final DescriptiveStatistics statistics =
485 ArmellinStatistics.getAlfanoProbabilityOfCollisionRelativeDifferenceStatistics(
486 armellinDataRowList);
487
488
489 Assertions.assertTrue(statistics.getMean() <= 1e-9);
490 Assertions.assertTrue(statistics.getStandardDeviation() <= 1e-9);
491 }
492
493 @Test
494 @DisplayName("Chan test case 01 Field version")
495 void ChanTestCase01Field() {
496
497 final Binary64 xm = new Binary64(0);
498 final Binary64 ym = new Binary64(10);
499 final Binary64 sigmaX = new Binary64(25);
500 final Binary64 sigmaY = new Binary64(50);
501 final Binary64 radius = new Binary64(5);
502
503
504 final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
505
506
507 Assertions.assertEquals(9.742e-3, result.getValue().getReal(), 1e-6);
508 }
509
510 @Test
511 @DisplayName("Chan test case 02 Field version")
512 void ChanTestCase02Field() {
513
514 final Binary64 xm = new Binary64(10);
515 final Binary64 ym = new Binary64(0);
516 final Binary64 sigmaX = new Binary64(25);
517 final Binary64 sigmaY = new Binary64(50);
518 final Binary64 radius = new Binary64(5);
519
520
521 final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
522
523
524 Assertions.assertEquals(9.181e-3, result.getValue().getReal(), 1e-6);
525 }
526
527 @Test
528 @DisplayName("Chan test case 03 Field version")
529 void ChanTestCase03Field() {
530
531 final Binary64 xm = new Binary64(0);
532 final Binary64 ym = new Binary64(10);
533 final Binary64 sigmaX = new Binary64(25);
534 final Binary64 sigmaY = new Binary64(75);
535 final Binary64 radius = new Binary64(5);
536
537
538 final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
539
540
541 Assertions.assertEquals(6.571e-3, result.getValue().getReal(), 1e-6);
542 }
543
544 @Test
545 @DisplayName("Chan test case 04 Field version")
546 void ChanTestCase04Field() {
547
548 final Binary64 xm = new Binary64(10);
549 final Binary64 ym = new Binary64(0);
550 final Binary64 sigmaX = new Binary64(25);
551 final Binary64 sigmaY = new Binary64(75);
552 final Binary64 radius = new Binary64(5);
553
554
555 final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
556
557
558 Assertions.assertEquals(6.125e-3, result.getValue().getReal(), 1e-6);
559 }
560
561 @Test
562 @DisplayName("Chan test case 05 Field version")
563 void ChanTestCase05Field() {
564
565 final Binary64 xm = new Binary64(0);
566 final Binary64 ym = new Binary64(1000);
567 final Binary64 sigmaX = new Binary64(1000);
568 final Binary64 sigmaY = new Binary64(3000);
569 final Binary64 radius = new Binary64(10);
570
571
572 final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
573
574
575 Assertions.assertEquals(1.577e-5, result.getValue().getReal(), 1e-8);
576 }
577
578 @Test
579 @DisplayName("Chan test case 06 Field version")
580 void ChanTestCase06Field() {
581
582 final Binary64 xm = new Binary64(1000);
583 final Binary64 ym = new Binary64(0);
584 final Binary64 sigmaX = new Binary64(1000);
585 final Binary64 sigmaY = new Binary64(3000);
586 final Binary64 radius = new Binary64(10);
587
588
589 final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
590
591
592 Assertions.assertEquals(1.011e-5, result.getValue().getReal(), 1e-8);
593 }
594
595 @Test
596 @DisplayName("Chan test case 07 Field version")
597 void ChanTestCase07Field() {
598
599 final Binary64 xm = new Binary64(0);
600 final Binary64 ym = new Binary64(10000);
601 final Binary64 sigmaX = new Binary64(1000);
602 final Binary64 sigmaY = new Binary64(3000);
603 final Binary64 radius = new Binary64(10);
604
605
606 final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
607
608
609 Assertions.assertEquals(6.443e-8, result.getValue().getReal(), 1e-11);
610 }
611
612 @Test
613 @DisplayName("Chan test case 08 Field version")
614 void ChanTestCase08Field() {
615
616 final Binary64 xm = new Binary64(10000);
617 final Binary64 ym = new Binary64(0);
618 final Binary64 sigmaX = new Binary64(1000);
619 final Binary64 sigmaY = new Binary64(3000);
620 final Binary64 radius = new Binary64(10);
621
622
623 final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
624
625
626 Assertions.assertEquals(3.219e-27, result.getValue().getReal(), 1e-30);
627 }
628
629 @Test
630 @DisplayName("Chan test case 09 Field version")
631 void ChanTestCase09Field() {
632
633 final Binary64 xm = new Binary64(0);
634 final Binary64 ym = new Binary64(10000);
635 final Binary64 sigmaX = new Binary64(1000);
636 final Binary64 sigmaY = new Binary64(10000);
637 final Binary64 radius = new Binary64(10);
638
639
640 final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
641
642
643 Assertions.assertEquals(3.033e-6, result.getValue().getReal(), 1e-9);
644 }
645
646 @Test
647 @DisplayName("Chan test case 10 Field version")
648 void ChanTestCase10Field() {
649
650 final Binary64 xm = new Binary64(10000);
651 final Binary64 ym = new Binary64(0);
652 final Binary64 sigmaX = new Binary64(1000);
653 final Binary64 sigmaY = new Binary64(10000);
654 final Binary64 radius = new Binary64(10);
655
656
657 final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
658
659
660 Assertions.assertEquals(9.656e-28, result.getValue().getReal(), 1e-31);
661 }
662
663 @Test
664 @DisplayName("Chan test case 11 Field version")
665 void ChanTestCase11Field() {
666
667 final Binary64 xm = new Binary64(0);
668 final Binary64 ym = new Binary64(5000);
669 final Binary64 sigmaX = new Binary64(1000);
670 final Binary64 sigmaY = new Binary64(3000);
671 final Binary64 radius = new Binary64(50);
672
673
674 final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
675
676
677 Assertions.assertEquals(1.039e-4, result.getValue().getReal(), 1e-7);
678 }
679
680 @Test
681 @DisplayName("Chan test case 12 Field version")
682 void ChanTestCase12Field() {
683
684 final Binary64 xm = new Binary64(5000);
685 final Binary64 ym = new Binary64(0);
686 final Binary64 sigmaX = new Binary64(1000);
687 final Binary64 sigmaY = new Binary64(3000);
688 final Binary64 radius = new Binary64(50);
689
690
691 final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
692
693
694 Assertions.assertEquals(1.564e-9, result.getValue().getReal(), 1e-12);
695 }
696
697 @Test
698 @DisplayName("CSM test case 1 Field version")
699 void CsmTestCase1Field() {
700
701 final Binary64 xm = new Binary64(84.875546);
702 final Binary64 ym = new Binary64(60.583685);
703 final Binary64 sigmaX = new Binary64(57.918666);
704 final Binary64 sigmaY = new Binary64(152.8814468);
705 final Binary64 radius = new Binary64(10.3);
706
707
708 final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
709
710
711 Assertions.assertEquals(1.9002e-3, result.getValue().getReal(), 1e-7);
712 }
713
714 @Test
715 @DisplayName("CSM test case 2 Field version")
716 void CsmTestCase2Field() {
717
718 final Binary64 xm = new Binary64(-81.618369);
719 final Binary64 ym = new Binary64(115.055899);
720 final Binary64 sigmaX = new Binary64(15.988242);
721 final Binary64 sigmaY = new Binary64(5756.840725);
722 final Binary64 radius = new Binary64(1.3);
723
724
725 final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
726
727
728 Assertions.assertEquals(2.0553e-11, result.getValue().getReal(), 1e-15);
729 }
730
731 @Test
732 @DisplayName("CSM test case 3 Field version")
733 void CsmTestCase3Field() {
734
735 final Binary64 xm = new Binary64(102.177247);
736 final Binary64 ym = new Binary64(693.405893);
737 final Binary64 sigmaX = new Binary64(94.230921);
738 final Binary64 sigmaY = new Binary64(643.409272);
739 final Binary64 radius = new Binary64(5.3);
740
741
742 final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
743
744
745 Assertions.assertEquals(7.2004e-5, result.getValue().getReal(), 1e-9);
746 }
747
748 @Test
749 @DisplayName("CDM test case 1 Field version")
750 void CdmTestCase1Field() {
751
752 final Binary64 xm = new Binary64(-752.672701);
753 final Binary64 ym = new Binary64(644.939441);
754 final Binary64 sigmaX = new Binary64(445.859950);
755 final Binary64 sigmaY = new Binary64(6095.858688);
756 final Binary64 radius = new Binary64(3.5);
757
758
759 final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
760
761
762 Assertions.assertEquals(5.3904e-7, result.getValue().getReal(), 1e-11);
763 }
764
765 @Test
766 @DisplayName("CDM test case 2 Field version")
767 void CdmTestCase2Field() {
768
769 final Binary64 xm = new Binary64(-692.362272);
770 final Binary64 ym = new Binary64(4475.456261);
771 final Binary64 sigmaX = new Binary64(193.454603);
772 final Binary64 sigmaY = new Binary64(562.027293);
773 final Binary64 radius = new Binary64(13.2);
774
775
776 final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
777
778
779 Assertions.assertEquals(2.1652e-20, result.getValue().getReal(), 1e-24);
780 }
781
782 @Test
783 @DisplayName("Alfano test case 3 Field version")
784 void AlfanoTestCase3Field() {
785
786 final Binary64 xm = new Binary64(-3.8872073);
787 final Binary64 ym = new Binary64(0.1591646);
788 final Binary64 sigmaX = new Binary64(1.4101830);
789 final Binary64 sigmaY = new Binary64(114.2585190);
790 final Binary64 radius = new Binary64(15);
791
792
793 final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
794
795
796 Assertions.assertEquals(1.0038e-1, result.getValue().getReal(), 1e-5);
797 }
798
799 @Test
800 @DisplayName("Alfano test case 5 Field version")
801 void AlfanoTestCase5Field() {
802
803 final Binary64 xm = new Binary64(-1.2217895);
804 final Binary64 ym = new Binary64(2.1230067);
805 final Binary64 sigmaX = new Binary64(0.0373279);
806 final Binary64 sigmaY = new Binary64(177.8109003);
807 final Binary64 radius = new Binary64(10);
808
809
810 final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
811
812
813 Assertions.assertEquals(4.4510e-2, result.getValue().getReal(), 1e-6);
814 }
815
816 @Test
817 @DisplayName("Test impact on the probability of collision of a slight difference in combined radius in Chan test case 4")
818 void testReturnExpectedValueWhenIntroducingSmallChangeOnCombinedRadius() {
819
820 final DSFactory factory = new DSFactory(1, 10);
821
822 final double xmNominal = 10;
823 final double ymNominal = 0;
824 final double sigmaXNominal = 25;
825 final double sigmaYNominal = 75;
826 final double radiusNominal = 5;
827 final double dRadius = 2.5;
828
829 final DerivativeStructure xm = factory.constant(xmNominal);
830 final DerivativeStructure ym = factory.constant(ymNominal);
831 final DerivativeStructure sigmaX = factory.constant(sigmaXNominal);
832 final DerivativeStructure sigmaY = factory.constant(sigmaYNominal);
833 final DerivativeStructure radius = factory.variable(0, radiusNominal);
834
835
836 final FieldProbabilityOfCollision<DerivativeStructure> resultNominal =
837 method.compute(xm, ym, sigmaX, sigmaY, radius);
838 final double taylorResult = resultNominal.getValue().taylor(dRadius);
839 final double exactResult =
840 method.compute(xmNominal, ymNominal, sigmaXNominal, sigmaYNominal, radiusNominal + dRadius).getValue();
841
842
843 Assertions.assertEquals(6.1e-3, resultNominal.getValue().getReal(), 1e-4);
844 Assertions.assertEquals(exactResult, taylorResult, 1e-15);
845 }
846
847 @Test
848 @DisplayName("Test that we can compute the Chan test case 01 from collision plane instance Field version")
849 void testComputeChanTestCase01ProbabilityOfCollisionFromCollisionInstanceField() {
850
851 @SuppressWarnings("unchecked")
852 final FieldShortTermEncounter2DDefinition<Binary64>
853 collision = Mockito.mock(FieldShortTermEncounter2DDefinition.class);
854
855
856 final Binary64[][] covarianceMatrixData = {
857 { new Binary64(25 * 25), new Binary64(0) },
858 { new Binary64(0), new Binary64(50 * 50) } };
859 Mockito.when(collision.computeProjectedAndDiagonalizedCombinedPositionalCovarianceMatrix())
860 .thenReturn(new BlockFieldMatrix<>(covarianceMatrixData));
861
862 final FieldVector2D<Binary64> otherPositionInRotatedCollisionPlane =
863 new FieldVector2D<>(new Binary64[] { new Binary64(0), new Binary64(10) });
864 Mockito.when(collision.computeOtherPositionInRotatedCollisionPlane(Mockito.anyDouble()))
865 .thenReturn(otherPositionInRotatedCollisionPlane);
866
867 Mockito.when(collision.getCombinedRadius()).thenReturn(new Binary64(5.));
868
869 final FieldProbabilityOfCollision<Binary64> result = method.compute(collision);
870
871
872 Assertions.assertEquals(9.742e-3, result.getValue().getReal(), 1e-6);
873 }
874
875 @Test
876 @DisplayName("Test that we can compute expected probability of collision from primary and secondary collision object instance Field version")
877 void testComputeExpectedProbabilityOfCollisionFromTimeAndCollisionObjectField() {
878
879
880
881 final Field<Binary64> field = Binary64Field.getInstance();
882 final FieldAbsoluteDate<Binary64> timeOfClosestApproach = new FieldAbsoluteDate<>(field);
883 final Binary64 mu = new Binary64(Constants.IERS2010_EARTH_MU);
884
885
886 final Frame primaryInertialFrame = FramesFactory.getEME2000();
887
888 final FieldOrbit<Binary64> primary = new FieldCartesianOrbit<>(
889 new FieldPVCoordinates<>(new FieldVector3D<>(new Binary64(6778000), new Binary64(0), new Binary64(0)),
890 new FieldVector3D<>(new Binary64(0), new Binary64(7668.631425), new Binary64(0))),
891 primaryInertialFrame, timeOfClosestApproach, mu);
892
893 final FieldMatrix<Binary64> primaryCovarianceMatrixInPrimaryRTN = new BlockFieldMatrix<>(
894 new Binary64[][] {
895 { new Binary64(100), new Binary64(100), new Binary64(100), new Binary64(100), new Binary64(100),
896 new Binary64(100) },
897 { new Binary64(100), new Binary64(100), new Binary64(100), new Binary64(100), new Binary64(100),
898 new Binary64(100) },
899 { new Binary64(100), new Binary64(100), new Binary64(200), new Binary64(100), new Binary64(100),
900 new Binary64(100) },
901 { new Binary64(100), new Binary64(100), new Binary64(100), new Binary64(100), new Binary64(100),
902 new Binary64(100) },
903 { new Binary64(100), new Binary64(100), new Binary64(100), new Binary64(100), new Binary64(100),
904 new Binary64(100) },
905 { new Binary64(100), new Binary64(100), new Binary64(100), new Binary64(100), new Binary64(100),
906 new Binary64(100) } });
907
908 final FieldStateCovariance<Binary64> primaryCovariance =
909 new FieldStateCovariance<>(primaryCovarianceMatrixInPrimaryRTN,
910 timeOfClosestApproach, LOFType.QSW_INERTIAL);
911 final Binary64 primaryRadius = new Binary64(8);
912
913
914 final Frame secondaryInertialFrame = FramesFactory.getEME2000();
915
916 final FieldOrbit<Binary64> secondary = new FieldCartesianOrbit<>(
917 new FieldPVCoordinates<>(
918 new FieldVector3D<>(new Binary64(6778000 + 1), new Binary64(0.), new Binary64(0.)),
919 new FieldVector3D<>(new Binary64(0.), new Binary64(0.), new Binary64(7668.631425))),
920 secondaryInertialFrame, timeOfClosestApproach, mu);
921
922 final FieldMatrix<Binary64> secondaryCovarianceMatrixInSecondaryRTN = new BlockFieldMatrix<>(
923 new Binary64[][] {
924 { new Binary64(100), new Binary64(100), new Binary64(100), new Binary64(100), new Binary64(100),
925 new Binary64(100) },
926 { new Binary64(100), new Binary64(100), new Binary64(100), new Binary64(100), new Binary64(100),
927 new Binary64(100) },
928 { new Binary64(100), new Binary64(100), new Binary64(200), new Binary64(100), new Binary64(100),
929 new Binary64(100) },
930 { new Binary64(100), new Binary64(100), new Binary64(100), new Binary64(100), new Binary64(100),
931 new Binary64(100) },
932 { new Binary64(100), new Binary64(100), new Binary64(100), new Binary64(100), new Binary64(100),
933 new Binary64(100) },
934 { new Binary64(100), new Binary64(100), new Binary64(100), new Binary64(100), new Binary64(100),
935 new Binary64(100) } });
936
937 final FieldStateCovariance<Binary64> secondaryCovariance =
938 new FieldStateCovariance<>(secondaryCovarianceMatrixInSecondaryRTN,
939 timeOfClosestApproach, LOFType.QSW_INERTIAL);
940
941 final Binary64 secondaryRadius = new Binary64(2);
942
943
944 final FieldProbabilityOfCollision<Binary64> result = method.compute(primary, primaryCovariance, primaryRadius,
945 secondary, secondaryCovariance, secondaryRadius);
946
947
948 Assertions.assertEquals(0.21464810889751232, result.getValue().getReal(), 1e-17);
949 }
950
951
952 @Test
953 @DisplayName("Alfano test case 1 (CDM) from NASA CARA")
954 void AlfanoCDMTestCase01() {
955
956 String cdmName = "AlfanoTestCase01.cdm";
957 final double combinedHbr = 15.;
958
959
960 final double expectedPc = 0.146749549;
961 final double tolerance = 1e-6;
962
963 computeAndCheckCollisionProbability(cdmName, combinedHbr, expectedPc, tolerance);
964 }
965
966 @Test
967 @DisplayName("Alfano test case 2 (CDM) from NASA CARA")
968 void AlfanoCDMTestCase02() {
969
970 String cdmName = "AlfanoTestCase02.cdm";
971 final double combinedHbr = 4.;
972
973
974 final double expectedPc = 0.006222267;
975 final double tolerance = 1e-6;
976
977 computeAndCheckCollisionProbability(cdmName, combinedHbr, expectedPc, tolerance);
978 }
979
980 @Test
981 @DisplayName("Alfano test case 3 (CDM) from NASA CARA")
982 void AlfanoCDMTestCase03() {
983
984 String cdmName = "AlfanoTestCase03.cdm";
985 final double combinedHbr = 15.;
986
987
988 final double expectedPc = 0.100351176;
989 final double tolerance = 1e-6;
990
991 computeAndCheckCollisionProbability(cdmName, combinedHbr, expectedPc, tolerance);
992 }
993
994 @Test
995 @DisplayName("Alfano test case 4 (CDM) from NASA CARA")
996 void AlfanoCDMTestCase04() {
997
998 String cdmName = "AlfanoTestCase04.cdm";
999 final double combinedHbr = 15.;
1000
1001
1002 final double expectedPc = 0.049323406;
1003 final double tolerance = 1e-5;
1004
1005 computeAndCheckCollisionProbability(cdmName, combinedHbr, expectedPc, tolerance);
1006 }
1007
1008 @Test
1009 @DisplayName("Alfano test case 5 (CDM) from NASA CARA")
1010 void AlfanoCDMTestCase05() {
1011
1012 String cdmName = "AlfanoTestCase05.cdm";
1013 final double combinedHbr = 10.;
1014
1015
1016 final double expectedPc = 0.044487386;
1017 final double tolerance = 1e-5;
1018
1019 computeAndCheckCollisionProbability(cdmName, combinedHbr, expectedPc, tolerance);
1020 }
1021
1022 @Test
1023 @DisplayName("Alfano test case 6 (CDM) from NASA CARA")
1024 void AlfanoCDMTestCase06() {
1025
1026 String cdmName = "AlfanoTestCase06.cdm";
1027 final double combinedHbr = 10.;
1028
1029
1030 final double expectedPc = 0.004335455;
1031 final double tolerance = 1e-8;
1032
1033 computeAndCheckCollisionProbability(cdmName, combinedHbr, expectedPc, tolerance);
1034 }
1035
1036 @Test
1037 @DisplayName("Alfano test case 7 (CDM) from NASA CARA")
1038 void AlfanoCDMTestCase07() {
1039
1040 String cdmName = "AlfanoTestCase07.cdm";
1041 final double combinedHbr = 10.;
1042
1043
1044 final double expectedPc = 0.000158147;
1045 final double tolerance = 1e-9;
1046
1047 computeAndCheckCollisionProbability(cdmName, combinedHbr, expectedPc, tolerance);
1048 }
1049
1050 @Test
1051 @DisplayName("Alfano test case 8 (CDM) from NASA CARA")
1052 void AlfanoCDMTestCase08() {
1053
1054 String cdmName = "AlfanoTestCase08.cdm";
1055 final double combinedHbr = 4.;
1056
1057
1058 final double expectedPc = 0.036948008;
1059 final double tolerance = 1e-5;
1060
1061 computeAndCheckCollisionProbability(cdmName, combinedHbr, expectedPc, tolerance);
1062 }
1063
1064 @Test
1065 @DisplayName("Alfano test case 9 (CDM) from NASA CARA")
1066 void AlfanoCDMTestCase09() {
1067
1068 String cdmName = "AlfanoTestCase09.cdm";
1069 final double combinedHbr = 6.;
1070
1071
1072 final double expectedPc = 0.290146291;
1073 final double tolerance = 2e-5;
1074
1075 computeAndCheckCollisionProbability(cdmName, combinedHbr, expectedPc, tolerance);
1076 }
1077
1078
1079
1080 @Test
1081 @DisplayName("Alfano test case 11 (CDM) from NASA CARA")
1082 void AlfanoCDMTestCase11() {
1083
1084 String cdmName = "AlfanoTestCase11.cdm";
1085 final double combinedHbr = 4.;
1086
1087
1088 final double expectedPc = 0.002672026;
1089 final double tolerance = 1e-7;
1090
1091 computeAndCheckCollisionProbability(cdmName, combinedHbr, expectedPc, tolerance);
1092 }
1093
1094 private void computeAndCheckCollisionProbability(
1095 final String cdmName,
1096 final double combinedHbr,
1097 final double expected,
1098 final double tolerance) {
1099
1100
1101 final String cdmPath = "/ccsds/cdm/" + cdmName;
1102 final DataSource data = new DataSource(cdmPath, () -> getClass().getResourceAsStream(cdmPath));
1103 final Cdm cdm = new ParserBuilder().buildCdmParser().parseMessage(data);
1104
1105
1106 final ProbabilityOfCollision result = method.compute(cdm, combinedHbr);
1107
1108
1109 Assertions.assertEquals(expected, result.getValue(), tolerance);
1110 }
1111 }