From 0b2aea7144b45d4fc3c7a26c8b107b93ac670fbd Mon Sep 17 00:00:00 2001
From: Luc Maisonobe <luc@orekit.org>
Date: Wed, 10 Sep 2014 17:09:34 +0200
Subject: [PATCH] Fixed errors in basic scan algorithm.

Near tile boundaries, the algorithm sometime failed to find intersection
points.
---
 .../intersection/BasicScanAlgorithm.java      | 20 +++---
 .../org/orekit/rugged/api/RuggedTest.java     | 68 ++++++++++++++++++-
 2 files changed, 76 insertions(+), 12 deletions(-)

diff --git a/core/src/main/java/org/orekit/rugged/intersection/BasicScanAlgorithm.java b/core/src/main/java/org/orekit/rugged/intersection/BasicScanAlgorithm.java
index 99e2cc4c..fca8fa28 100644
--- a/core/src/main/java/org/orekit/rugged/intersection/BasicScanAlgorithm.java
+++ b/core/src/main/java/org/orekit/rugged/intersection/BasicScanAlgorithm.java
@@ -87,7 +87,7 @@ public class BasicScanAlgorithm implements IntersectionAlgorithm {
                 exitPoint = ellipsoid.transform(ellipsoid.pointAtAltitude(position, los, Double.isInfinite(hMin) ? 0.0 : hMin),
                                                 ellipsoid.getBodyFrame(), null);
                 final SimpleTile exitTile = cache.getTile(exitPoint.getLatitude(), exitPoint.getLongitude());
-                addIfNotPresent(scannedTiles, entryTile);
+                addIfNotPresent(scannedTiles, exitTile);
 
                 minLatitude  = FastMath.min(entryPoint.getLatitude(),  exitPoint.getLatitude());
                 maxLatitude  = FastMath.max(entryPoint.getLatitude(),  exitPoint.getLatitude());
@@ -116,7 +116,7 @@ public class BasicScanAlgorithm implements IntersectionAlgorithm {
             for (final SimpleTile tile : scannedTiles) {
                 for (int i = latitudeIndex(tile, minLatitude); i <= latitudeIndex(tile, maxLatitude); ++i) {
                     for (int j = longitudeIndex(tile, minLongitude); j <= longitudeIndex(tile, maxLongitude); ++j) {
-                        GeodeticPoint gp = tile.pixelIntersection(entryPoint, ellipsoid.convertLos(entryPoint, los), i, j);
+                        final GeodeticPoint gp = tile.pixelIntersection(entryPoint, ellipsoid.convertLos(entryPoint, los), i, j);
                         if (gp != null) {
 
                             // improve the point, by projecting it back on the 3D line, fixing the small body curvature at pixel level
@@ -124,13 +124,15 @@ public class BasicScanAlgorithm implements IntersectionAlgorithm {
                             final double        s         = Vector3D.dotProduct(delta, los) / los.getNormSq();
                             final GeodeticPoint projected = ellipsoid.transform(new Vector3D(1, position, s, los),
                                                                                 ellipsoid.getBodyFrame(), null);
-                            gp = tile.pixelIntersection(projected, ellipsoid.convertLos(projected, los), i, j);
-
-                            final Vector3D point = ellipsoid.transform(gp);
-                            final double dot = Vector3D.dotProduct(point.subtract(position), los);
-                            if (dot < intersectionDot) {
-                                intersectionGP  = gp;
-                                intersectionDot = dot;
+                            final GeodeticPoint gpImproved = tile.pixelIntersection(projected, ellipsoid.convertLos(projected, los), i, j);
+
+                            if (gpImproved != null) {
+                                final Vector3D point = ellipsoid.transform(gpImproved);
+                                final double dot = Vector3D.dotProduct(point.subtract(position), los);
+                                if (dot < intersectionDot) {
+                                    intersectionGP  = gpImproved;
+                                    intersectionDot = dot;
+                                }
                             }
 
                         }
diff --git a/core/src/test/java/org/orekit/rugged/api/RuggedTest.java b/core/src/test/java/org/orekit/rugged/api/RuggedTest.java
index 63bc4b18..8697873b 100644
--- a/core/src/test/java/org/orekit/rugged/api/RuggedTest.java
+++ b/core/src/test/java/org/orekit/rugged/api/RuggedTest.java
@@ -148,7 +148,7 @@ public class RuggedTest {
                                                  FastMath.toRadians(1.0), 1201);
 
         Rugged rugged = new Rugged(updater, 8, AlgorithmId.DUVENHAGE,
-                                   EllipsoidId.WGS84, InertialFrameId.EME2000, BodyRotatingFrameId.ITRF,
+                                   EllipsoidId.GRS80, InertialFrameId.EME2000, BodyRotatingFrameId.ITRF,
                                    pv.get(0).getDate(), pv.get(pv.size() - 1).getDate(), 5.0,
                                    pv, 8, CartesianDerivativesFilter.USE_PV, q, 8, AngularDerivativesFilter.USE_R, 0.001);
 
@@ -178,7 +178,7 @@ public class RuggedTest {
                                                  FastMath.toRadians(1.0), 1201);
 
         Rugged rugged = new Rugged(updater, 8, AlgorithmId.DUVENHAGE,
-                                   EllipsoidId.WGS84, InertialFrameId.EME2000, BodyRotatingFrameId.ITRF,
+                                   EllipsoidId.IERS96, InertialFrameId.EME2000, BodyRotatingFrameId.ITRF,
                                    orbit.getDate().shiftedBy(-10.0), orbit.getDate().shiftedBy(+10.0), 5.0,
                                    1.0, 8, CartesianDerivativesFilter.USE_PV, AngularDerivativesFilter.USE_R, 0.001, propagator);
 
@@ -392,7 +392,7 @@ public class RuggedTest {
         AbsoluteDate maxDate = lineSensor.getDate(lastLine);
 
         Rugged rugged = new Rugged(null, -1, AlgorithmId.IGNORE_DEM_USE_ELLIPSOID,
-                                   EllipsoidId.WGS84, InertialFrameId.EME2000, BodyRotatingFrameId.ITRF,
+                                   EllipsoidId.IERS2003, InertialFrameId.EME2000, BodyRotatingFrameId.ITRF,
                                    minDate, maxDate, 5.0,
                                    orbitToPV(orbit, earth, minDate.shiftedBy(-1.0), maxDate.shiftedBy(+1.0), 0.25), 8,
                                    CartesianDerivativesFilter.USE_PV,
@@ -538,6 +538,68 @@ public class RuggedTest {
 
     }
 
+    @Test
+    public void testBasicScan()
+        throws RuggedException, OrekitException, URISyntaxException {
+
+        int dimension = 200;
+
+        String path = getClass().getClassLoader().getResource("orekit-data").toURI().getPath();
+        DataProvidersManager.getInstance().addProvider(new DirectoryCrawler(new File(path)));
+        final BodyShape  earth = createEarth();
+        final Orbit      orbit = createOrbit(Constants.EIGEN5C_EARTH_MU);
+
+        AbsoluteDate crossing = new AbsoluteDate("2012-01-01T12:30:00.000", TimeScalesFactory.getUTC());
+
+        // one line sensor
+        // position: 1.5m in front (+X) and 20 cm above (-Z) of the S/C center of mass
+        // los: swath in the (YZ) plane, looking at 50° roll, ±1° aperture
+        Vector3D position = new Vector3D(1.5, 0, -0.2);
+        List<Vector3D> los = createLOS(new Rotation(Vector3D.PLUS_I, FastMath.toRadians(50.0)).applyTo(Vector3D.PLUS_K),
+                                       Vector3D.PLUS_I,
+                                       FastMath.toRadians(1.0), dimension);
+
+        // linear datation model: at reference time we get line 100, and the rate is one line every 1.5ms
+        LineDatation lineDatation = new LinearLineDatation(crossing, dimension / 2, 1.0 / 1.5e-3);
+        int firstLine = 0;
+        int lastLine  = dimension;
+        LineSensor lineSensor = new LineSensor("line", lineDatation, position, los);
+        AbsoluteDate minDate = lineSensor.getDate(firstLine);
+        AbsoluteDate maxDate = lineSensor.getDate(lastLine);
+
+        TileUpdater updater =
+                new RandomLandscapeUpdater(0.0, 9000.0, 0.5, 0xe12ef744f224cf43l,
+                                           FastMath.toRadians(1.0), 257);
+
+        Rugged ruggedDuvenhage = new Rugged(updater, 8, AlgorithmId.DUVENHAGE,
+                                            EllipsoidId.WGS84, InertialFrameId.EME2000, BodyRotatingFrameId.ITRF,
+                                            minDate, maxDate, 5.0,
+                                            orbitToPV(orbit, earth, minDate.shiftedBy(-1.0), maxDate.shiftedBy(+1.0), 0.25), 8,
+                                            CartesianDerivativesFilter.USE_PV,
+                                            orbitToQ(orbit, earth, minDate.shiftedBy(-1.0), maxDate.shiftedBy(+1.0), 0.25), 2,
+                                            AngularDerivativesFilter.USE_R, 0.001);
+        ruggedDuvenhage.addLineSensor(lineSensor);
+        GeodeticPoint[] gpDuvenhage = ruggedDuvenhage.directLocalization("line", 100);
+
+        Rugged ruggedBasicScan = new Rugged(updater, 8, AlgorithmId.BASIC_SLOW_EXHAUSTIVE_SCAN_FOR_TESTS_ONLY,
+                                            EllipsoidId.WGS84, InertialFrameId.EME2000, BodyRotatingFrameId.ITRF,
+                                            minDate, maxDate, 5.0,
+                                            orbitToPV(orbit, earth, minDate.shiftedBy(-1.0), maxDate.shiftedBy(+1.0), 0.25), 8,
+                                            CartesianDerivativesFilter.USE_PV,
+                                            orbitToQ(orbit, earth, minDate.shiftedBy(-1.0), maxDate.shiftedBy(+1.0), 0.25), 2,
+                                            AngularDerivativesFilter.USE_R, 0.001);
+        ruggedBasicScan.addLineSensor(lineSensor);
+        GeodeticPoint[] gpBasicScan = ruggedBasicScan.directLocalization("line", 100);
+
+
+        for (int i = 0; i < gpDuvenhage.length; ++i) {
+            Vector3D pDuvenhage = earth.transform(gpDuvenhage[i]);
+            Vector3D pBasicScan = earth.transform(gpBasicScan[i]);
+            Assert.assertEquals(0.0, Vector3D.distance(pDuvenhage, pBasicScan), 4.0e-5);
+        }
+
+    }
+
     @Test
     public void testInterpolatorDump()
         throws RuggedException, OrekitException, URISyntaxException {
-- 
GitLab