Added ExtendedEllipsoid to chop off line-of-sight according to DEM.

rugged-api/src/main/java/org/orekit/rugged/api/RuggedMessages.java

@@ -54,7 +54,10 @@ public enum RuggedMessages implements Localizable {
     OUT_OF_TILE_INDICES("no data at indices [{0}, {1}], tile only covers from [0, 0] to [{2}, {3}] (inclusive)"),
     UNINITIALIZED_CONTEXT("general context has not been initialized"),
     EMPTY_TILE("tile is empty: {0} ⨉ {1}"),
-    UNKNOWN_SENSOR("unknown sensor {0}");
+    UNKNOWN_SENSOR("unknown sensor {0}"),
+    LINE_OF_SIGHT_NEVER_CROSSES_LATITUDE("line-of-sight never crosses latitude {0}"),
+    LINE_OF_SIGHT_NEVER_CROSSES_LONGITUDE("line-of-sight never crosses longitude {0}"),
+    LINE_OF_SIGHT_NEVER_CROSSES_ALTITUDE("line-of-sight never crosses altitude {0}");
 
     // CHECKSTYLE: resume JavadocVariable check
 

rugged-api/src/main/resources/assets/org/orekit/rugged/RuggedMessages_en.utf8

@@ -12,3 +12,12 @@ EMPTY_TILE = tile is empty: {0} ⨉ {1}
 
 # unknown sensor {0}
 UNKNOWN_SENSOR = unknown sensor {0}
+
+# line-of-sight never crosses latitude {0}
+LINE_OF_SIGHT_NEVER_CROSSES_LATITUDE = line-of-sight never crosses latitude {0}
+
+# line-of-sight never crosses longitude {0}
+LINE_OF_SIGHT_NEVER_CROSSES_LONGITUDE = line-of-sight never crosses longitude {0}
+
+# line never crosses altitude {0}
+LINE_OF_SIGHT_NEVER_CROSSES_ALTITUDE = line-of-sight never crosses altitude {0}

rugged-api/src/main/resources/assets/org/orekit/rugged/RuggedMessages_fr.utf8

@@ -12,3 +12,12 @@ EMPTY_TILE = la tuile est vide : {0} ⨉ {1}
 
 # unknown sensor {0}
 UNKNOWN_SENSOR = capteur {0} inconnu
+
+# line-of-sight never crosses latitude {0}
+LINE_OF_SIGHT_NEVER_CROSSES_LATITUDE = la ligne de visée ne franchit jamais la latitude {0}
+
+# line-of-sight never crosses longitude {0}
+LINE_OF_SIGHT_NEVER_CROSSES_LONGITUDE = la ligne de visée ne franchit jamais la longitude {0}
+
+# line never crosses altitude {0}
+LINE_OF_SIGHT_NEVER_CROSSES_ALTITUDE = la ligne de visée ne franchit jamais l''altitude {0}

rugged-api/src/test/java/org/orekit/rugged/api/RuggedMessagesTest.java

@@ -29,7 +29,7 @@ public class RuggedMessagesTest {
 
     @Test
     public void testMessageNumber() {
-        Assert.assertEquals(5, RuggedMessages.values().length);
+        Assert.assertEquals(8, RuggedMessages.values().length);
     }
 
     @Test

rugged-core/src/main/java/org/orekit/rugged/core/ExtendedEllipsoid.java

+/* Copyright 2013-2014 CS Systèmes d'Information
+ * Licensed to CS Systèmes d'Information (CS) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * CS licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.orekit.rugged.core;
+
+import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
+import org.apache.commons.math3.util.FastMath;
+import org.orekit.bodies.GeodeticPoint;
+import org.orekit.bodies.OneAxisEllipsoid;
+import org.orekit.errors.OrekitException;
+import org.orekit.frames.Frame;
+import org.orekit.rugged.api.RuggedException;
+import org.orekit.rugged.api.RuggedMessages;
+
+/** Transform provider from Spacecraft frame to observed body frame.
+ * @author Luc Maisonobe
+ */
+public class ExtendedEllipsoid extends OneAxisEllipsoid {
+
+    /** Serializable UID. */
+    private static final long serialVersionUID = 20140312L;
+
+    /** Convergence threshold for {@link #pointAtAltitude(Vector3D, Vector3D, double)}. */
+    private static final double ALTITUDE_CONVERGENCE = 1.0e-3;
+
+    /** Simple constructor.
+     * @param ae equatorial radius
+     * @param f the flattening (f = (a-b)/a)
+     * @param bodyFrame body frame related to body shape
+     * @see org.orekit.frames.FramesFactory#getITRF(org.orekit.utils.IERSConventions, boolean)
+     */
+    public ExtendedEllipsoid(final double ae, final double f, final Frame bodyFrame) {
+        super(ae, f, bodyFrame);
+    }
+
+    /** Get point at some latitude along a pixel line of sight.
+     * @param position pixel position (in body frame)
+     * @param los pixel line-of-sight, not necessarily normalized (in body frame)
+     * @param latitude latitude with respect to ellipsoid
+     * @return point at altitude
+     * @exception RuggedException if no such point exists
+     */
+    public Vector3D pointAtLatitude(final Vector3D position, final Vector3D los, final double latitude)
+        throws RuggedException {
+
+        // find apex of iso-latitude cone, somewhere along polar axis
+        final GeodeticPoint groundPoint = new GeodeticPoint(latitude, 0, 0);
+        final Vector3D gpCartesian = transform(groundPoint);
+        final double r             = FastMath.sqrt(gpCartesian.getX() * gpCartesian.getX() +
+                                                   gpCartesian.getY() * gpCartesian.getY());
+        final Vector3D zenith      = groundPoint.getZenith();
+        final Vector3D apex        = new Vector3D(1, gpCartesian, -r / FastMath.cos(latitude), zenith);
+
+        // quadratic equation representing line intersection with iso-latitude cone
+        // a k² + 2 b k + c = 0
+        final Vector3D delta = position.subtract(apex);
+        final double zz2     = zenith.getZ() * zenith.getZ();
+        final double   a     = zz2 * los.getNormSq()                 - los.getZ()   * los.getZ();
+        final double   b     = zz2 * Vector3D.dotProduct(delta, los) - delta.getZ() * los.getZ();
+        final double   c     = zz2 * delta.getNormSq()               - delta.getZ() * delta.getZ();
+
+        // find the two intersections along the line
+        final double   b2    = b * b;
+        final double   ac    = a * c;
+        if (b2 < ac) {
+            throw new RuggedException(RuggedMessages.LINE_OF_SIGHT_NEVER_CROSSES_LATITUDE,
+                                      FastMath.toDegrees(latitude));
+        }
+        final double s  = FastMath.sqrt(b2 - ac);
+        final double k1 = (b > 0) ? -(s + b) / a : c / (s - b);
+        final double k2 = c / (a * k1);
+
+        // the quadratic equation may find spurious solutions in the wrong cone nappe
+        final boolean  k1IsOK = (delta.getZ() + k1 * los.getZ()) * latitude >= 0;
+        final boolean  k2IsOK = (delta.getZ() + k2 * los.getZ()) * latitude >= 0;
+        final double selectedK;
+        if (k1IsOK) {
+            if (k2IsOK) {
+                // both solutions are in the good nappe,
+                // we arbitrarily select the one closest to the initial position
+                selectedK = FastMath.abs(k1) <= FastMath.abs(k2) ? k1 : k2;
+            } else {
+                // only k1 is in the good nappe
+                selectedK = k1;
+            }
+        } else {  
+            if (k2IsOK) {
+                // only k2 is in the good nappe
+                selectedK = k2;
+            } else {
+                // both solutions are in the wrong nappe,
+                // there are no solutions
+                throw new RuggedException(RuggedMessages.LINE_OF_SIGHT_NEVER_CROSSES_LATITUDE,
+                                          FastMath.toDegrees(latitude));
+            }
+        }
+
+        // compute point
+        return new Vector3D(1, position, selectedK, los);
+
+    }
+
+    /** Get point at some longitude along a pixel line of sight.
+     * @param position pixel position (in body frame)
+     * @param los pixel line-of-sight, not necessarily normalized (in body frame)
+     * @param longitude longitude with respect to ellipsoid
+     * @return point at altitude
+     * @exception RuggedException if no such point exists
+     */
+    public Vector3D pointAtLongitude(final Vector3D position, final Vector3D los, final double longitude)
+        throws RuggedException {
+
+        // normal to meridian
+        final Vector3D normal = new Vector3D(-FastMath.sin(longitude), FastMath.cos(longitude), 0);
+        final double d = Vector3D.dotProduct(los, normal);
+        if (FastMath.abs(d) < 1.0e-12) {
+            throw new RuggedException(RuggedMessages.LINE_OF_SIGHT_NEVER_CROSSES_LONGITUDE,
+                                      FastMath.toDegrees(longitude));
+        }
+
+        // compute point
+       return new Vector3D(1, position, -Vector3D.dotProduct(position, normal) / d, los);
+
+    }
+
+    /** Get point at some altitude along a pixel line of sight.
+     * @param position pixel position (in body frame)
+     * @param los pixel line-of-sight, not necessarily normalized (in body frame)
+     * @param altitude altitude with respect to ellipsoid
+     * @return point at altitude
+     * @exception RuggedException if no such point exists (typically too negative altitude)
+     */
+    public Vector3D pointAtAltitude(final Vector3D position, final Vector3D los, final double altitude)
+        throws RuggedException {
+        try {
+
+            // point on line closest to origin
+            final double   los2   = los.getNormSq();
+            final double   dot    = Vector3D.dotProduct(position, los);
+            final double   k0     = -dot / los2;
+            final Vector3D close0 = new Vector3D(1, position, k0, los);
+
+            // very rough guess: if body is spherical, the desired point on line
+            // is at distance ae + altitude from origin
+            final double r        = getEquatorialRadius() + altitude;
+            final double delta2   = r * r - close0.getNormSq();
+            if (delta2 < 0) {
+                throw new RuggedException(RuggedMessages.LINE_OF_SIGHT_NEVER_CROSSES_ALTITUDE, altitude);
+            }
+            final double deltaK   = FastMath.sqrt(delta2 / los2);
+            final double k1       = k0 + deltaK;
+            final double k2       = k0 - deltaK;
+            double k              = (FastMath.abs(k1) <= FastMath.abs(k2)) ? k1 : k2;
+
+            // this loop generally converges in 3 iterations
+            for (int i = 0; i < 100; ++i) {
+
+                final Vector3D      point   = new Vector3D(1, position, k, los);
+                final GeodeticPoint gpK     = transform(point, getBodyFrame(), null);
+                final double        deltaH  = altitude - gpK.getAltitude();
+                if (FastMath.abs(deltaH) <= ALTITUDE_CONVERGENCE) {
+                    return point;
+                }
+
+                // improve the offset using linear ratio between
+                // altitude variation and displacement along line-of-sight
+                k += deltaH / Vector3D.dotProduct(gpK.getZenith(), los);
+
+            }
+
+            // this should never happen
+            throw new RuggedException(RuggedMessages.LINE_OF_SIGHT_NEVER_CROSSES_ALTITUDE, altitude);
+
+        } catch (OrekitException oe) {
+            // this should never happen
+            throw new RuggedException(oe, oe.getSpecifier(), oe.getParts());
+        }
+    }
+
+}

rugged-core/src/test/java/org/orekit/rugged/core/ExtendedEllipsoidTest.java

+/* Copyright 2013-2014 CS Systèmes d'Information
+ * Licensed to CS Systèmes d'Information (CS) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * CS licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.orekit.rugged.core;
+
+import java.io.File;
+import java.net.URISyntaxException;
+
+import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
+import org.apache.commons.math3.util.FastMath;
+import org.junit.After;
+import org.junit.Assert;
+import org.junit.Before;
+import org.junit.Test;
+import org.orekit.bodies.GeodeticPoint;
+import org.orekit.data.DataProvidersManager;
+import org.orekit.data.DirectoryCrawler;
+import org.orekit.errors.OrekitException;
+import org.orekit.frames.Frame;
+import org.orekit.frames.FramesFactory;
+import org.orekit.rugged.api.RuggedException;
+import org.orekit.rugged.api.RuggedMessages;
+import org.orekit.rugged.core.ExtendedEllipsoid;
+import org.orekit.utils.Constants;
+import org.orekit.utils.IERSConventions;
+
+public class ExtendedEllipsoidTest {
+
+    @Test
+    public void testPointAtLongitude() throws RuggedException, OrekitException {
+
+        Vector3D p = new Vector3D(3220103.0, 69623.0, -6449822.0);
+        Vector3D d = new Vector3D(1.0, 2.0, 3.0);
+
+        for (double longitude = -1.0; longitude < 1.0; longitude += 0.01) {
+            GeodeticPoint gp = ellipsoid.transform(ellipsoid.pointAtLongitude(p, d, longitude),
+                                                   ellipsoid.getBodyFrame(), null);
+            Assert.assertEquals(longitude, gp.getLongitude(), 1.0e-15);
+        }
+
+    }
+
+    @Test
+    public void testPointAtLongitudeError() throws RuggedException, OrekitException {
+
+        Vector3D p = new Vector3D(3220103.0, 69623.0, -6449822.0);
+        double longitude = 1.25;
+        Vector3D parallelToLongitudePlane = new Vector3D(FastMath.cos(longitude),
+                                                         FastMath.sin(longitude),
+                                                         -2.4);
+        try {
+            ellipsoid.pointAtLongitude(p, parallelToLongitudePlane, longitude);
+            Assert.fail("an error should have been triggered");
+        } catch (RuggedException re) {
+            Assert.assertEquals(RuggedMessages.LINE_OF_SIGHT_NEVER_CROSSES_LONGITUDE, re.getSpecifier());
+            Assert.assertEquals(FastMath.toDegrees(longitude), re.getParts()[0]);
+        }
+
+    }
+
+    @Test
+    public void testPointAtLatitude() throws RuggedException, OrekitException {
+
+        Vector3D p = new Vector3D(3220103.0, 69623.0, -6449822.0);
+        Vector3D d = new Vector3D(1.0, 2.0, 3.0);
+
+        for (double latitude = -d.getDelta() + 1.0e-6; latitude < d.getDelta(); latitude += 0.1) {
+            GeodeticPoint gp = ellipsoid.transform(ellipsoid.pointAtLatitude(p, d, latitude),
+                                                   ellipsoid.getBodyFrame(), null);
+            Assert.assertEquals(latitude, gp.getLatitude(), 1.0e-12);
+        }
+
+    }
+
+    @Test
+    public void testPointAtLatitudeTwoPointsInGoodNappe() throws RuggedException, OrekitException {
+
+        Vector3D p = new Vector3D(3220103.0, 69623.0, -6449822.0);
+        Vector3D d = new Vector3D(1.0, 2.0, 0.1);
+        double latitude = -0.5;
+
+        GeodeticPoint gpPlus = ellipsoid.transform(ellipsoid.pointAtLatitude(p, d, latitude),
+                                                   ellipsoid.getBodyFrame(), null);
+        Assert.assertEquals(latitude, gpPlus.getLatitude(), 1.0e-12);
+
+        GeodeticPoint gpMinus = ellipsoid.transform(ellipsoid.pointAtLatitude(p, d.negate(), latitude),
+                                                    ellipsoid.getBodyFrame(), null);
+        Assert.assertEquals(latitude, gpMinus.getLatitude(), 1.0e-12);
+
+    }
+
+    @Test
+    public void testPointAtLatitudeErrorQuadraticEquation() throws RuggedException, OrekitException {
+
+        Vector3D p = new Vector3D(3220103.0, 69623.0, -6449822.0);
+        Vector3D d = new Vector3D(1.0, 2.0, 3.0);
+        double latitude = -1.4;
+
+        try {
+            ellipsoid.pointAtLatitude(p, d, latitude);
+            Assert.fail("an error should have been triggered");
+        } catch (RuggedException re) {
+            Assert.assertEquals(RuggedMessages.LINE_OF_SIGHT_NEVER_CROSSES_LATITUDE, re.getSpecifier());
+            Assert.assertEquals(FastMath.toDegrees(latitude), re.getParts()[0]);
+        }
+
+    }
+
+    @Test
+    public void testPointAtLatitudeErrorNappe() throws RuggedException, OrekitException {
+
+        Vector3D p = new Vector3D(3220103.0, 69623.0, -6449822.0);
+        Vector3D d = new Vector3D(1.0, 2.0, 0.1);
+        double latitude = 0.5;
+
+        try {
+            ellipsoid.pointAtLatitude(p, d, latitude);
+            Assert.fail("an error should have been triggered");
+        } catch (RuggedException re) {
+            Assert.assertEquals(RuggedMessages.LINE_OF_SIGHT_NEVER_CROSSES_LATITUDE, re.getSpecifier());
+            Assert.assertEquals(FastMath.toDegrees(latitude), re.getParts()[0]);
+        }
+
+    }
+
+    @Test
+    public void testPointAtAltitude() throws RuggedException, OrekitException {
+
+        Vector3D p = new Vector3D(3220103.0, 69623.0, -6449822.0);
+        Vector3D d = new Vector3D(1.0, 2.0, 3.0);
+        for (double altitude = -500000; altitude < 800000.0; altitude += 100) {
+            GeodeticPoint gp = ellipsoid.transform(ellipsoid.pointAtAltitude(p, d, altitude),
+                                                   ellipsoid.getBodyFrame(), null);
+            Assert.assertEquals(altitude, gp.getAltitude(), 1.0e-3);
+        }
+
+    }
+
+    @Test
+    public void testPointAtAltitudeStartInside() throws RuggedException, OrekitException {
+
+        Vector3D p = new Vector3D(322010.30, 6962.30, -644982.20);
+        Vector3D d = new Vector3D(-1.0, -2.0, -3.0);
+        for (double altitude = -500000; altitude < 800000.0; altitude += 100) {
+            GeodeticPoint gp = ellipsoid.transform(ellipsoid.pointAtAltitude(p, d, altitude),
+                                                   ellipsoid.getBodyFrame(), null);
+            Assert.assertEquals(altitude, gp.getAltitude(), 1.0e-3);
+        }
+
+    }
+
+    @Test
+    public void testPointAtAltitudeError() throws RuggedException, OrekitException {
+
+        Vector3D p = new Vector3D(3220103.0, 69623.0, -6449822.0);
+        Vector3D d = new Vector3D(1.0, 2.0, 3.0);
+        double altitude = -580000.0;
+        try {
+            ellipsoid.pointAtAltitude(p, d, altitude);
+            Assert.fail("an error should have been triggered");
+        } catch (RuggedException re) {
+            Assert.assertEquals(RuggedMessages.LINE_OF_SIGHT_NEVER_CROSSES_ALTITUDE, re.getSpecifier());
+            Assert.assertEquals(altitude, re.getParts()[0]);
+        }
+
+    }
+
+    @Before
+    public void setUp() {
+        try {
+
+            String path = getClass().getClassLoader().getResource("orekit-data").toURI().getPath();
+            DataProvidersManager.getInstance().addProvider(new DirectoryCrawler(new File(path)));
+
+            Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
+            ellipsoid = new ExtendedEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
+                                              Constants.WGS84_EARTH_FLATTENING,
+                                              itrf);
+        } catch (OrekitException oe) {
+            Assert.fail(oe.getLocalizedMessage());
+        } catch (URISyntaxException use) {
+            Assert.fail(use.getLocalizedMessage());
+        }
+    }
+
+    @After
+    public void tearDown() {
+        ellipsoid = null;
+    }
+
+    private ExtendedEllipsoid ellipsoid;
+
+}

rugged-core/src/test/resources/orekit-data/UTC-TAI.history

+
+ ---------------
+ UTC-TAI.history
+ ---------------
+ RELATIONSHIP BETWEEN TAI AND UTC
+ ------------------------------------------------------------------------------- 
+ Limits of validity(at 0h UTC)       TAI - UTC  
+ 
+ 1961  Jan.  1 - 1961  Aug.  1     1.422 818 0s + (MJD - 37 300) x 0.001 296s
+       Aug.  1 - 1962  Jan.  1     1.372 818 0s +        ""
+ 1962  Jan.  1 - 1963  Nov.  1     1.845 858 0s + (MJD - 37 665) x 0.001 123 2s
+ 1963  Nov.  1 - 1964  Jan.  1     1.945 858 0s +        ""
+ 1964  Jan.  1 -       April 1     3.240 130 0s + (MJD - 38 761) x 0.001 296s
+       April 1 -       Sept. 1     3.340 130 0s +        ""
+       Sept. 1 - 1965  Jan.  1     3.440 130 0s +        ""
+ 1965  Jan.  1 -       March 1     3.540 130 0s +        ""
+       March 1 -       Jul.  1     3.640 130 0s +        ""
+       Jul.  1 -       Sept. 1     3.740 130 0s +        ""
+       Sept. 1 - 1966  Jan.  1     3.840 130 0s +        ""
+ 1966  Jan.  1 - 1968  Feb.  1     4.313 170 0s + (MJD - 39 126) x 0.002 592s
+ 1968  Feb.  1 - 1972  Jan.  1     4.213 170 0s +        ""
+ 1972  Jan.  1 -       Jul.  1    10s            
+       Jul.  1 - 1973  Jan.  1    11s
+ 1973  Jan.  1 - 1974  Jan.  1    12s	
+ 1974  Jan.  1 - 1975  Jan.  1    13s	
+ 1975  Jan.  1 - 1976  Jan.  1    14s	  
+ 1976  Jan.  1 - 1977  Jan.  1    15s	    
+ 1977  Jan.  1 - 1978  Jan.  1    16s
+ 1978  Jan.  1 - 1979  Jan.  1    17s	
+ 1979  Jan.  1 - 1980  Jan.  1    18s	
+ 1980  Jan.  1 - 1981  Jul.  1    19s	
+ 1981  Jul.  1 - 1982  Jul.  1    20s	
+ 1982  Jul.  1 - 1983  Jul.  1    21s
+ 1983  Jul.  1 - 1985  Jul.  1    22s
+ 1985  Jul.  1 - 1988  Jan.  1    23s
+ 1988  Jan.  1 - 1990  Jan.  1    24s
+ 1990  Jan.  1 - 1991  Jan.  1    25s
+ 1991  Jan.  1 - 1992  Jul.  1    26s
+ 1992  Jul.  1.- 1993  Jul   1    27s
+ 1993  Jul.  1 - 1994  Jul.  1    28s
+ 1994  Jul.  1 - 1996  Jan.  1    29s
+ 1996  Jan.  1 - 1997  Jul.  1    30s
+ 1997  Jul.  1.- 1999  Jan.  1    31s
+ 1999  Jan.  1.- 2006  Jan.  1    32s
+ 2006  Jan.  1.- 2009  Jan.  1    33s
+ 2009  Jan.  1.- 2012  Jul   1    34s
+ 2012  Jul   1 -                  35s
+ ----------------------------------------------------------------------