From 046af11e3e816e20371184865d8939c34c14715e Mon Sep 17 00:00:00 2001
From: Luc Maisonobe <luc@orekit.org>
Date: Wed, 15 Jul 2015 12:07:38 +0200
Subject: [PATCH] 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.
---
.../java/org/orekit/rugged/api/Rugged.java | 56 ++++++
.../utils/RoughVisibilityEstimator.java | 183 ++++++++++++++++++
src/site/xdoc/changes.xml | 3 +
.../utils/RoughVisibilityEstimatorTest.java | 182 +++++++++++++++++
4 files changed, 424 insertions(+)
create mode 100644 src/main/java/org/orekit/rugged/utils/RoughVisibilityEstimator.java
create mode 100644 src/test/java/org/orekit/rugged/utils/RoughVisibilityEstimatorTest.java
diff --git a/src/main/java/org/orekit/rugged/api/Rugged.java b/src/main/java/org/orekit/rugged/api/Rugged.java
index fe9b8dcf..ba8554ae 100644
--- a/src/main/java/org/orekit/rugged/api/Rugged.java
+++ b/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)
diff --git a/src/main/java/org/orekit/rugged/utils/RoughVisibilityEstimator.java b/src/main/java/org/orekit/rugged/utils/RoughVisibilityEstimator.java
new file mode 100644
index 00000000..cfd7d055
--- /dev/null
+++ b/src/main/java/org/orekit/rugged/utils/RoughVisibilityEstimator.java
@@ -0,0 +1,183 @@
+/* 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);
+
+ }
+
+}
diff --git a/src/site/xdoc/changes.xml b/src/site/xdoc/changes.xml
index 85943e20..97deeea4 100644
--- a/src/site/xdoc/changes.xml
+++ b/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>
diff --git a/src/test/java/org/orekit/rugged/utils/RoughVisibilityEstimatorTest.java b/src/test/java/org/orekit/rugged/utils/RoughVisibilityEstimatorTest.java
new file mode 100644
index 00000000..2245083c
--- /dev/null
+++ b/src/test/java/org/orekit/rugged/utils/RoughVisibilityEstimatorTest.java
@@ -0,0 +1,182 @@
+/* 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;
+
+}
--
GitLab