Added a utility to help estimate visibility.

This can be used to set up min/max search lines.

This is still preliminary work, convergence seems slightly wrong.

src/main/java/org/orekit/rugged/api/Rugged.java

@@ -348,6 +348,19 @@ public class Rugged {
      * The point is given only by its latitude and longitude, the elevation is
      * computed from the Digital Elevation Model.
      * </p>
+     * <p>
+     * Note that for each sensor name, the {@code minLine} and {@code maxLine} settings
+     * are cached, because they induce costly frames computation. So these settings
+     * should not be tuned very finely and changed at each call, but should rather be
+     * a few thousand lines wide and refreshed only when needed. If for example an
+     * inverse location is roughly estimated to occur near line 53764 (for example
+     * using {@link org.orekit.rugged.utils.RoughVisibilityEstimator), {@code minLine}
+     * and {@code maxLine} could be set for example to 50000 and 60000, which would
+     * be OK also if next line inverse location is expected to occur near line 53780,
+     * and next one ... The setting could be changed for example to 55000 and 65000 when
+     * an inverse location is expected to occur after 55750. Of course, these values
+     * are only an example and should be adjusted depending on mission needs.
+     * </p>
      * @param sensorName name of the line  sensor
      * @param latitude ground point latitude
      * @param longitude ground point longitude
@@ -356,6 +369,7 @@ public class Rugged {
      * @return date at which ground point is seen by line sensor
      * @exception RuggedException if line cannot be localized, or sensor is unknown
      * @see #inverseLocation(String, double, double, int, int)
+     * @see org.orekit.rugged.utils.RoughVisibilityEstimator
      */
     public AbsoluteDate dateLocation(final String sensorName,
                                      final double latitude, final double longitude,
@@ -371,6 +385,19 @@ public class Rugged {
      * This method is a partial {@link #inverseLocation(String,
      * GeodeticPoint, int, int) inverse location} focusing only on date.
      * </p>
+     * <p>
+     * Note that for each sensor name, the {@code minLine} and {@code maxLine} settings
+     * are cached, because they induce costly frames computation. So these settings
+     * should not be tuned very finely and changed at each call, but should rather be
+     * a few thousand lines wide and refreshed only when needed. If for example an
+     * inverse location is roughly estimated to occur near line 53764 (for example
+     * using {@link org.orekit.rugged.utils.RoughVisibilityEstimator), {@code minLine}
+     * and {@code maxLine} could be set for example to 50000 and 60000, which would
+     * be OK also if next line inverse location is expected to occur near line 53780,
+     * and next one ... The setting could be changed for example to 55000 and 65000 when
+     * an inverse location is expected to occur after 55750. Of course, these values
+     * are only an example and should be adjusted depending on mission needs.
+     * </p>
      * @param sensorName name of the line  sensor
      * @param point point to localize
      * @param minLine minimum line number
@@ -378,6 +405,7 @@ public class Rugged {
      * @return date at which ground point is seen by line sensor
      * @exception RuggedException if line cannot be localized, or sensor is unknown
      * @see #inverseLocation(String, GeodeticPoint, int, int)
+     * @see org.orekit.rugged.utils.RoughVisibilityEstimator
      */
     public AbsoluteDate dateLocation(final String sensorName, final GeodeticPoint point,
                                      final int minLine, final int maxLine)
@@ -403,6 +431,19 @@ public class Rugged {
      * The point is given only by its latitude and longitude, the elevation is
      * computed from the Digital Elevation Model.
      * </p>
+     * <p>
+     * Note that for each sensor name, the {@code minLine} and {@code maxLine} settings
+     * are cached, because they induce costly frames computation. So these settings
+     * should not be tuned very finely and changed at each call, but should rather be
+     * a few thousand lines wide and refreshed only when needed. If for example an
+     * inverse location is roughly estimated to occur near line 53764 (for example
+     * using {@link org.orekit.rugged.utils.RoughVisibilityEstimator), {@code minLine}
+     * and {@code maxLine} could be set for example to 50000 and 60000, which would
+     * be OK also if next line inverse location is expected to occur near line 53780,
+     * and next one ... The setting could be changed for example to 55000 and 65000 when
+     * an inverse location is expected to occur after 55750. Of course, these values
+     * are only an example and should be adjusted depending on mission needs.
+     * </p>
      * @param sensorName name of the line  sensor
      * @param latitude ground point latitude
      * @param longitude ground point longitude
@@ -411,6 +452,7 @@ public class Rugged {
      * @return sensor pixel seeing ground point, or null if ground point cannot
      * be seen between the prescribed line numbers
      * @exception RuggedException if line cannot be localized, or sensor is unknown
+     * @see org.orekit.rugged.utils.RoughVisibilityEstimator
      */
     public SensorPixel inverseLocation(final String sensorName,
                                        final double latitude, final double longitude,
@@ -422,6 +464,19 @@ public class Rugged {
     }
 
     /** Inverse location of a point.
+     * <p>
+     * Note that for each sensor name, the {@code minLine} and {@code maxLine} settings
+     * are cached, because they induce costly frames computation. So these settings
+     * should not be tuned very finely and changed at each call, but should rather be
+     * a few thousand lines wide and refreshed only when needed. If for example an
+     * inverse location is roughly estimated to occur near line 53764 (for example
+     * using {@link org.orekit.rugged.utils.RoughVisibilityEstimator), {@code minLine}
+     * and {@code maxLine} could be set for example to 50000 and 60000, which would
+     * be OK also if next line inverse location is expected to occur near line 53780,
+     * and next one ... The setting could be changed for example to 55000 and 65000 when
+     * an inverse location is expected to occur after 55750. Of course, these values
+     * are only an example and should be adjusted depending on mission needs.
+     * </p>
      * @param sensorName name of the line  sensor
      * @param point point to localize
      * @param minLine minimum line number
@@ -430,6 +485,7 @@ public class Rugged {
      * prescribed line numbers
      * @exception RuggedException if line cannot be localized, or sensor is unknown
      * @see #dateLocation(String, GeodeticPoint, int, int)
+     * @see org.orekit.rugged.utils.RoughVisibilityEstimator
      */
     public SensorPixel inverseLocation(final String sensorName, final GeodeticPoint point,
                                        final int minLine, final int maxLine)

src/main/java/org/orekit/rugged/utils/RoughVisibilityEstimator.java

+/* Copyright 2013-2015 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.utils;
+
+import java.util.ArrayList;
+import java.util.List;
+
+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.frames.Transform;
+import org.orekit.time.AbsoluteDate;
+import org.orekit.utils.TimeStampedPVCoordinates;
+
+/** Class estimating very roughly when a point may be visible from spacecraft.
+ * <p>
+ * The class only uses spacecraft position to compute a very rough sub-satellite
+ * point. It assumes the position-velocities are regular enough and without holes.
+ * It is intended only only has a quick estimation in order to set up search
+ * boundaries in inverse location.
+ * </p>
+ * @see org.orekit.rugged.api.Rugged#dateLocation(String, org.orekit.bodies.GeodeticPoint, int, int)
+ * @see org.orekit.rugged.api.Rugged#dateLocation(String, double, double, int, int)
+ * @see org.orekit.rugged.api.Rugged#inverseLocation(String, org.orekit.bodies.GeodeticPoint, int, int)
+ * @see org.orekit.rugged.api.Rugged#inverseLocation(String, double, double, int, int)
+ * @author Luc Maisonobe
+ */
+public class RoughVisibilityEstimator {
+
+    /** Ground ellipsoid. */
+    private final OneAxisEllipsoid ellipsoid;
+
+    /** Sub-satellite point. */
+    private final List<TimeStampedPVCoordinates> pvGround;
+
+    /** Mean angular rate with respect to position-velocity indices. */
+    private final double rateVSIndices;
+
+    /** Mean angular rate with respect to time. */
+    private final double rateVSTime;
+
+    /** Last found index. */
+    private int last;
+
+    /**
+     * Simple constructor.
+     * @param ellipsoid ground ellipsoid
+     * @param frame frame in which position and velocity are defined (may be inertial or body frame)
+     * @param positionsVelocities satellite position and velocity (m and m/s in specified frame)
+     * @exception OrekitException if position-velocity cannot be converted to body frame
+     */
+    public RoughVisibilityEstimator(final OneAxisEllipsoid ellipsoid, final Frame frame,
+                                    final List<TimeStampedPVCoordinates> positionsVelocities)
+        throws OrekitException {
+
+        this.ellipsoid = ellipsoid;
+
+        // project spacecraft position-velocity to ground
+        final Frame bodyFrame = ellipsoid.getBodyFrame();
+        final int n = positionsVelocities.size();
+        this.pvGround = new ArrayList<TimeStampedPVCoordinates>(n);
+        for (final TimeStampedPVCoordinates pv : positionsVelocities) {
+            final Transform t = frame.getTransformTo(bodyFrame, pv.getDate());
+            pvGround.add(ellipsoid.projectToGround(t.transformPVCoordinates(pv), bodyFrame));
+        }
+
+        // initialize first search at mid point
+        this.last = n / 2;
+
+        // estimate mean angular rate with respect to indices
+        double alpha = 0;
+        for (int i = 0; i < n - 1; ++i) {
+            // angular motion between points i and i+1
+            alpha += Vector3D.angle(pvGround.get(i).getPosition(),
+                                    pvGround.get(i + 1).getPosition());
+        }
+        this.rateVSIndices = alpha / n;
+
+        // estimate mean angular rate with respect to time
+        final AbsoluteDate firstDate = pvGround.get(0).getDate();
+        final AbsoluteDate lastDate  = pvGround.get(pvGround.size() - 1).getDate();
+        this.rateVSTime              = alpha / lastDate.durationFrom(firstDate);
+
+    }
+
+    /** Estimate <em>very roughly</em> when spacecraft comes close to a ground point.
+     * @param groundPoint ground point to check
+     * @return rough date at which spacecraft comes close to ground point (never null,
+     * but may be really far from reality if ground point is away from trajectory)
+     */
+    public AbsoluteDate estimateVisibility(final GeodeticPoint groundPoint) {
+
+        final Vector3D point = ellipsoid.transform(groundPoint);
+        int closeIndex = findClose(last, point);
+
+        // check if there are closer points in previous half periods
+        final int repeat = (int) FastMath.rint(FastMath.PI / rateVSIndices);
+        for (int index = closeIndex - repeat; index > 0; index -= repeat) {
+            final int otherIndex = findClose(index, point);
+            if (otherIndex != closeIndex &&
+                Vector3D.distance(pvGround.get(otherIndex).getPosition(), point) <
+                Vector3D.distance(pvGround.get(closeIndex).getPosition(), point)) {
+                closeIndex = otherIndex;
+            }
+        }
+
+        // check if there are closer points in next half periods
+        for (int index = closeIndex + repeat; index < pvGround.size(); index += repeat) {
+            final int otherIndex = findClose(index, point);
+            if (otherIndex != closeIndex &&
+                Vector3D.distance(pvGround.get(otherIndex).getPosition(), point) <
+                Vector3D.distance(pvGround.get(closeIndex).getPosition(), point)) {
+                closeIndex = otherIndex;
+            }
+        }
+
+        // we have found the closest sub-satellite point index
+        last = closeIndex;
+
+        // final adjustment
+        final TimeStampedPVCoordinates closest = pvGround.get(closeIndex);
+        final double alpha = neededMotion(closest, point);
+        return closest.getDate().shiftedBy(alpha / rateVSTime);
+
+    }
+
+    /** Find the index of a close sub-satellite point.
+     * @param start start index for the search
+     * @param point test point
+     * @return index of a sub-satellite point close to the test point
+     */
+    private int findClose(final int start, final Vector3D point) {
+        int current  = start;
+        int previous = Integer.MIN_VALUE;
+        int maxLoop  = 1000;
+        while (maxLoop-- > 0 && FastMath.abs(current - previous) > 1) {
+            previous = current;
+            final double alpha = neededMotion(pvGround.get(current), point);
+            final int    shift = (int) FastMath.rint(alpha / rateVSIndices);
+            current = FastMath.max(0, FastMath.min(pvGround.size() - 1, current + shift));
+        }
+        return current;
+    }
+
+    /** Estimate angular motion needed to go past test point.
+     * @param subSatellite current sub-satellite position-velocity
+     * @param point test point
+     * @return angular motion to go past test point (positive is
+     * test point is ahead, negative if it is behind)
+     */
+    private double neededMotion(final TimeStampedPVCoordinates subSatellite,
+                                final Vector3D point) {
+
+        final Vector3D ssP      = subSatellite.getPosition();
+        final Vector3D momentum = subSatellite.getMomentum();
+
+        // find phase angle around momentum
+        final Vector3D yAxis = Vector3D.crossProduct(momentum, ssP).normalize();
+        final Vector3D xAxis = Vector3D.crossProduct(yAxis, momentum).normalize();
+        final double   y     = Vector3D.dotProduct(point, yAxis);
+        final double   x     = Vector3D.dotProduct(point, xAxis);
+        return FastMath.atan2(y, x);
+
+    }
+
+}

src/site/xdoc/changes.xml

@@ -22,6 +22,9 @@
   <body>
     <release version="1.0" date="TBD"
              description="TBD">
+      <action dev="luc" type="add" >
+        Added a utility to help estimate visibility and set up min/max search lines.
+      </action>
       <action dev="luc" type="fix" >
         Improved latitude crossing robustness.
       </action>

src/test/java/org/orekit/rugged/utils/RoughVisibilityEstimatorTest.java

+/* Copyright 2013-2015 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.utils;
+
+import java.io.File;
+import java.net.URISyntaxException;
+import java.util.ArrayList;
+import java.util.List;
+
+import org.apache.commons.math3.ode.nonstiff.DormandPrince853Integrator;
+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.attitudes.AttitudeProvider;
+import org.orekit.attitudes.NadirPointing;
+import org.orekit.attitudes.YawCompensation;
+import org.orekit.bodies.BodyShape;
+import org.orekit.bodies.CelestialBodyFactory;
+import org.orekit.bodies.GeodeticPoint;
+import org.orekit.bodies.OneAxisEllipsoid;
+import org.orekit.data.DataProvidersManager;
+import org.orekit.data.DirectoryCrawler;
+import org.orekit.errors.OrekitException;
+import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
+import org.orekit.forces.gravity.ThirdBodyAttraction;
+import org.orekit.forces.gravity.potential.GravityFieldFactory;
+import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
+import org.orekit.frames.Frame;
+import org.orekit.frames.FramesFactory;
+import org.orekit.orbits.CircularOrbit;
+import org.orekit.orbits.Orbit;
+import org.orekit.orbits.OrbitType;
+import org.orekit.orbits.PositionAngle;
+import org.orekit.propagation.Propagator;
+import org.orekit.propagation.SpacecraftState;
+import org.orekit.propagation.numerical.NumericalPropagator;
+import org.orekit.propagation.sampling.OrekitFixedStepHandler;
+import org.orekit.rugged.errors.RuggedException;
+import org.orekit.time.AbsoluteDate;
+import org.orekit.time.TimeScalesFactory;
+import org.orekit.utils.Constants;
+import org.orekit.utils.IERSConventions;
+import org.orekit.utils.TimeStampedPVCoordinates;
+
+public class RoughVisibilityEstimatorTest {
+
+    @Test
+    public void testThreeOrbitsSpan()
+        throws RuggedException, OrekitException, URISyntaxException {
+
+        String path = getClass().getClassLoader().getResource("orekit-data").toURI().getPath();
+        DataProvidersManager.getInstance().addProvider(new DirectoryCrawler(new File(path)));
+        BodyShape  earth                                  = createEarth();
+        NormalizedSphericalHarmonicsProvider gravityField = createGravityField();
+        Orbit      orbit                                  = createOrbit(gravityField.getMu());
+        Propagator propagator                             = createPropagator(earth, gravityField, orbit);
+        final List<TimeStampedPVCoordinates> pv = new ArrayList<TimeStampedPVCoordinates>();
+        propagator.setMasterMode(1.0, new OrekitFixedStepHandler() {
+            public void init(SpacecraftState s0, AbsoluteDate t) {
+            }
+            public void handleStep(SpacecraftState currentState, boolean isLast) {
+                pv.add(currentState.getPVCoordinates());
+            }
+        });
+        propagator.propagate(orbit.getDate().shiftedBy(3 * orbit.getKeplerianPeriod()));
+
+        RoughVisibilityEstimator estimator = new RoughVisibilityEstimator(ellipsoid, orbit.getFrame(), pv);
+        AbsoluteDate d = estimator.estimateVisibility(new GeodeticPoint(FastMath.toRadians(-81.5),
+                                                                        FastMath.toRadians(-2.0),
+                                                                        0.0));
+        Assert.assertEquals("2012-01-01T03:47:08.814", d.toString(TimeScalesFactory.getUTC()));
+
+    }
+
+    private BodyShape createEarth()
+                    throws OrekitException {
+        return new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
+                                    Constants.WGS84_EARTH_FLATTENING,
+                                    FramesFactory.getITRF(IERSConventions.IERS_2010, true));
+    }
+
+    private NormalizedSphericalHarmonicsProvider createGravityField()
+                    throws OrekitException {
+        return GravityFieldFactory.getNormalizedProvider(12, 12);
+    }
+
+    private Orbit createOrbit(double mu)
+                    throws OrekitException {
+        // the following orbital parameters have been computed using
+        // Orekit tutorial about phasing, using the following configuration:
+        //
+        //  orbit.date                          = 2012-01-01T00:00:00.000
+        //  phasing.orbits.number               = 143
+        //  phasing.days.number                 =  10
+        //  sun.synchronous.reference.latitude  = 0
+        //  sun.synchronous.reference.ascending = false
+        //  sun.synchronous.mean.solar.time     = 10:30:00
+        //  gravity.field.degree                = 12
+        //  gravity.field.order                 = 12
+        AbsoluteDate date = new AbsoluteDate("2012-01-01T00:00:00.000", TimeScalesFactory.getUTC());
+        Frame eme2000 = FramesFactory.getEME2000();
+        return new CircularOrbit(7173352.811913891,
+                                 -4.029194321683225E-4, 0.0013530362644647786,
+                                 FastMath.toRadians(98.63218182243709),
+                                 FastMath.toRadians(77.55565567747836),
+                                 FastMath.PI, PositionAngle.TRUE,
+                                 eme2000, date, mu);
+    }
+
+    private Propagator createPropagator(BodyShape earth,
+                                        NormalizedSphericalHarmonicsProvider gravityField,
+                                        Orbit orbit)
+                                                        throws OrekitException {
+
+        AttitudeProvider yawCompensation = new YawCompensation(orbit.getFrame(), new NadirPointing(orbit.getFrame(), earth));
+        SpacecraftState state = new SpacecraftState(orbit,
+                                                    yawCompensation.getAttitude(orbit,
+                                                                                orbit.getDate(),
+                                                                                orbit.getFrame()),
+                                                    1180.0);
+
+        // numerical model for improving orbit
+        OrbitType type = OrbitType.CIRCULAR;
+        double[][] tolerances = NumericalPropagator.tolerances(0.1, orbit, type);
+        DormandPrince853Integrator integrator =
+                        new DormandPrince853Integrator(1.0e-4 * orbit.getKeplerianPeriod(),
+                                                       1.0e-1 * orbit.getKeplerianPeriod(),
+                                                       tolerances[0], tolerances[1]);
+        integrator.setInitialStepSize(1.0e-2 * orbit.getKeplerianPeriod());
+        NumericalPropagator numericalPropagator = new NumericalPropagator(integrator);
+        numericalPropagator.addForceModel(new HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), gravityField));
+        numericalPropagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
+        numericalPropagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
+        numericalPropagator.setOrbitType(type);
+        numericalPropagator.setInitialState(state);
+        numericalPropagator.setAttitudeProvider(yawCompensation);
+        return numericalPropagator;
+
+    }
+
+    @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;
+
+}