1   /* Copyright 2002-2023 CS GROUP
2    * Licensed to CS GROUP (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *   http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.orekit.ssa.collision.shorttermencounter.probability.twod;
18  
19  import org.hipparchus.analysis.differentiation.DSFactory;
20  import org.hipparchus.analysis.differentiation.DerivativeStructure;
21  import org.hipparchus.util.Binary64;
22  import org.junit.jupiter.api.Assertions;
23  import org.junit.jupiter.api.BeforeAll;
24  import org.junit.jupiter.api.DisplayName;
25  import org.junit.jupiter.api.Test;
26  import org.orekit.Utils;
27  import org.orekit.data.DataSource;
28  import org.orekit.files.ccsds.ndm.ParserBuilder;
29  import org.orekit.files.ccsds.ndm.cdm.Cdm;
30  import org.orekit.ssa.metrics.FieldProbabilityOfCollision;
31  import org.orekit.ssa.metrics.ProbabilityOfCollision;
32  
33  class Laas2015Test {
34  
35      /** Systems analysis and architecture laboratory's method to compute probability of collision. */
36      private final ShortTermEncounter2DPOCMethod method = new Laas2015();
37  
38      @BeforeAll
39      static void initializeOrekitData() {
40          Utils.setDataRoot("regular-data");
41      }
42  
43      @Test
44      @DisplayName("Chan test case 01")
45      void ChanTestCase01() {
46          // GIVEN
47          final double xm     = 0;
48          final double ym     = 10;
49          final double sigmaX = 25;
50          final double sigmaY = 50;
51          final double radius = 5;
52  
53          // WHEN
54          final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
55  
56          // THEN
57          Assertions.assertEquals(9.742e-3, result.getValue(), 1e-6);
58          Assertions.assertEquals(9.704e-3, result.getLowerLimit(), 1e-6);
59          Assertions.assertEquals(9.742e-3, result.getUpperLimit(), 1e-6);
60      }
61  
62      @Test
63      @DisplayName("Chan test case 02")
64      void ChanTestCase02() {
65          // GIVEN
66          final double xm     = 10;
67          final double ym     = 0;
68          final double sigmaX = 25;
69          final double sigmaY = 50;
70          final double radius = 5;
71  
72          // WHEN
73          final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
74  
75          // THEN
76          Assertions.assertEquals(9.181e-3, result.getValue(), 1e-6);
77          Assertions.assertEquals(9.139e-3, result.getLowerLimit(), 1e-6);
78          Assertions.assertEquals(9.182e-3, result.getUpperLimit(), 1e-6);
79      }
80  
81      @Test
82      @DisplayName("Chan test case 03")
83      void ChanTestCase03() {
84          // GIVEN
85          final double xm     = 0;
86          final double ym     = 10;
87          final double sigmaX = 25;
88          final double sigmaY = 75;
89          final double radius = 5;
90  
91          // WHEN
92          final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
93  
94          // THEN
95          Assertions.assertEquals(6.571e-3, result.getValue(), 1e-6);
96          Assertions.assertEquals(6.542e-3, result.getLowerLimit(), 1e-6);
97          Assertions.assertEquals(6.572e-3, result.getUpperLimit(), 1e-6);
98      }
99  
100     @Test
101     @DisplayName("Chan test case 04")
102     void ChanTestCase04() {
103         // GIVEN
104         final double xm     = 10;
105         final double ym     = 0;
106         final double sigmaX = 25;
107         final double sigmaY = 75;
108         final double radius = 5;
109 
110         // WHEN
111         final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
112 
113         // THEN
114         Assertions.assertEquals(6.125e-3, result.getValue(), 1e-6);
115         Assertions.assertEquals(6.09e-3, result.getLowerLimit(), 1e-5);
116         Assertions.assertEquals(6.13e-3, result.getUpperLimit(), 1e-5);
117     }
118 
119     @Test
120     @DisplayName("Chan test case 05")
121     void ChanTestCase05() {
122         // GIVEN
123         final double xm     = 0;
124         final double ym     = 1000;
125         final double sigmaX = 1000;
126         final double sigmaY = 3000;
127         final double radius = 10;
128 
129         // WHEN
130         final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
131 
132         // THEN
133         Assertions.assertEquals(1.577e-5, result.getValue(), 1e-8);
134         Assertions.assertEquals(1.576561e-5, result.getLowerLimit(), 1e-9);
135         Assertions.assertEquals(1.576576e-5, result.getUpperLimit(), 1e-9);
136     }
137 
138     @Test
139     @DisplayName("Chan test case 06")
140     void ChanTestCase06() {
141         // GIVEN
142         final double xm     = 1000;
143         final double ym     = 0;
144         final double sigmaX = 1000;
145         final double sigmaY = 3000;
146         final double radius = 10;
147 
148         // WHEN
149         final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
150 
151         // THEN
152         Assertions.assertEquals(1.011e-5, result.getValue(), 1e-8);
153         Assertions.assertEquals(1.010860e-5, result.getLowerLimit(), 1e-11);
154         Assertions.assertEquals(1.010883e-5, result.getUpperLimit(), 1e-11);
155     }
156 
157     @Test
158     @DisplayName("Chan test case 07")
159     void ChanTestCase07() {
160         // GIVEN
161         final double xm     = 0;
162         final double ym     = 10000;
163         final double sigmaX = 1000;
164         final double sigmaY = 3000;
165         final double radius = 10;
166 
167         // WHEN
168         final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
169 
170         // THEN
171         Assertions.assertEquals(6.443e-8, result.getValue(), 1e-11);
172         Assertions.assertEquals(6.44304e-8, result.getLowerLimit(), 1e-13);
173         Assertions.assertEquals(6.44321e-8, result.getUpperLimit(), 1e-13);
174     }
175 
176     @Test
177     @DisplayName("Chan test case 08")
178     void ChanTestCase08() {
179         // GIVEN
180         final double xm     = 10000;
181         final double ym     = 0;
182         final double sigmaX = 1000;
183         final double sigmaY = 3000;
184         final double radius = 10;
185 
186         // WHEN
187         final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
188 
189         // THEN
190         Assertions.assertEquals(3.219e-27, result.getValue(), 1e-30);
191         Assertions.assertEquals(3.2145e-27, result.getLowerLimit(), 1e-31);
192         Assertions.assertEquals(3.2186e-27, result.getUpperLimit(), 1e-31);
193     }
194 
195     @Test
196     @DisplayName("Chan test case 09")
197     void ChanTestCase09() {
198         // GIVEN
199         final double xm     = 0;
200         final double ym     = 10000;
201         final double sigmaX = 1000;
202         final double sigmaY = 10000;
203         final double radius = 10;
204 
205         // WHEN
206         final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
207 
208         // THEN
209         Assertions.assertEquals(3.033e-6, result.getValue(), 1e-9);
210         Assertions.assertEquals(3.03258e-6, result.getLowerLimit(), 1e-11);
211         Assertions.assertEquals(3.03261e-6, result.getUpperLimit(), 1e-11);
212     }
213 
214     @Test
215     @DisplayName("Chan test case 10")
216     void ChanTestCase10() {
217         // GIVEN
218         final double xm     = 10000;
219         final double ym     = 0;
220         final double sigmaX = 1000;
221         final double sigmaY = 10000;
222         final double radius = 10;
223 
224         // WHEN
225         final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
226 
227         // THEN
228         Assertions.assertEquals(9.656e-28, result.getValue(), 1e-31);
229         Assertions.assertEquals(9.643e-28, result.getLowerLimit(), 1e-31);
230         Assertions.assertEquals(9.656e-28, result.getUpperLimit(), 1e-31);
231     }
232 
233     @Test
234     @DisplayName("Chan test case 11")
235     void ChanTestCase11() {
236         // GIVEN
237         final double xm     = 0;
238         final double ym     = 5000;
239         final double sigmaX = 1000;
240         final double sigmaY = 3000;
241         final double radius = 50;
242 
243         // WHEN
244         final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
245 
246         // THEN
247         Assertions.assertEquals(1.039e-4, result.getValue(), 1e-7);
248         Assertions.assertEquals(1.03831e-4, result.getLowerLimit(), 1e-9);
249         Assertions.assertEquals(1.03871e-4, result.getUpperLimit(), 1e-9);
250     }
251 
252     @Test
253     @DisplayName("Chan test case 12")
254     void ChanTestCase12() {
255         // GIVEN
256         final double xm     = 5000;
257         final double ym     = 0;
258         final double sigmaX = 1000;
259         final double sigmaY = 3000;
260         final double radius = 50;
261 
262         // WHEN
263         final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
264 
265         // THEN
266         Assertions.assertEquals(1.564e-9, result.getValue(), 1e-12);
267         Assertions.assertEquals(1.552e-9, result.getLowerLimit(), 1e-12);
268         Assertions.assertEquals(1.565e-9, result.getUpperLimit(), 1e-12);
269     }
270 
271     @Test
272     @DisplayName("CSM test case 1")
273     void CsmTestCase1() {
274         // GIVEN
275         final double xm     = 84.875546;
276         final double ym     = 60.583685;
277         final double sigmaX = 57.918666;
278         final double sigmaY = 152.8814468;
279         final double radius = 10.3;
280 
281         // WHEN
282         final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
283 
284         // THEN
285         Assertions.assertEquals(1.9002e-3, result.getValue(), 1e-7);
286         Assertions.assertEquals(1.878e-3, result.getLowerLimit(), 1e-6);
287         Assertions.assertEquals(1.900e-3, result.getUpperLimit(), 1e-6);
288     }
289 
290     @Test
291     @DisplayName("CSM test case 2")
292     void CsmTestCase2() {
293         // GIVEN
294         final double xm     = -81.618369;
295         final double ym     = 115.055899;
296         final double sigmaX = 15.988242;
297         final double sigmaY = 5756.840725;
298         final double radius = 1.3;
299 
300         // WHEN
301         final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
302 
303         // THEN
304         Assertions.assertEquals(2.0553e-11, result.getValue(), 1e-15);
305         Assertions.assertEquals(2.0101e-11, result.getLowerLimit(), 1e-15);
306         Assertions.assertEquals(2.0557e-11, result.getUpperLimit(), 1e-15);
307     }
308 
309     @Test
310     @DisplayName("CSM test case 3")
311     void CsmTestCase3() {
312         // GIVEN
313         final double xm     = 102.177247;
314         final double ym     = 693.405893;
315         final double sigmaX = 94.230921;
316         final double sigmaY = 643.409272;
317         final double radius = 5.3;
318 
319         // WHEN
320         final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
321 
322         // THEN
323         Assertions.assertEquals(7.2003e-5, result.getValue(), 1e-9);
324         Assertions.assertEquals(7.194e-5, result.getLowerLimit(), 1e-8);
325         Assertions.assertEquals(7.200e-5, result.getUpperLimit(), 1e-8);
326     }
327 
328     @Test
329     @DisplayName("CDM test case 1")
330     void CdmTestCase1() {
331         // GIVEN
332         final double xm     = -752.672701;
333         final double ym     = 644.939441;
334         final double sigmaX = 445.859950;
335         final double sigmaY = 6095.858688;
336         final double radius = 3.5;
337 
338         // WHEN
339         final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
340 
341         // THEN
342         Assertions.assertEquals(5.3904e-7, result.getValue(), 1e-11);
343         Assertions.assertEquals(5.3902e-7, result.getLowerLimit(), 1e-11);
344         Assertions.assertEquals(5.3904e-7, result.getUpperLimit(), 1e-11);
345     }
346 
347     @Test
348     @DisplayName("CDM test case 2")
349     void CdmTestCase2() {
350         // GIVEN
351         final double xm     = -692.362272;
352         final double ym     = 4475.456261;
353         final double sigmaX = 193.454603;
354         final double sigmaY = 562.027293;
355         final double radius = 13.2;
356 
357         // WHEN
358         final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
359 
360         // THEN
361         Assertions.assertEquals(2.2796e-20, result.getValue(), 1e-24);
362         Assertions.assertEquals(2.2517e-20, result.getLowerLimit(), 1e-24);
363         Assertions.assertEquals(2.2797e-20, result.getUpperLimit(), 1e-24);
364     }
365 
366     @Test
367     @DisplayName("Alfano test case 3")
368     void AlfanoTestCase3() {
369         // GIVEN
370         final double xm     = -3.8872073;
371         final double ym     = 0.1591646;
372         final double sigmaX = 1.4101830;
373         final double sigmaY = 114.2585190;
374         final double radius = 15;
375 
376         // WHEN
377         final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
378 
379         // THEN
380         Assertions.assertEquals(1.0038e-1, result.getValue(), 1e-5);
381     }
382 
383     /**
384      * Alfano test case 5.
385      * <p>
386      * Initial expected result of 4.4509e-2 has been slightly modified to 4.4507e-2 as this test case is numerically
387      * unstable. This has been validated with Romain SERRA.
388      * <p>
389      * Result obtained here is (slightly) different that what is obtained with the field version because the non-field uses
390      * the {@link org.hipparchus.util.MathArrays} linearCombination method to be more precise. However, the field equivalent
391      * does not exist.
392      */
393     @Test
394     @DisplayName("Alfano test case 5")
395     void AlfanoTestCase5() {
396         // GIVEN
397         final double xm     = -1.2217895;
398         final double ym     = 2.1230067;
399         final double sigmaX = 0.0373279;
400         final double sigmaY = 177.8109003;
401         final double radius = 10;
402 
403         // WHEN
404         final ProbabilityOfCollision result = method.compute(xm, ym, sigmaX, sigmaY, radius);
405 
406         // THEN
407         Assertions.assertEquals(4.4507e-2, result.getValue(), 1e-6);
408     }
409 
410     @Test
411     @DisplayName("Chan test case 04 with custom method constructor and relative boundaries close enough to the desired accuracy")
412     void ChanTestCase04WithCustomLaasMethod() {
413         // GIVEN
414         final ShortTermEncounter2DPOCMethod customMethod = new Laas2015(1e-4, 1000);
415 
416         final double xm     = 10;
417         final double ym     = 0;
418         final double sigmaX = 25;
419         final double sigmaY = 75;
420         final double radius = 5;
421 
422         // WHEN
423         final ProbabilityOfCollision result = customMethod.compute(xm, ym, sigmaX, sigmaY, radius);
424 
425         // THEN
426         Assertions.assertEquals(6.1e-3, result.getValue(), 1e-4);
427     }
428 
429     @Test
430     @DisplayName(
431             "Test probability from a real CDM. It has been mentioned in https://forum.orekit.org/t/collision-package-in-development/2638/16"
432                     + " and define an expected probability of collision of 0.004450713 which is higher than the value found using"
433                     + " the methods provided in this ssa package because the CSpOC use the projection of a square containing the"
434                     + " collision disk instead of the collision disk directly. Hence an overestimation of the probability of"
435                     + " collision. This test serves as a non regression test")
436     void testComputeProbabilityFromACdm() {
437 
438         // GIVEN
439         final String     cdmPath = "/ccsds/cdm/ION_SCV8_vs_STARLINK_1233.txt";
440         final DataSource data    = new DataSource(cdmPath, () -> getClass().getResourceAsStream(cdmPath));
441         final Cdm        cdm     = new ParserBuilder().buildCdmParser().parseMessage(data);
442 
443         // Radii taken from comments in the conjunction data message
444         final double primaryRadius   = 5;
445         final double secondaryRadius = 5;
446 
447         // WHEN
448         final ProbabilityOfCollision result = method.compute(cdm, primaryRadius, secondaryRadius);
449 
450         // THEN
451         Assertions.assertEquals(0.0034965176443840897, result.getValue(), 1e-18);
452     }
453 
454     @Test
455     @DisplayName("Chan test case 01 Field version")
456     void ChanTestCase01Field() {
457         // GIVEN
458         final Binary64 xm     = new Binary64(0);
459         final Binary64 ym     = new Binary64(10);
460         final Binary64 sigmaX = new Binary64(25);
461         final Binary64 sigmaY = new Binary64(50);
462         final Binary64 radius = new Binary64(5);
463 
464         // WHEN
465         final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
466 
467         // THEN
468         Assertions.assertEquals(9.742e-3, result.getValue().getReal(), 1e-6);
469         Assertions.assertEquals(9.704e-3, result.getLowerLimit().getReal(), 1e-6);
470         Assertions.assertEquals(9.742e-3, result.getUpperLimit().getReal(), 1e-6);
471     }
472 
473     @Test
474     @DisplayName("Chan test case 02 Field version")
475     void ChanTestCase02Field() {
476         // GIVEN
477         final Binary64 xm     = new Binary64(10);
478         final Binary64 ym     = new Binary64(0);
479         final Binary64 sigmaX = new Binary64(25);
480         final Binary64 sigmaY = new Binary64(50);
481         final Binary64 radius = new Binary64(5);
482 
483         // WHEN
484         final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
485 
486         // THEN
487         Assertions.assertEquals(9.181e-3, result.getValue().getReal(), 1e-6);
488         Assertions.assertEquals(9.139e-3, result.getLowerLimit().getReal(), 1e-6);
489         Assertions.assertEquals(9.182e-3, result.getUpperLimit().getReal(), 1e-6);
490     }
491 
492     @Test
493     @DisplayName("Chan test case 03 Field version")
494     void ChanTestCase03Field() {
495         // GIVEN
496         final Binary64 xm     = new Binary64(0);
497         final Binary64 ym     = new Binary64(10);
498         final Binary64 sigmaX = new Binary64(25);
499         final Binary64 sigmaY = new Binary64(75);
500         final Binary64 radius = new Binary64(5);
501 
502         // WHEN
503         final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
504 
505         // THEN
506         Assertions.assertEquals(6.571e-3, result.getValue().getReal(), 1e-6);
507         Assertions.assertEquals(6.542e-3, result.getLowerLimit().getReal(), 1e-6);
508         Assertions.assertEquals(6.572e-3, result.getUpperLimit().getReal(), 1e-6);
509     }
510 
511     @Test
512     @DisplayName("Chan test case 04 Field version")
513     void ChanTestCase04Field() {
514         // GIVEN
515         final Binary64 xm     = new Binary64(10);
516         final Binary64 ym     = new Binary64(0);
517         final Binary64 sigmaX = new Binary64(25);
518         final Binary64 sigmaY = new Binary64(75);
519         final Binary64 radius = new Binary64(5);
520 
521         // WHEN
522         final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
523 
524         // THEN
525         Assertions.assertEquals(6.125e-3, result.getValue().getReal(), 1e-6);
526         Assertions.assertEquals(6.09e-3, result.getLowerLimit().getReal(), 1e-5);
527         Assertions.assertEquals(6.13e-3, result.getUpperLimit().getReal(), 1e-5);
528     }
529 
530     @Test
531     @DisplayName("Chan test case 05 Field version")
532     void ChanTestCase05Field() {
533         // GIVEN
534         final Binary64 xm     = new Binary64(0);
535         final Binary64 ym     = new Binary64(1000);
536         final Binary64 sigmaX = new Binary64(1000);
537         final Binary64 sigmaY = new Binary64(3000);
538         final Binary64 radius = new Binary64(10);
539 
540         // WHEN
541         final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
542 
543         // THEN
544         Assertions.assertEquals(1.577e-5, result.getValue().getReal(), 1e-8);
545         Assertions.assertEquals(1.576561e-5, result.getLowerLimit().getReal(), 1e-9);
546         Assertions.assertEquals(1.576576e-5, result.getUpperLimit().getReal(), 1e-9);
547     }
548 
549     @Test
550     @DisplayName("Chan test case 06 Field version")
551     void ChanTestCase06Field() {
552         // GIVEN
553         final Binary64 xm     = new Binary64(1000);
554         final Binary64 ym     = new Binary64(0);
555         final Binary64 sigmaX = new Binary64(1000);
556         final Binary64 sigmaY = new Binary64(3000);
557         final Binary64 radius = new Binary64(10);
558 
559         // WHEN
560         final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
561 
562         // THEN
563         Assertions.assertEquals(1.011e-5, result.getValue().getReal(), 1e-8);
564         Assertions.assertEquals(1.010860e-5, result.getLowerLimit().getReal(), 1e-11);
565         Assertions.assertEquals(1.010883e-5, result.getUpperLimit().getReal(), 1e-11);
566     }
567 
568     @Test
569     @DisplayName("Chan test case 07 Field version")
570     void ChanTestCase07Field() {
571         // GIVEN
572         final Binary64 xm     = new Binary64(0);
573         final Binary64 ym     = new Binary64(10000);
574         final Binary64 sigmaX = new Binary64(1000);
575         final Binary64 sigmaY = new Binary64(3000);
576         final Binary64 radius = new Binary64(10);
577 
578         // WHEN
579         final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
580 
581         // THEN
582         Assertions.assertEquals(6.443e-8, result.getValue().getReal(), 1e-11);
583         Assertions.assertEquals(6.44304e-8, result.getLowerLimit().getReal(), 1e-13);
584         Assertions.assertEquals(6.44321e-8, result.getUpperLimit().getReal(), 1e-13);
585     }
586 
587     @Test
588     @DisplayName("Chan test case 08 Field version")
589     void ChanTestCase08Field() {
590         // GIVEN
591         final Binary64 xm     = new Binary64(10000);
592         final Binary64 ym     = new Binary64(0);
593         final Binary64 sigmaX = new Binary64(1000);
594         final Binary64 sigmaY = new Binary64(3000);
595         final Binary64 radius = new Binary64(10);
596 
597         // WHEN
598         final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
599 
600         // THEN
601         Assertions.assertEquals(3.219e-27, result.getValue().getReal(), 1e-30);
602         Assertions.assertEquals(3.2145e-27, result.getLowerLimit().getReal(), 1e-31);
603         Assertions.assertEquals(3.2186e-27, result.getUpperLimit().getReal(), 1e-31);
604     }
605 
606     @Test
607     @DisplayName("Chan test case 09 Field version")
608     void ChanTestCase09Field() {
609         // GIVEN
610         final Binary64 xm     = new Binary64(0);
611         final Binary64 ym     = new Binary64(10000);
612         final Binary64 sigmaX = new Binary64(1000);
613         final Binary64 sigmaY = new Binary64(10000);
614         final Binary64 radius = new Binary64(10);
615 
616         // WHEN
617         final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
618 
619         // THEN
620         Assertions.assertEquals(3.033e-6, result.getValue().getReal(), 1e-9);
621         Assertions.assertEquals(3.03258e-6, result.getLowerLimit().getReal(), 1e-11);
622         Assertions.assertEquals(3.03261e-6, result.getUpperLimit().getReal(), 1e-11);
623     }
624 
625     @Test
626     @DisplayName("Chan test case 10 Field version")
627     void ChanTestCase10Field() {
628         // GIVEN
629         final Binary64 xm     = new Binary64(10000);
630         final Binary64 ym     = new Binary64(0);
631         final Binary64 sigmaX = new Binary64(1000);
632         final Binary64 sigmaY = new Binary64(10000);
633         final Binary64 radius = new Binary64(10);
634 
635         // WHEN
636         final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
637 
638         // THEN
639         Assertions.assertEquals(9.656e-28, result.getValue().getReal(), 1e-31);
640         Assertions.assertEquals(9.643e-28, result.getLowerLimit().getReal(), 1e-31);
641         Assertions.assertEquals(9.656e-28, result.getUpperLimit().getReal(), 1e-31);
642     }
643 
644     @Test
645     @DisplayName("Chan test case 11 Field version")
646     void ChanTestCase11Field() {
647         // GIVEN
648         final Binary64 xm     = new Binary64(0);
649         final Binary64 ym     = new Binary64(5000);
650         final Binary64 sigmaX = new Binary64(1000);
651         final Binary64 sigmaY = new Binary64(3000);
652         final Binary64 radius = new Binary64(50);
653 
654         // WHEN
655         final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
656 
657         // THEN
658         Assertions.assertEquals(1.039e-4, result.getValue().getReal(), 1e-7);
659         Assertions.assertEquals(1.03831e-4, result.getLowerLimit().getReal(), 1e-9);
660         Assertions.assertEquals(1.03871e-4, result.getUpperLimit().getReal(), 1e-9);
661     }
662 
663     @Test
664     @DisplayName("Chan test case 12 Field version")
665     void ChanTestCase12Field() {
666         // GIVEN
667         final Binary64 xm     = new Binary64(5000);
668         final Binary64 ym     = new Binary64(0);
669         final Binary64 sigmaX = new Binary64(1000);
670         final Binary64 sigmaY = new Binary64(3000);
671         final Binary64 radius = new Binary64(50);
672 
673         // WHEN
674         final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
675 
676         // THEN
677         Assertions.assertEquals(1.564e-9, result.getValue().getReal(), 1e-12);
678         Assertions.assertEquals(1.552e-9, result.getLowerLimit().getReal(), 1e-12);
679         Assertions.assertEquals(1.565e-9, result.getUpperLimit().getReal(), 1e-12);
680     }
681 
682     @Test
683     @DisplayName("CSM test case 1 Field version")
684     void CsmTestCase1Field() {
685         // GIVEN
686         final Binary64 xm     = new Binary64(84.875546);
687         final Binary64 ym     = new Binary64(60.583685);
688         final Binary64 sigmaX = new Binary64(57.918666);
689         final Binary64 sigmaY = new Binary64(152.8814468);
690         final Binary64 radius = new Binary64(10.3);
691 
692         // WHEN
693         final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
694 
695         // THEN
696         Assertions.assertEquals(1.9002e-3, result.getValue().getReal(), 1e-7);
697         Assertions.assertEquals(1.878e-3, result.getLowerLimit().getReal(), 1e-6);
698         Assertions.assertEquals(1.900e-3, result.getUpperLimit().getReal(), 1e-6);
699     }
700 
701     @Test
702     @DisplayName("CSM test case 2 Field version")
703     void CsmTestCase2Field() {
704         // GIVEN
705         final Binary64 xm     = new Binary64(-81.618369);
706         final Binary64 ym     = new Binary64(115.055899);
707         final Binary64 sigmaX = new Binary64(15.988242);
708         final Binary64 sigmaY = new Binary64(5756.840725);
709         final Binary64 radius = new Binary64(1.3);
710 
711         // WHEN
712         final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
713 
714         // THEN
715         Assertions.assertEquals(2.0553e-11, result.getValue().getReal(), 1e-15);
716         Assertions.assertEquals(2.0101e-11, result.getLowerLimit().getReal(), 1e-15);
717         Assertions.assertEquals(2.0557e-11, result.getUpperLimit().getReal(), 1e-15);
718     }
719 
720     @Test
721     @DisplayName("CSM test case 3 Field version")
722     void CsmTestCase3Field() {
723         // GIVEN
724         final Binary64 xm     = new Binary64(102.177247);
725         final Binary64 ym     = new Binary64(693.405893);
726         final Binary64 sigmaX = new Binary64(94.230921);
727         final Binary64 sigmaY = new Binary64(643.409272);
728         final Binary64 radius = new Binary64(5.3);
729 
730         // WHEN
731         final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
732 
733         // THEN
734         Assertions.assertEquals(7.2003e-5, result.getValue().getReal(), 1e-9);
735         Assertions.assertEquals(7.194e-5, result.getLowerLimit().getReal(), 1e-8);
736         Assertions.assertEquals(7.200e-5, result.getUpperLimit().getReal(), 1e-8);
737     }
738 
739     @Test
740     @DisplayName("CDM test case 1 Field version")
741     void CdmTestCase1Field() {
742         // GIVEN
743         final Binary64 xm     = new Binary64(-752.672701);
744         final Binary64 ym     = new Binary64(644.939441);
745         final Binary64 sigmaX = new Binary64(445.859950);
746         final Binary64 sigmaY = new Binary64(6095.858688);
747         final Binary64 radius = new Binary64(3.5);
748 
749         // WHEN
750         final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
751 
752         // THEN
753         Assertions.assertEquals(5.3904e-7, result.getValue().getReal(), 1e-11);
754         Assertions.assertEquals(5.3902e-7, result.getLowerLimit().getReal(), 1e-11);
755         Assertions.assertEquals(5.3904e-7, result.getUpperLimit().getReal(), 1e-11);
756     }
757 
758     @Test
759     @DisplayName("CDM test case 2 Field version")
760     void CdmTestCase2Field() {
761         // GIVEN
762         final Binary64 xm     = new Binary64(-692.362272);
763         final Binary64 ym     = new Binary64(4475.456261);
764         final Binary64 sigmaX = new Binary64(193.454603);
765         final Binary64 sigmaY = new Binary64(562.027293);
766         final Binary64 radius = new Binary64(13.2);
767 
768         // WHEN
769         final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
770 
771         // THEN
772         Assertions.assertEquals(2.2796e-20, result.getValue().getReal(), 1e-24);
773         Assertions.assertEquals(2.2517e-20, result.getLowerLimit().getReal(), 1e-24);
774         Assertions.assertEquals(2.2797e-20, result.getUpperLimit().getReal(), 1e-24);
775     }
776 
777     @Test
778     @DisplayName("Alfano test case 3 Field version.")
779     void AlfanoTestCase3Field() {
780         // GIVEN
781         final Binary64 xm     = new Binary64(-3.8872073);
782         final Binary64 ym     = new Binary64(0.1591646);
783         final Binary64 sigmaX = new Binary64(1.4101830);
784         final Binary64 sigmaY = new Binary64(114.2585190);
785         final Binary64 radius = new Binary64(15);
786 
787         // WHEN
788         final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
789 
790         // THEN
791         Assertions.assertEquals(1.0038e-1, result.getValue().getReal(), 1e-5);
792     }
793 
794     /**
795      * Alfano test case 5.
796      * <p>
797      * Initial expected result of 4.4509e-2 has been slightly modified to 4.4515e-2 as this test case is numerically
798      * unstable. This has been validated with Romain SERRA.
799      * <p>
800      * Result obtained here is (slightly) different that what is obtained with the non field version because the non-field
801      * uses the {@link org.hipparchus.util.MathArrays} linearCombination method to be more precise. However, the field
802      * equivalent does not exist.
803      */
804     @Test
805     @DisplayName("Alfano test case 5 Field version")
806     void AlfanoTestCase5Field() {
807         // GIVEN
808         final Binary64 xm     = new Binary64(-1.2217895);
809         final Binary64 ym     = new Binary64(2.1230067);
810         final Binary64 sigmaX = new Binary64(0.0373279);
811         final Binary64 sigmaY = new Binary64(177.8109003);
812         final Binary64 radius = new Binary64(10);
813 
814         // WHEN
815         final FieldProbabilityOfCollision<Binary64> result = method.compute(xm, ym, sigmaX, sigmaY, radius);
816 
817         // THEN
818         Assertions.assertEquals(4.4515e-2, result.getValue().getReal(), 1e-6);
819     }
820 
821     @Test
822     @DisplayName("Chan test case 04 with custom method constructor and relative boundaries close enough to the desired accuracy Field version")
823     void ChanTestCase04WithCustomLaasMethodField() {
824         // GIVEN
825         final ShortTermEncounter2DPOCMethod customMethod = new Laas2015(1e-4, 1000);
826 
827         final Binary64 xm     = new Binary64(10);
828         final Binary64 ym     = new Binary64(0);
829         final Binary64 sigmaX = new Binary64(25);
830         final Binary64 sigmaY = new Binary64(75);
831         final Binary64 radius = new Binary64(5);
832 
833         // WHEN
834         final FieldProbabilityOfCollision<Binary64> result = customMethod.compute(xm, ym, sigmaX, sigmaY, radius);
835 
836         // THEN
837         Assertions.assertEquals(6.1e-3, result.getValue().getReal(), 1e-4);
838     }
839 
840     @Test
841     @DisplayName(
842             "Test probability from a real CDM. It has been mentioned in https://forum.orekit.org/t/collision-package-in-development/2638/16"
843                     + " and define an expected probability of collision of 0.004450713 which is higher than the value found using"
844                     + " the methods provided in this ssa package because the CSpOC use the projection of a square containing the"
845                     + " collision disk instead of the collision disk directly. Hence an overestimation of the probability of"
846                     + " collision. This test serves as a non regression test")
847     void testReturnExpectedProbabilityFromACdmField() {
848 
849         // GIVEN
850         final String     cdmPath = "/ccsds/cdm/ION_SCV8_vs_STARLINK_1233.txt";
851         final DataSource data    = new DataSource(cdmPath, () -> getClass().getResourceAsStream(cdmPath));
852         final Cdm        cdm     = new ParserBuilder().buildCdmParser().parseMessage(data);
853 
854         // Radii taken from comments in the conjunction data message
855         final Binary64 primaryRadius   = new Binary64(5);
856         final Binary64 secondaryRadius = new Binary64(5);
857 
858         // WHEN
859         final FieldProbabilityOfCollision<Binary64> result = method.compute(cdm, primaryRadius, secondaryRadius);
860 
861         // THEN
862         Assertions.assertEquals(0.0034965176443840893, result.getValue().getReal(), 1e-19);
863     }
864 
865     @Test
866     @DisplayName("Test impact on the probability of collision of a slight difference in combined radius in Chan test case 4")
867     void testReturnExpectedValueWhenIntroducingSmallChangeOnCombinedRadius() {
868         // GIVEN
869         final DSFactory factory = new DSFactory(1, 10);
870 
871         final double xmNominal     = 10;
872         final double ymNominal     = 0;
873         final double sigmaXNominal = 25;
874         final double sigmaYNominal = 75;
875         final double radiusNominal = 5;
876         final double dRadius       = 2.5;
877 
878         final DerivativeStructure xm     = factory.constant(xmNominal);
879         final DerivativeStructure ym     = factory.constant(ymNominal);
880         final DerivativeStructure sigmaX = factory.constant(sigmaXNominal);
881         final DerivativeStructure sigmaY = factory.constant(sigmaYNominal);
882         final DerivativeStructure radius = factory.variable(0, radiusNominal);
883 
884         // WHEN
885         final FieldProbabilityOfCollision<DerivativeStructure> resultNominal =
886                 method.compute(xm, ym, sigmaX, sigmaY, radius);
887         final double taylorResult = resultNominal.getValue().taylor(dRadius);
888         final double exactResult =
889                 method.compute(xmNominal, ymNominal, sigmaXNominal, sigmaYNominal, radiusNominal + dRadius).getValue();
890 
891         // THEN
892         Assertions.assertEquals(6.1e-3, resultNominal.getValue().getReal(), 1e-4);
893         Assertions.assertEquals(exactResult, taylorResult, 1e-16);
894     }
895 
896     @Test
897     @DisplayName("Alfano test case 1 (CDM) from NASA CARA")
898     void AlfanoCDMTestCase01() {
899         // Inputs from NASA CARA
900         String cdmName = "AlfanoTestCase01.cdm";
901         final double combinedHbr = 15.;
902 
903         // Excepted outcome
904         final double expectedPc = 0.146749549;
905         final double tolerance = 1e-6;
906 
907         computeAndCheckCollisionProbability(cdmName, combinedHbr, expectedPc, tolerance);
908     }
909 
910     @Test
911     @DisplayName("Alfano test case 2 (CDM) from NASA CARA")
912     void AlfanoCDMTestCase02() {
913         // Inputs from NASA CARA
914         String cdmName = "AlfanoTestCase02.cdm";
915         final double combinedHbr = 4.;
916 
917         // Excepted outcome
918         final double expectedPc = 0.006222267;
919         final double tolerance = 1e-6;
920 
921         computeAndCheckCollisionProbability(cdmName, combinedHbr, expectedPc, tolerance);
922     }
923 
924     @Test
925     @DisplayName("Alfano test case 3 (CDM) from NASA CARA")
926     void AlfanoCDMTestCase03() {
927         // Inputs from NASA CARA
928         String cdmName = "AlfanoTestCase03.cdm";
929         final double combinedHbr = 15.;
930 
931         // Excepted outcome
932         final double expectedPc = 0.100351176;
933         final double tolerance = 1e-6;
934 
935         computeAndCheckCollisionProbability(cdmName, combinedHbr, expectedPc, tolerance);
936     }
937 
938     @Test
939     @DisplayName("Alfano test case 4 (CDM) from NASA CARA")
940     void AlfanoCDMTestCase04() {
941         // Inputs from NASA CARA
942         String cdmName = "AlfanoTestCase04.cdm";
943         final double combinedHbr = 15.;
944 
945         // Excepted outcome
946         final double expectedPc = 0.049323406;
947         final double tolerance = 1e-5;
948 
949         computeAndCheckCollisionProbability(cdmName, combinedHbr, expectedPc, tolerance);
950     }
951 
952     @Test
953     @DisplayName("Alfano test case 5 (CDM) from NASA CARA")
954     void AlfanoCDMTestCase05() {
955         // Inputs from NASA CARA
956         String cdmName = "AlfanoTestCase05.cdm";
957         final double combinedHbr = 10.;
958 
959         // Excepted outcome
960         final double expectedPc = 0.044487386;
961         final double tolerance = 1e-5;
962 
963         computeAndCheckCollisionProbability(cdmName, combinedHbr, expectedPc, tolerance);
964     }
965 
966     @Test
967     @DisplayName("Alfano test case 6 (CDM) from NASA CARA")
968     void AlfanoCDMTestCase06() {
969         // Inputs from NASA CARA
970         String cdmName = "AlfanoTestCase06.cdm";
971         final double combinedHbr = 10.;
972 
973         // Excepted outcome
974         final double expectedPc = 0.004335455;
975         final double tolerance = 1e-8;
976 
977         computeAndCheckCollisionProbability(cdmName, combinedHbr, expectedPc, tolerance);
978     }
979 
980     @Test
981     @DisplayName("Alfano test case 7 (CDM) from NASA CARA")
982     void AlfanoCDMTestCase07() {
983         // Inputs from NASA CARA
984         String cdmName = "AlfanoTestCase07.cdm";
985         final double combinedHbr = 10.;
986 
987         // Excepted outcome
988         final double expectedPc = 0.000158147;
989         final double tolerance = 1e-9;
990 
991         computeAndCheckCollisionProbability(cdmName, combinedHbr, expectedPc, tolerance);
992     }
993 
994     @Test
995     @DisplayName("Alfano test case 8 (CDM) from NASA CARA")
996     void AlfanoCDMTestCase08() {
997         // Inputs from NASA CARA
998         String cdmName = "AlfanoTestCase08.cdm";
999         final double combinedHbr = 4.;
1000 
1001         // Excepted outcome
1002         final double expectedPc = 0.036948008;
1003         final double tolerance = 1e-5;
1004 
1005         computeAndCheckCollisionProbability(cdmName, combinedHbr, expectedPc, tolerance);
1006     }
1007 
1008     @Test
1009     @DisplayName("Alfano test case 9 (CDM) from NASA CARA")
1010     void AlfanoCDMTestCase09() {
1011         // Inputs from NASA CARA
1012         String cdmName = "AlfanoTestCase09.cdm";
1013         final double combinedHbr = 6.;
1014 
1015         // Excepted outcome
1016         final double expectedPc = 0.290146291;
1017         final double tolerance = 2e-5;
1018 
1019         computeAndCheckCollisionProbability(cdmName, combinedHbr, expectedPc, tolerance);
1020     }
1021 
1022     // Alfano test case 10 is a duplicate of test 9 in the NASA CARA unit tests.
1023 
1024     @Test
1025     @DisplayName("Alfano test case 11 (CDM) from NASA CARA")
1026     void AlfanoCDMTestCase11() {
1027         // Inputs from NASA CARA
1028         String cdmName = "AlfanoTestCase11.cdm";
1029         final double combinedHbr = 4.;
1030 
1031         // Excepted outcome
1032         final double expectedPc = 0.002672026;
1033         final double tolerance = 1e-7;
1034 
1035         computeAndCheckCollisionProbability(cdmName, combinedHbr, expectedPc, tolerance);
1036     }
1037 
1038     private void computeAndCheckCollisionProbability(
1039             final String cdmName,
1040             final double combinedHbr,
1041             final double expected,
1042             final double tolerance) {
1043 
1044         // Given
1045         final String     cdmPath = "/ccsds/cdm/" + cdmName;
1046         final DataSource data    = new DataSource(cdmPath, () -> getClass().getResourceAsStream(cdmPath));
1047         final Cdm        cdm     = new ParserBuilder().buildCdmParser().parseMessage(data);
1048 
1049         // When
1050         final ProbabilityOfCollision result = method.compute(cdm, combinedHbr);
1051 
1052         // Then
1053         Assertions.assertEquals(expected, result.getValue(), tolerance);
1054     }
1055 }