1   /* Copyright 2002-2025 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.models.earth.tessellation;
18  
19  import java.io.IOException;
20  import java.util.List;
21  
22  import org.hipparchus.geometry.euclidean.threed.Vector3D;
23  import org.hipparchus.geometry.partitioning.Region.Location;
24  import org.hipparchus.geometry.spherical.twod.S2Point;
25  import org.hipparchus.geometry.spherical.twod.SphericalPolygonsSet;
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.frames.FramesFactory;
35  import org.orekit.utils.Constants;
36  import org.orekit.utils.IERSConventions;
37  
38  public class DivertedSingularityAimingTest {
39  
40      @Test
41      public void testSingularityOutside() {
42          Assertions.assertEquals(1, aiming.getSingularPoints().size());
43          final GeodeticPoint singularity = aiming.getSingularPoints().get(0);
44          Assertions.assertEquals(Location.OUTSIDE, aoi.checkPoint(toS2Point(singularity)));
45      }
46  
47      @Test
48      public void testAroundSingularity() throws IOException {
49  
50          GeodeticPoint singularityGP = aiming.getSingularPoints().get(0);
51          Vector3D singularity = earth.transform(singularityGP);
52          Assertions.assertEquals(GeodeticPoint.SOUTH_POLE.getLatitude(), singularityGP.getLatitude(), 1.0e-10);
53  
54          // in a small disk (less than 1cm radius) around singularity, aiming direction changes a lot
55          double lat = GeodeticPoint.SOUTH_POLE.getLatitude() + 1.0e-9;
56          GeodeticPoint gp000  = new GeodeticPoint(lat, 0.0, 0.0);
57          Vector3D      p000   = earth.transform(gp000);
58          Vector3D      dir000 = aiming.alongTileDirection(p000, gp000);
59          GeodeticPoint gp090  = new GeodeticPoint(lat, 0.5 * FastMath.PI, 0.0);
60          Vector3D      p090   = earth.transform(gp090);
61          Vector3D      dir090 = aiming.alongTileDirection(p090, gp090);
62          Assertions.assertEquals(0.0064, Vector3D.distance(singularity, p000), 1.0e-4);
63          Assertions.assertEquals(0.0064, Vector3D.distance(singularity, p090), 1.0e-4);
64          Assertions.assertEquals(FastMath.PI, Vector3D.angle(dir000, dir090), 5.0e-7);
65  
66      }
67  
68      @Test
69      public void testOppositeSingularity() throws IOException {
70  
71          GeodeticPoint singularityGP = aiming.getSingularPoints().get(0);
72          Vector3D singularity = earth.transform(singularityGP);
73          Vector3D opposite    = singularity.negate();
74          GeodeticPoint oppositeGP = earth.transform(opposite, earth.getBodyFrame(), null);
75          Assertions.assertEquals(GeodeticPoint.NORTH_POLE.getLatitude(), oppositeGP.getLatitude(), 1.0e-10);
76  
77          // around opposite of singularity, aiming direction is almost constant
78          // (as we use dipole field to model aiming direction, there is only one singularity)
79          Vector3D refDir = aiming.alongTileDirection(opposite, oppositeGP);
80          double lat = GeodeticPoint.NORTH_POLE.getLatitude() - 1.0e-9;
81          for (double lon = 0; lon < 2 * FastMath.PI; lon += 0.001) {
82              GeodeticPoint gp = new GeodeticPoint(lat, lon, 0.0);
83              Vector3D      p  = earth.transform(gp);
84              Vector3D      dir = aiming.alongTileDirection(p, gp);
85              Assertions.assertEquals(0.0064, Vector3D.distance(opposite, p), 1.0e-4);
86              Assertions.assertEquals(0.0, Vector3D.angle(refDir, dir), 1.1e-9);
87          }
88  
89      }
90  
91      /** Test issue 969 on computation of singularity point.
92       * Issue was due to an Hipparchus bug (see https://github.com/Hipparchus-Math/hipparchus/issues/208).
93       * It was fixed by upgrading to Hipparchus 2.3. 
94       */
95      @Test
96      public void testIssue969() throws IOException {
97  
98          // Given
99          // -----
100         
101         // Zone on Earth
102         final GeodeticPoint northWest = new GeodeticPoint(FastMath.toRadians(30.), FastMath.toRadians(-30.), 10.);
103         final GeodeticPoint southWest = new GeodeticPoint(FastMath.toRadians(-10.), FastMath.toRadians(-30.), 3000.);
104         final GeodeticPoint southEast = new GeodeticPoint(FastMath.toRadians(-10.), FastMath.toRadians(20.), -2000.);
105         final GeodeticPoint northEast = new GeodeticPoint(FastMath.toRadians(30.), FastMath.toRadians(20.), -30.);
106 
107         // When
108         // ----
109         
110         // Counter clockwise zone definition
111         final SphericalPolygonsSet targetZone = EllipsoidTessellator.buildSimpleZone(1.0e-10, northWest, southWest,
112                                                                                      southEast, northEast);
113         // Build DivertedSingularityAiming
114         final TileAiming tileAiming = new DivertedSingularityAiming(targetZone);
115 
116         // Then
117         // ----
118         
119         // Check center and singularity
120         final S2Point centerZone = targetZone.getEnclosingCap().getCenter();
121         final GeodeticPoint centerGP = new GeodeticPoint(0.5 * FastMath.PI - centerZone.getPhi(), centerZone.getTheta(), 0.0);       
122 
123         // Get singularity point (there should be just one)
124         List<GeodeticPoint> singularGPs = tileAiming.getSingularPoints();
125         final GeodeticPoint singularGP = singularGPs.get(0);
126 
127         // Singular list size
128         Assertions.assertEquals(1, singularGPs.size());
129 
130         // Check center
131         Assertions.assertEquals(9.0794674733, FastMath.toDegrees(centerGP.getLatitude()), 1.0e-10);
132         Assertions.assertEquals(-5., FastMath.toDegrees(centerGP.getLongitude()), 1.0e-14);
133         Assertions.assertEquals(0., FastMath.toDegrees(centerGP.getAltitude()), 0.);
134 
135         // Check singularity (should be at the antipodes of center)
136         Assertions.assertEquals(-9.0794674733, FastMath.toDegrees(singularGP.getLatitude()), 1.0e-10);
137         Assertions.assertEquals(175., FastMath.toDegrees(singularGP.getLongitude()), 1.0e-13);
138         Assertions.assertEquals(0., FastMath.toDegrees(singularGP.getAltitude()), 0.);
139 
140     }
141 
142     private S2Point toS2Point(final GeodeticPoint point) {
143         return new S2Point(point.getLongitude(), 0.5 * FastMath.PI - point.getLatitude());
144     }
145 
146     @BeforeEach
147     public void setUp() {
148         Utils.setDataRoot("regular-data");
149         earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
150                                      Constants.WGS84_EARTH_FLATTENING,
151                                      FramesFactory.getITRF(IERSConventions.IERS_2010, true));
152         aoi = new SphericalPolygonsSet(1.e-9, new S2Point[] {
153             new S2Point(FastMath.toRadians(-120.0), FastMath.toRadians(5.0)),
154             new S2Point(FastMath.toRadians(   0.0), FastMath.toRadians(5.0)),
155             new S2Point(FastMath.toRadians( 120.0), FastMath.toRadians(5.0))
156         });
157         aiming = new DivertedSingularityAiming(aoi);
158     }
159 
160     @AfterEach
161     public void tearDown() {
162         earth  = null;
163         aoi    = null;
164         aiming = null;
165     }
166 
167     private OneAxisEllipsoid earth;
168     private SphericalPolygonsSet aoi;
169     private TileAiming aiming;
170 
171 }