diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index 4d0f90cc2dd1b7e7012f79fce19cc0cd8d1bd50b..4750eb52f912aee992a8ce8dfd4fa205b589494d 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -21,6 +21,9 @@
</properties>
<body>
<release version="3.1" date="TBD" description="TBD">
+ <action dev="guylaine" type="update" issue="392">
+ Fixed light time correction in atmospheric refraction computation.
+ </action>
</release>
<release version="3.0" date="2022-07-05" description="This is a major release.
It fixes a few bugs.
diff --git a/src/main/java/org/orekit/rugged/api/Rugged.java b/src/main/java/org/orekit/rugged/api/Rugged.java
index 32ae412739ea5c9dd30036aa299a99938884529d..6a89830acbc5d91621b9351e8bc3d913fa2a5e4b 100644
--- a/src/main/java/org/orekit/rugged/api/Rugged.java
+++ b/src/main/java/org/orekit/rugged/api/Rugged.java
@@ -300,10 +300,35 @@ public class Rugged {
// compute with atmospheric refraction correction if necessary
if (atmosphericRefraction != null && atmosphericRefraction.mustBeComputed()) {
+ final Vector3D pBody;
+ final Vector3D lBody;
+
+ // Take into account the light time correction
+ // @since 3.1
+ if (lightTimeCorrection) {
+ // Transform sensor position in inertial frame to observed body
+ final Vector3D sP = inertToBody.transformPosition(pInert);
+ // Convert ground location of the pixel in cartesian coordinates
+ final Vector3D eP = ellipsoid.transform(gp[i]);
+ // Compute the light time correction (s)
+ final double deltaT = eP.distance(sP) / Constants.SPEED_OF_LIGHT;
+
+ // Apply shift due to light time correction
+ final Transform shiftedInertToBody = inertToBody.shiftedBy(-deltaT);
+
+ pBody = shiftedInertToBody.transformPosition(pInert);
+ lBody = shiftedInertToBody.transformVector(lInert);
+
+ } else { // Light time correction NOT to be taken into account
+
+ pBody = inertToBody.transformPosition(pInert);
+ lBody = inertToBody.transformVector(lInert);
+
+ } // end test on lightTimeCorrection
+
// apply atmospheric refraction correction
- final Vector3D pBody = inertToBody.transformPosition(pInert);
- final Vector3D lBody = inertToBody.transformVector(lInert);
gp[i] = atmosphericRefraction.applyCorrection(pBody, lBody, (NormalizedGeodeticPoint) gp[i], algorithm);
+
}
DumpManager.dumpDirectLocationResult(gp[i]);
}
@@ -346,7 +371,7 @@ public class Rugged {
lInert = obsLInert;
}
- // Compute ground location of specified pixel
+ // Compute ground location of specified pixel according to light time correction flag
final NormalizedGeodeticPoint gp;
if (lightTimeCorrection) {
@@ -362,14 +387,39 @@ public class Rugged {
algorithm.intersection(ellipsoid, pBody, lBody));
}
+ // Compute ground location of specified pixel according to atmospheric refraction correction flag
NormalizedGeodeticPoint result = gp;
// compute the ground location with atmospheric correction if asked for
if (atmosphericRefraction != null && atmosphericRefraction.mustBeComputed()) {
+ final Vector3D pBody;
+ final Vector3D lBody;
+
+ // Take into account the light time correction
+ // @since 3.1
+ if (lightTimeCorrection) {
+ // Transform sensor position in inertial frame to observed body
+ final Vector3D sP = inertToBody.transformPosition(pInert);
+ // Convert ground location of the pixel in cartesian coordinates
+ final Vector3D eP = ellipsoid.transform(gp);
+ // Compute the light time correction (s)
+ final double deltaT = eP.distance(sP) / Constants.SPEED_OF_LIGHT;
+
+ // Apply shift due to light time correction
+ final Transform shiftedInertToBody = inertToBody.shiftedBy(-deltaT);
+
+ pBody = shiftedInertToBody.transformPosition(pInert);
+ lBody = shiftedInertToBody.transformVector(lInert);
+
+ } else { // Light time correction NOT to be taken into account
+
+ pBody = inertToBody.transformPosition(pInert);
+ lBody = inertToBody.transformVector(lInert);
+
+ } // end test on lightTimeCorrection
+
// apply atmospheric refraction correction
- final Vector3D pBody = inertToBody.transformPosition(pInert);
- final Vector3D lBody = inertToBody.transformVector(lInert);
result = atmosphericRefraction.applyCorrection(pBody, lBody, gp, algorithm);
} // end test on atmosphericRefraction != null
@@ -569,7 +619,7 @@ public class Rugged {
* @param scToInert transform for the date from spacecraft to inertial
* @param inertToBody transform for the date from inertial to body
* @param pInert sensor position in inertial frame
- * @param lInert line of sight in inertial frame
+ * @param lInert line of sight in inertial frame (with light time correction if asked for)
* @return geodetic point with light time correction
*/
private NormalizedGeodeticPoint computeWithLightTimeCorrection(final AbsoluteDate date,
@@ -577,28 +627,41 @@ public class Rugged {
final Transform scToInert, final Transform inertToBody,
final Vector3D pInert, final Vector3D lInert) {
- // compute the approximate transform between spacecraft and observed body
+ // Compute the transform between spacecraft and observed body
final Transform approximate = new Transform(date, scToInert, inertToBody);
+ // Transform LOS in spacecraft frame to observed body
final Vector3D sL = approximate.transformVector(los);
+ // Transform sensor position in spacecraft frame to observed body
final Vector3D sP = approximate.transformPosition(sensorPosition);
+ // Compute point intersecting ground (= the ellipsoid) along the pixel LOS
final Vector3D eP1 = ellipsoid.transform(ellipsoid.pointOnGround(sP, sL, 0.0));
+
+ // Compute light time time correction (vs the ellipsoid) (s)
final double deltaT1 = eP1.distance(sP) / Constants.SPEED_OF_LIGHT;
+
+ // Apply shift due to light time correction (vs the ellipsoid)
final Transform shifted1 = inertToBody.shiftedBy(-deltaT1);
+
+ // Search the intersection of LOS (taking into account the light time correction if asked for) with DEM
final NormalizedGeodeticPoint gp1 = algorithm.intersection(ellipsoid,
shifted1.transformPosition(pInert),
shifted1.transformVector(lInert));
+ // Convert the geodetic point (intersection of LOS with DEM) in cartesian coordinates
final Vector3D eP2 = ellipsoid.transform(gp1);
+
+ // Compute the light time correction (vs DEM) (s)
final double deltaT2 = eP2.distance(sP) / Constants.SPEED_OF_LIGHT;
+
+ // Apply shift due to light time correction (vs DEM)
final Transform shifted2 = inertToBody.shiftedBy(-deltaT2);
- // At this stage the light time correction is negligible so not taken into account even if asked for
- // For more details: https://forum.orekit.org/t/light-time-correction-coupled-with-atmospheric-refraction-issue-in-direct-location/1963
+
return algorithm.refineIntersection(ellipsoid,
- shifted2.transformPosition(pInert),
- shifted2.transformVector(lInert),
- gp1);
+ shifted2.transformPosition(pInert),
+ shifted2.transformVector(lInert),
+ gp1);
}
/**
diff --git a/src/test/java/org/orekit/rugged/TestUtils.java b/src/test/java/org/orekit/rugged/TestUtils.java
index c40a8d7c4e9ec55717d8e8512222edbb1f4df648..fa0f3400a6737def026a538a0be68269b5730747 100644
--- a/src/test/java/org/orekit/rugged/TestUtils.java
+++ b/src/test/java/org/orekit/rugged/TestUtils.java
@@ -217,6 +217,8 @@ public class TestUtils {
// gravity.field.order = 12
AbsoluteDate date = new AbsoluteDate("2012-01-01T00:00:00.000", TimeScalesFactory.getUTC());
Frame eme2000 = FramesFactory.getEME2000();
+
+ // Observation satellite about 800km above ground
return new CircularOrbit(7173352.811913891,
-4.029194321683225E-4, 0.0013530362644647786,
FastMath.toRadians(98.63218182243709),
diff --git a/src/test/java/org/orekit/rugged/refraction/AtmosphericRefractionTest.java b/src/test/java/org/orekit/rugged/refraction/AtmosphericRefractionTest.java
index 5a270b7afa75f59fdd12526f2e432ce34e9a0944..9ebf993c520eceb9db12f01e7081b867afb9641a 100644
--- a/src/test/java/org/orekit/rugged/refraction/AtmosphericRefractionTest.java
+++ b/src/test/java/org/orekit/rugged/refraction/AtmosphericRefractionTest.java
@@ -67,6 +67,7 @@ public class AtmosphericRefractionTest {
int dimension = 4000;
RuggedBuilder builder = initRuggedForAtmosphericTests(dimension, sensorName);
+
// Build Rugged without atmospheric refraction model
Rugged ruggedWithout = builder.build();
@@ -93,7 +94,7 @@ public class AtmosphericRefractionTest {
final double epsilonPixel = pixelThreshold;
final double epsilonLine = lineThreshold;
- double earthRadius = ruggedWithout.getEllipsoid().getEquatorialRadius();
+ final double earthRadius = ruggedWithout.getEllipsoid().getEquatorialRadius();
// Direct loc on a line WITHOUT and WITH atmospheric correction
// ============================================================
@@ -316,6 +317,72 @@ public class AtmosphericRefractionTest {
Assert.assertTrue(pixel != null);
}
+
+ /**
+ * Test for issue #392
+ */
+ @Test
+ public void testLightTimeCorrection() throws URISyntaxException {
+
+ String sensorName = "line";
+ int dimension = 4000;
+
+ RuggedBuilder builder = initRuggedForAtmosphericTests(dimension, sensorName);
+
+ // Build Rugged without atmospheric refraction but with light time correction
+ builder.setLightTimeCorrection(true);
+ Rugged ruggedWithLightTimeWithoutRefraction = builder.build();
+
+ // Defines atmospheric refraction model (with the default multi layers model)
+ AtmosphericRefraction atmosphericRefraction = new MultiLayerModel(builder.getEllipsoid());
+ int pixelStep = 100;
+ int lineStep = 100;
+ atmosphericRefraction.setGridSteps(pixelStep, lineStep);
+
+ // Add atmospheric refraction model
+ builder.setRefractionCorrection(atmosphericRefraction);
+
+ // Build Rugged without light time correction (with atmospheric refraction)
+ builder.setLightTimeCorrection(false);
+ Rugged ruggedWithoutLightTime = builder.build();
+
+ // Build Rugged with light time correction (with atmospheric refraction)
+ builder.setLightTimeCorrection(true);
+ Rugged ruggedWithLightTime = builder.build();
+
+
+ // Compare direct loc on a line :
+ // * with atmospheric refraction, WITHOUT and WITH light time correction:
+ // distance on ground must be not null and < 1.2 m (max shift at equator for orbit at 800km)
+ // * with light time correction, WITHOUT and WITH atmospheric refraction
+ // distance on ground must be not null and < 2 m (max shift due to atmospheric refraction)
+ // =========================================================================================
+ double chosenLine = 200.;
+ GeodeticPoint[] gpWithoutLightTime = ruggedWithoutLightTime.directLocation(sensorName, chosenLine);
+ GeodeticPoint[] gpWithLightTime = ruggedWithLightTime.directLocation(sensorName, chosenLine);
+
+ GeodeticPoint[] gpWithLightTimeWithoutRefraction = ruggedWithLightTimeWithoutRefraction.directLocation(sensorName, chosenLine);
+
+ double earthRadius = builder.getEllipsoid().getEquatorialRadius();
+
+ // Check the shift on the ground
+ for (int i = 0; i < gpWithLightTime.length; i++) {
+
+ double currentRadius = earthRadius + (gpWithLightTime[i].getAltitude() + gpWithoutLightTime[i].getAltitude())/2.;
+ // Compute distance between point (with atmospheric refraction) with light time correction and without
+ double distance = GeodeticUtilities.computeDistanceInMeter(currentRadius, gpWithLightTime[i], gpWithoutLightTime[i]);
+
+ // Check if the distance is not 0 and < 1.2m (at equator max of shift)
+ Assert.assertTrue(distance > 0.0);
+ Assert.assertTrue(distance <= 1.2);
+
+ // Compute distance between point (with light time correction) with refraction and without refraction
+ distance = GeodeticUtilities.computeDistanceInMeter(currentRadius, gpWithLightTime[i], gpWithLightTimeWithoutRefraction[i]);
+ // Check if the distance is not 0 and < 2m
+ Assert.assertTrue(distance > 0.0);
+ Assert.assertTrue(distance < 2.);
+ }
+ }
@Test
public void testBadConfig() {