1   /* Copyright 2013-2020 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.rugged.utils;
18  
19  import java.io.File;
20  import java.net.URISyntaxException;
21  
22  import org.hipparchus.geometry.euclidean.threed.Line;
23  import org.hipparchus.geometry.euclidean.threed.Vector3D;
24  import org.hipparchus.util.FastMath;
25  import org.junit.After;
26  import org.junit.Assert;
27  import org.junit.Before;
28  import org.junit.Test;
29  import org.orekit.bodies.GeodeticPoint;
30  import org.orekit.data.DataContext;
31  import org.orekit.data.DirectoryCrawler;
32  import org.orekit.errors.OrekitException;
33  import org.orekit.frames.Frame;
34  import org.orekit.frames.FramesFactory;
35  import org.orekit.rugged.errors.RuggedException;
36  import org.orekit.rugged.errors.RuggedMessages;
37  import org.orekit.utils.Constants;
38  import org.orekit.utils.IERSConventions;
39  
40  public class ExtendedEllipsoidTest {
41  
42      @Test
43      public void testPointAtLongitude() {
44  
45          Vector3D p = new Vector3D(3220103.0, 69623.0, -6449822.0);
46          Vector3D d = new Vector3D(1.0, 2.0, 3.0);
47  
48          for (double longitude = -1.0; longitude < 1.0; longitude += 0.01) {
49              GeodeticPoint gp = ellipsoid.transform(ellipsoid.pointAtLongitude(p, d, longitude),
50                                                     ellipsoid.getBodyFrame(), null);
51              Assert.assertEquals(longitude, gp.getLongitude(), 1.0e-15);
52          }
53  
54      }
55  
56      @Test
57      public void testPointAtLongitudeError() {
58  
59          Vector3D p = new Vector3D(3220103.0, 69623.0, -6449822.0);
60          double longitude = 1.25;
61          Vector3D parallelToLongitudePlane = new Vector3D(FastMath.cos(longitude),
62                                                           FastMath.sin(longitude),
63                                                           -2.4);
64          try {
65              ellipsoid.pointAtLongitude(p, parallelToLongitudePlane, longitude);
66              Assert.fail("an error should have been triggered");
67          } catch (RuggedException re) {
68              Assert.assertEquals(RuggedMessages.LINE_OF_SIGHT_NEVER_CROSSES_LONGITUDE, re.getSpecifier());
69              Assert.assertEquals(FastMath.toDegrees(longitude), re.getParts()[0]);
70          }
71  
72      }
73  
74      @Test
75      public void testPointAtLatitude() {
76  
77          Vector3D p = new Vector3D(3220103.0, 69623.0, -6449822.0);
78          Vector3D d = new Vector3D(1.0, 2.0, 3.0);
79  
80          double epsilon = 5.0e-15;
81          for (double latitude = -d.getDelta() + 1.0e-5; latitude < d.getDelta(); latitude += 0.1) {
82              GeodeticPoint gp = ellipsoid.transform(ellipsoid.pointAtLatitude(p, d, latitude, p),
83                                                     ellipsoid.getBodyFrame(), null);
84              Assert.assertEquals(latitude, gp.getLatitude(), epsilon);
85          }
86  
87      }
88  
89      @Test
90      public void testPointAtLatitudeTwoPointsSameSide() {
91  
92          // the line of sight is almost parallel an iso-latitude cone generatrix
93          // the spacecraft is at latitude lTarget - 0.951", and altitude 794.6km
94          // so at a latitude slightly less than the target
95          // the line of sight crosses the latitude cone first about 70km along line of sight
96          // (so still at a very high altitude) and a second time about 798km along line of sight,
97          // only a few hundreds meters above allipsoid
98          // Note that this happens despite the line of sight is not along nadir, the longitudes
99          // of the spacecraft and crossing points span in a 0.88° wide longitude range
100         Vector3D position     = new Vector3D(-748528.2769999998, -5451658.432000002, 4587158.354);
101         Vector3D los          = new Vector3D(0.010713435156834539, 0.7688536080293823, -0.6393350856376809);
102         double h              = ellipsoid.transform(position, ellipsoid.getBodyFrame(), null).getAltitude();
103         double lTarget        = 0.6978408125890662;
104 
105         // spacecraft is in LEO
106         Assert.assertEquals(h, 794652.782, 0.001);
107 
108         Vector3D pHigh        = ellipsoid.pointAtLatitude(position, los, lTarget, position);
109         GeodeticPoint gpHigh  = ellipsoid.transform(pHigh, ellipsoid.getBodyFrame(), null);
110         Assert.assertEquals(lTarget, gpHigh.getLatitude(), 1.0e-12);
111         // first crossing point is high, but below spacecraft and along positive line of sight
112         Assert.assertEquals(724335.409, gpHigh.getAltitude(), 0.001);
113         Assert.assertTrue(Vector3D.dotProduct(pHigh.subtract(position), los) > 0);
114 
115         Vector3D pLow         = ellipsoid.pointAtLatitude(position, los, lTarget,
116                                                           new Vector3D(1, position, 900000, los));
117         GeodeticPoint gpLow   = ellipsoid.transform(pLow, ellipsoid.getBodyFrame(), null);
118         Assert.assertEquals(lTarget, gpLow.getLatitude(), 1.0e-12);
119         // second crossing point is almost on ground, also along positive line of sight
120         Assert.assertEquals(492.804, gpLow.getAltitude(), 0.001);
121         Assert.assertTrue(Vector3D.dotProduct(pLow.subtract(position), los) > 0);
122 
123     }
124 
125     @Test
126     public void testPointAtLatitudeTwoPointsOppositeSides() {
127 
128         Vector3D p = new Vector3D(3220103.0, 69623.0, -6449822.0);
129         Vector3D d = new Vector3D(1.0, 2.0, 0.1);
130         double latitude = -0.5;
131 
132         Vector3D pPlus = ellipsoid.pointAtLatitude(p, d, latitude, new Vector3D(1, p, +2.0e7, d));
133         GeodeticPoint gpPlus = ellipsoid.transform(pPlus, ellipsoid.getBodyFrame(), null);
134         Assert.assertEquals(latitude, gpPlus.getLatitude(), 4.0e-16);
135         Assert.assertEquals(20646364.047, Vector3D.dotProduct(d, pPlus.subtract(p)), 0.001);
136 
137         Vector3D pMinus = ellipsoid.pointAtLatitude(p, d, latitude, new Vector3D(1, p, -3.0e7, d));
138         GeodeticPoint gpMinus = ellipsoid.transform(pMinus, ellipsoid.getBodyFrame(), null);
139         Assert.assertEquals(latitude, gpMinus.getLatitude(), 3.0e-16);
140         Assert.assertEquals(-31797895.234, Vector3D.dotProduct(d, pMinus.subtract(p)), 0.001);
141 
142     }
143 
144     @Test
145     public void testPointAtLatitudeAlmostEquator() {
146         Vector3D      p              = new Vector3D(5767483.098580201, 4259689.325372237, -41553.67750784925);
147         Vector3D      d              = new Vector3D(-0.7403523952347795, -0.6701811835520302, 0.05230212180799747);
148         double        latitude       = -3.469446951953614E-18;
149         Vector3D      closeReference = new Vector3D(5177991.74844521, 3726070.452427455, 90.88067547897226);
150         Vector3D      intersection   = ellipsoid.pointAtLatitude(p, d, latitude, closeReference);
151         GeodeticPoint gp             = ellipsoid.transform(intersection, ellipsoid.getBodyFrame(), null);
152         Assert.assertEquals(latitude, gp.getLatitude(), 1.0e-10);
153         Assert.assertEquals(2866.297, gp.getAltitude(), 1.0e-3);
154     }
155 
156     @Test
157     public void testPointAtLatitudeErrorQuadraticEquation() {
158 
159         Vector3D p = new Vector3D(3220103.0, 69623.0, -6449822.0);
160         Vector3D d = new Vector3D(1.0, 2.0, 3.0);
161         double latitude = -1.4;
162 
163         try {
164             ellipsoid.pointAtLatitude(p, d, latitude, p);
165             Assert.fail("an error should have been triggered");
166         } catch (RuggedException re) {
167             Assert.assertEquals(RuggedMessages.LINE_OF_SIGHT_NEVER_CROSSES_LATITUDE, re.getSpecifier());
168             Assert.assertEquals(FastMath.toDegrees(latitude), re.getParts()[0]);
169         }
170 
171     }
172 
173     @Test
174     public void testPointAtLatitudeErrorNappe() {
175 
176         Vector3D p = new Vector3D(3220103.0, 69623.0, -6449822.0);
177         Vector3D d = new Vector3D(1.0, 2.0, 0.1);
178         double latitude = 0.5;
179 
180         try {
181             ellipsoid.pointAtLatitude(p, d, latitude, p);
182             Assert.fail("an error should have been triggered");
183         } catch (RuggedException re) {
184             Assert.assertEquals(RuggedMessages.LINE_OF_SIGHT_NEVER_CROSSES_LATITUDE, re.getSpecifier());
185             Assert.assertEquals(FastMath.toDegrees(latitude), re.getParts()[0]);
186         }
187 
188     }
189 
190     @Test
191     public void testPointAtAltitude() {
192 
193         Vector3D p = new Vector3D(3220103.0, 69623.0, -6449822.0);
194         Vector3D d = new Vector3D(1.0, 2.0, 3.0);
195         for (double altitude = -500000; altitude < 800000.0; altitude += 100) {
196             GeodeticPoint gp = ellipsoid.transform(ellipsoid.pointAtAltitude(p, d, altitude),
197                                                    ellipsoid.getBodyFrame(), null);
198             Assert.assertEquals(altitude, gp.getAltitude(), 1.0e-3);
199         }
200 
201     }
202 
203     @Test
204     public void testPointAtAltitudeStartInside() {
205 
206         Vector3D p = new Vector3D(322010.30, 6962.30, -644982.20);
207         Vector3D d = new Vector3D(-1.0, -2.0, -3.0);
208         for (double altitude = -500000; altitude < 800000.0; altitude += 100) {
209             GeodeticPoint gp = ellipsoid.transform(ellipsoid.pointAtAltitude(p, d, altitude),
210                                                    ellipsoid.getBodyFrame(), null);
211             Assert.assertEquals(altitude, gp.getAltitude(), 1.0e-3);
212         }
213 
214     }
215 
216     @Test
217     public void testPointAtAltitudeError() {
218 
219         Vector3D p = new Vector3D(3220103.0, 69623.0, -6449822.0);
220         Vector3D d = new Vector3D(1.0, 2.0, 3.0);
221         double altitude = -580000.0;
222         try {
223             ellipsoid.pointAtAltitude(p, d, altitude);
224             Assert.fail("an error should have been triggered");
225         } catch (RuggedException re) {
226             Assert.assertEquals(RuggedMessages.LINE_OF_SIGHT_NEVER_CROSSES_ALTITUDE, re.getSpecifier());
227             Assert.assertEquals(altitude, re.getParts()[0]);
228         }
229 
230     }
231 
232     @Test
233     public void testConvertLOS() {
234 
235         GeodeticPoint gp   = new GeodeticPoint(-0.2, 1.8, 2400.0);
236         Vector3D p         = ellipsoid.transform(gp);
237         Vector3D los       = new Vector3D(-1, -2, -3);
238         Vector3D converted = ellipsoid.convertLos(gp, los);
239         Line line = new Line(p, new Vector3D(1.0, p, 1000, los), 1.0e-10);
240 
241         for (double delta = 0.1; delta < 100.0; delta += 0.1) {
242             GeodeticPoint shifted = new GeodeticPoint(gp.getLatitude()  + delta * converted.getY(),
243                                                       gp.getLongitude() + delta * converted.getX(),
244                                                       gp.getAltitude()  + delta * converted.getZ());
245             Vector3D converted2 = ellipsoid.convertLos(p, ellipsoid.transform(shifted));
246             Assert.assertEquals(0.0, Vector3D.distance(converted, converted2), 3.0e-5 * converted.getNorm());
247             Assert.assertEquals(0.0, line.distance(ellipsoid.transform(shifted)), 8.0e-4);
248         }
249 
250     }
251 
252     @Test
253     public void testPointAtLatitudeError() {
254 
255         Vector3D p = new Vector3D(-3052690.88784496, 6481300.309857268, 25258.7478104745);
256         Vector3D d = new Vector3D(0.6, -0.8, 0.0);
257         double latitude = 0.1;
258         Vector3D c = new Vector3D(-2809972.5765414005, 5727461.020250551, 26.163518446261833);
259 
260         try {
261             ellipsoid.pointAtLatitude(p, d, latitude, c);
262             Assert.fail("an error should have been triggered");
263         } catch (RuggedException re) {
264             Assert.assertEquals(RuggedMessages.LINE_OF_SIGHT_NEVER_CROSSES_LATITUDE, re.getSpecifier());
265             Assert.assertEquals(FastMath.toDegrees(latitude), re.getParts()[0]);
266         }
267 
268     }
269 
270     @Test
271     public void testPointAtLatitudeIssue1() {
272 
273         Vector3D position = new Vector3D(-1988136.619268088, -2905373.394638188, 6231185.484365295);
274         Vector3D los = new Vector3D(0.3489121277213534, 0.3447806500507106, -0.8714279261531437);
275         Vector3D close = new Vector3D(-1709383.0948608494, -2630206.8820586684, 5535282.169189105);
276         double latitude =  1.0581058590215624;
277 
278         Vector3D s = ellipsoid.pointAtLatitude(position, los, latitude, close);
279         GeodeticPoint gp = ellipsoid.transform(s, ellipsoid.getBodyFrame(), null);
280         Assert.assertEquals(latitude, gp.getLatitude(), 1.0e-15);
281 
282     }
283 
284     @Before
285     public void setUp() {
286         try {
287 
288             String path = getClass().getClassLoader().getResource("orekit-data").toURI().getPath();
289             DataContext.getDefault().getDataProvidersManager().addProvider(new DirectoryCrawler(new File(path)));
290 
291             Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
292             ellipsoid = new ExtendedEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
293                                               Constants.WGS84_EARTH_FLATTENING,
294                                               itrf);
295         } catch (OrekitException oe) {
296             Assert.fail(oe.getLocalizedMessage());
297         } catch (URISyntaxException use) {
298             Assert.fail(use.getLocalizedMessage());
299         }
300     }
301 
302     @After
303     public void tearDown() {
304         ellipsoid = null;
305     }
306 
307     private ExtendedEllipsoid ellipsoid;
308 
309 }