src/main/java/org/orekit/rugged/intersection/package-info.java

+/* Copyright 2013-2022 CS GROUP
+ * Licensed to CS GROUP (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.
+ */
+/**
+ *
+ * This package provides the interface for Digital Elevation Model intersection
+ * algorithm, as well as some simple implementations.
+ *
+ *
+ * @author Luc Maisonobe
+ * @author Guylaine Prat
+ *
+ */
+package org.orekit.rugged.intersection;

src/main/java/org/orekit/rugged/api/LineDatation.java → src/main/java/org/orekit/rugged/linesensor/LineDatation.java

-/* Copyright 2013-2014 CS Systèmes d'Information
- * Licensed to CS Systèmes d'Information (CS) under one or more
+/* Copyright 2013-2022 CS GROUP
+ * Licensed to CS GROUP (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
@@ -14,7 +14,9 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
-package org.orekit.rugged.api;
+package org.orekit.rugged.linesensor;
+
+import org.orekit.time.AbsoluteDate;
 
 /** Interface representing line datation model.
  * @see LinearLineDatation
@@ -24,14 +26,20 @@ public interface LineDatation {
 
     /** Get the date for a given line.
      * @param lineNumber line number
-     * @return date, as an offset in seconds from reference date
+     * @return date at which line is acquired
      */
-    double getDate(double lineNumber);
+    AbsoluteDate getDate(double lineNumber);
 
     /** Get the line for a given date.
-     * @param date date, as an offset in seconds from reference date
+     * @param date date
      * @return line number
      */
-    double getLine(double date);
+    double getLine(AbsoluteDate date);
+
+    /** Get the rate of lines scanning.
+     * @param lineNumber line number
+     * @return rate of lines scanning (lines / seconds)
+     */
+    double getRate(double lineNumber);
 
 }

src/main/java/org/orekit/rugged/linesensor/LineSensor.java

+/* Copyright 2013-2022 CS GROUP
+ * Licensed to CS GROUP (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.linesensor;
+
+import java.util.stream.Stream;
+
+import org.hipparchus.analysis.differentiation.Derivative;
+import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
+import org.hipparchus.geometry.euclidean.threed.Vector3D;
+import org.hipparchus.util.FastMath;
+import org.orekit.rugged.errors.DumpManager;
+import org.orekit.rugged.los.TimeDependentLOS;
+import org.orekit.rugged.utils.DerivativeGenerator;
+import org.orekit.time.AbsoluteDate;
+import org.orekit.utils.ParameterDriver;
+
+/** Line sensor model.
+ * @author Luc Maisonobe
+ * @author Guylaine Prat
+ */
+public class LineSensor {
+
+    /** Name of the sensor. */
+    private final String name;
+
+    /** Datation model. */
+    private final LineDatation datationModel;
+
+    /** Sensor position. */
+    private final Vector3D position;
+
+    /** Pixels lines-of-sight. */
+    private final TimeDependentLOS los;
+
+    /** Simple constructor.
+     * @param name name of the sensor
+     * @param datationModel datation model
+     * @param position sensor position in spacecraft frame
+     * @param los pixels lines-of-sight in spacecraft frame
+     * @see org.orekit.rugged.los.LOSBuilder
+     */
+    public LineSensor(final String name, final LineDatation datationModel,
+                      final Vector3D position, final TimeDependentLOS los) {
+
+        this.name          = name;
+        this.datationModel = datationModel;
+        this.position      = position;
+        this.los           = los;
+
+    }
+
+    /** Get the name of the sensor.
+     * @return name of the sensor
+     */
+    public String getName() {
+        return name;
+    }
+
+    /** Get the number of pixels.
+     * @return number of pixels
+     */
+    public int getNbPixels() {
+        return los.getNbPixels();
+    }
+
+    /** Get the drivers for LOS parameters.
+     * @return drivers for LOS parameters
+     * @since 2.0
+     */
+    public Stream<ParameterDriver> getParametersDrivers() {
+        return los.getParametersDrivers();
+    }
+
+    /** Get the pixel normalized line-of-sight at some date.
+     * @param date current date
+     * @param i pixel index (must be between 0 and {@link #getNbPixels()} - 1
+     * @return pixel normalized line-of-sight
+     */
+    public Vector3D getLOS(final AbsoluteDate date, final int i) {
+        final Vector3D l = los.getLOS(i, date);
+        DumpManager.dumpSensorLOS(this, date, i, l);
+        return l;
+    }
+
+    /** Get the pixel normalized interpolated line-of-sight at some date.
+     * @param date current date
+     * @param i pixel index (must be between 0 and {@link #getNbPixels()} - 1
+     * @return pixel normalized line-of-sight
+     * @since 2.0
+     */
+    public Vector3D getLOS(final AbsoluteDate date, final double i) {
+
+        final int iInf = FastMath.max(0, FastMath.min(getNbPixels() - 2, (int) FastMath.floor(i)));
+        final int iSup = iInf + 1;
+        final Vector3D interpolatedLos     = new Vector3D(iSup - i, los.getLOS(iInf, date),
+                                                  i - iInf, los.getLOS(iSup, date));
+        return interpolatedLos;
+    }
+
+    /** Get the pixel normalized line-of-sight at some date,
+     * and their derivatives with respect to estimated parameters.
+     * @param <T> derivative type
+     * @param date current date
+     * @param i pixel index (must be between 0 and {@link #getNbPixels()} - 1
+     * @param generator generator to use for building {@link Derivative} instances
+     * @return pixel normalized line-of-sight
+     */
+    public <T extends Derivative<T>> FieldVector3D<T> getLOSDerivatives(final AbsoluteDate date, final int i,
+                                                                        final DerivativeGenerator<T> generator) {
+        return los.getLOSDerivatives(i, date, generator);
+    }
+
+    /** Get the pixel normalized line-of-sight at some date,
+     * and their derivatives with respect to estimated parameters.
+     * @param <T> derivative type
+     * @param date current date
+     * @param i pixel index (must be between 0 and {@link #getNbPixels()} - 1
+     * @param generator generator to use for building {@link Derivative} instances
+     * @return pixel normalized line-of-sight
+     * @since 2.0
+     */
+    public <T extends Derivative<T>> FieldVector3D<T> getLOSDerivatives(final AbsoluteDate date, final double i,
+                                                                        final DerivativeGenerator<T> generator) {
+
+        // find surrounding pixels of pixelB (in order to interpolate LOS from pixelB (that is not an integer)
+        final int iInf = FastMath.max(0, FastMath.min(getNbPixels() - 2, (int) FastMath.floor(i)));
+        final int iSup = iInf + 1;
+
+        final FieldVector3D<T> interpolatedLos =
+                        new FieldVector3D<> (iSup - i, los.getLOSDerivatives(iInf, date, generator),
+                                             i - iInf, los.getLOSDerivatives(iSup, date, generator)).normalize();
+        return interpolatedLos;
+    }
+
+    /** Get the date.
+     * @param lineNumber line number
+     * @return date corresponding to line number
+     */
+    public AbsoluteDate getDate(final double lineNumber) {
+        final AbsoluteDate date = datationModel.getDate(lineNumber);
+        DumpManager.dumpSensorDatation(this, lineNumber, date);
+        return date;
+    }
+
+    /** Get the line number.
+     * @param date date
+     * @return line number corresponding to date
+     */
+    public double getLine(final AbsoluteDate date) {
+        final double lineNumber = datationModel.getLine(date);
+        DumpManager.dumpSensorDatation(this, lineNumber, date);
+        return lineNumber;
+    }
+
+    /** Get the rate of lines scanning.
+     * @param lineNumber line number
+     * @return rate of lines scanning (lines / seconds)
+     */
+    public double getRate(final double lineNumber) {
+        final double rate = datationModel.getRate(lineNumber);
+        DumpManager.dumpSensorRate(this, lineNumber, rate);
+        return rate;
+    }
+
+    /** Get the sensor position.
+     * @return position
+     */
+    public Vector3D getPosition() {
+        return position;
+    }
+
+    /** Dump the rate for the current line number.
+     * @param lineNumber line number
+     */
+    public void dumpRate(final double lineNumber) {
+        final double rate = datationModel.getRate(lineNumber);
+        DumpManager.dumpSensorRate(this, lineNumber, rate);
+    }
+}

src/main/java/org/orekit/rugged/api/LinearLineDatation.java → src/main/java/org/orekit/rugged/linesensor/LinearLineDatation.java

-/* Copyright 2013-2014 CS Systèmes d'Information
- * Licensed to CS Systèmes d'Information (CS) under one or more
+/* Copyright 2013-2022 CS GROUP
+ * Licensed to CS GROUP (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
@@ -14,7 +14,9 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
-package org.orekit.rugged.api;
+package org.orekit.rugged.linesensor;
+
+import org.orekit.time.AbsoluteDate;
 
 
 /** Linear model for {@link LineDatation line datation}.
@@ -25,31 +27,43 @@ package org.orekit.rugged.api;
  */
 public class LinearLineDatation implements LineDatation {
 
+    /** Reference date. */
+    private final AbsoluteDate referenceDate;
+
     /** Line number at reference date. */
-    private final double line0;
+    private final double referenceLine;
 
     /** Rate of lines scanning (lines / seconds). */
     private final double rate;
 
     /** Simple constructor.
-     * @param line0 line number at reference date
+     * @param referenceDate reference date
+     * @param referenceLine line number at reference date
      * @param rate rate of lines scanning (lines / seconds)
      */
-    public LinearLineDatation(final double line0, final double rate) {
-        this.line0 = line0;
-        this.rate  = rate;
+    public LinearLineDatation(final AbsoluteDate referenceDate, final double referenceLine,
+                              final double rate) {
+        this.referenceDate = referenceDate;
+        this.referenceLine = referenceLine;
+        this.rate          = rate;
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public AbsoluteDate getDate(final double lineNumber) {
+        return referenceDate.shiftedBy((lineNumber - referenceLine) / rate);
     }
 
     /** {@inheritDoc} */
     @Override
-    public double getDate(final double lineNumber) {
-        return (lineNumber - line0) / rate;
+    public double getLine(final AbsoluteDate date) {
+        return referenceLine + rate * date.durationFrom(referenceDate);
     }
 
     /** {@inheritDoc} */
     @Override
-    public double getLine(final double date) {
-        return line0 + rate * date;
+    public double getRate(final double lineNumber) {
+        return rate;
     }
 
 }

src/main/java/org/orekit/rugged/linesensor/SensorMeanPlaneCrossing.java

+/* Copyright 2013-2022 CS GROUP
+ * Licensed to CS GROUP (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.linesensor;
+
+import java.util.ArrayList;
+import java.util.List;
+import java.util.stream.Stream;
+
+import org.hipparchus.analysis.UnivariateFunction;
+import org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver;
+import org.hipparchus.analysis.solvers.UnivariateSolver;
+import org.hipparchus.exception.MathIllegalArgumentException;
+import org.hipparchus.exception.MathIllegalStateException;
+import org.hipparchus.geometry.euclidean.threed.Vector3D;
+import org.hipparchus.linear.Array2DRowRealMatrix;
+import org.hipparchus.linear.ArrayRealVector;
+import org.hipparchus.linear.DecompositionSolver;
+import org.hipparchus.linear.MatrixUtils;
+import org.hipparchus.linear.QRDecomposition;
+import org.hipparchus.linear.RealMatrix;
+import org.hipparchus.linear.RealVector;
+import org.hipparchus.linear.SingularValueDecomposition;
+import org.hipparchus.util.FastMath;
+import org.hipparchus.util.Precision;
+import org.orekit.frames.Transform;
+import org.orekit.rugged.errors.RuggedException;
+import org.orekit.rugged.errors.RuggedInternalError;
+import org.orekit.rugged.utils.SpacecraftToObservedBody;
+import org.orekit.time.AbsoluteDate;
+import org.orekit.utils.Constants;
+import org.orekit.utils.PVCoordinates;
+
+/** Class dedicated to find when ground point crosses mean sensor plane.
+ * <p>
+ * This class is used in the first stage of inverse location.
+ * </p>
+ * @author Luc Maisonobe
+ */
+public class SensorMeanPlaneCrossing {
+
+    /** Number of cached results for guessing the line number. */
+    private static final int CACHED_RESULTS = 6;
+
+    /** Converter between spacecraft and body. */
+    private final SpacecraftToObservedBody scToBody;
+
+    /** Middle line. */
+    private final double midLine;
+
+    /** Body to inertial transform for middle line. */
+    private final Transform midBodyToInert;
+
+    /** Spacecraft to inertial transform for middle line. */
+    private final Transform midScToInert;
+
+    /** Minimum line number in the search interval. */
+    private final int minLine;
+
+    /** Maximum line number in the search interval. */
+    private final int maxLine;
+
+    /** Flag for light time correction. */
+    private final boolean lightTimeCorrection;
+
+    /** Flag for aberration of light correction. */
+    private final boolean aberrationOfLightCorrection;
+
+    /** Line sensor. */
+    private final LineSensor sensor;
+
+    /** Sensor mean plane normal. */
+    private final Vector3D meanPlaneNormal;
+
+    /** Maximum number of evaluations. */
+    private final int maxEval;
+
+    /** Accuracy to use for finding crossing line number. */
+    private final double accuracy;
+
+    /** Cached results. */
+    private final List<CrossingResult> cachedResults;
+
+    /** Simple constructor.
+     * @param sensor sensor to consider
+     * @param scToBody converter between spacecraft and body
+     * @param minLine minimum line number
+     * @param maxLine maximum line number
+     * @param lightTimeCorrection flag for light time correction
+     * @param aberrationOfLightCorrection flag for aberration of light correction.
+     * @param maxEval maximum number of evaluations
+     * @param accuracy accuracy to use for finding crossing line number
+     */
+    public SensorMeanPlaneCrossing(final LineSensor sensor,
+                                   final SpacecraftToObservedBody scToBody,
+                                   final int minLine, final int maxLine,
+                                   final boolean lightTimeCorrection,
+                                   final boolean aberrationOfLightCorrection,
+                                   final int maxEval, final double accuracy) {
+        this(sensor, scToBody, minLine, maxLine, lightTimeCorrection, aberrationOfLightCorrection,
+             maxEval, accuracy, computeMeanPlaneNormal(sensor, minLine, maxLine),
+             Stream.<CrossingResult>empty());
+    }
+
+    /** Simple constructor.
+     * @param sensor sensor to consider
+     * @param scToBody converter between spacecraft and body
+     * @param minLine minimum line number
+     * @param maxLine maximum line number
+     * @param lightTimeCorrection flag for light time correction
+     * @param aberrationOfLightCorrection flag for aberration of light correction.
+     * @param maxEval maximum number of evaluations
+     * @param accuracy accuracy to use for finding crossing line number
+     * @param meanPlaneNormal mean plane normal
+     * @param cachedResults cached results
+     */
+    public SensorMeanPlaneCrossing(final LineSensor sensor,
+                                   final SpacecraftToObservedBody scToBody,
+                                   final int minLine, final int maxLine,
+                                   final boolean lightTimeCorrection,
+                                   final boolean aberrationOfLightCorrection,
+                                   final int maxEval, final double accuracy,
+                                   final Vector3D meanPlaneNormal,
+                                   final Stream<CrossingResult> cachedResults) {
+
+        this.sensor                      = sensor;
+        this.minLine                     = minLine;
+        this.maxLine                     = maxLine;
+        this.lightTimeCorrection         = lightTimeCorrection;
+        this.aberrationOfLightCorrection = aberrationOfLightCorrection;
+        this.maxEval                     = maxEval;
+        this.accuracy                    = accuracy;
+        this.scToBody                    = scToBody;
+
+        this.midLine                     = 0.5 * (minLine + maxLine);
+        final AbsoluteDate midDate       = sensor.getDate(midLine);
+        this.midBodyToInert              = scToBody.getBodyToInertial(midDate);
+        this.midScToInert                = scToBody.getScToInertial(midDate);
+
+        this.meanPlaneNormal             = meanPlaneNormal;
+
+        this.cachedResults               = new ArrayList<>(CACHED_RESULTS);
+        cachedResults.forEach(crossingResult -> {
+            if (crossingResult != null && this.cachedResults.size() < CACHED_RESULTS) {
+                this.cachedResults.add(crossingResult);
+            }
+        });
+
+    }
+
+    /** Compute the plane containing origin that best fits viewing directions point cloud.
+     * @param sensor line sensor
+     * @param minLine minimum line number
+     * @param maxLine maximum line number
+     * <p>
+     * The normal is oriented such that traversing pixels in increasing indices
+     * order corresponds to trigonometric order (i.e. counterclockwise).
+     * </p>
+     * @return normal of the mean plane
+     */
+    private static Vector3D computeMeanPlaneNormal(final LineSensor sensor, final int minLine, final int maxLine) {
+
+        final AbsoluteDate midDate = sensor.getDate(0.5 * (minLine + maxLine));
+
+        // build a centered data matrix
+        // (for each viewing direction, we add both the direction and its
+        //  opposite, thus ensuring the plane will contain origin)
+        final RealMatrix matrix = MatrixUtils.createRealMatrix(3, 2 * sensor.getNbPixels());
+        for (int i = 0; i < sensor.getNbPixels(); ++i) {
+            final Vector3D l = sensor.getLOS(midDate, i);
+            matrix.setEntry(0, 2 * i,      l.getX());
+            matrix.setEntry(1, 2 * i,      l.getY());
+            matrix.setEntry(2, 2 * i,      l.getZ());
+            matrix.setEntry(0, 2 * i + 1, -l.getX());
+            matrix.setEntry(1, 2 * i + 1, -l.getY());
+            matrix.setEntry(2, 2 * i + 1, -l.getZ());
+        }
+
+        // compute Singular Value Decomposition
+        final SingularValueDecomposition svd = new SingularValueDecomposition(matrix);
+
+        // extract the left singular vector corresponding to least singular value
+        // (i.e. last vector since Hipparchus returns the values
+        //  in non-increasing order)
+        final Vector3D singularVector = new Vector3D(svd.getU().getColumn(2)).normalize();
+
+        // check rotation order
+        final Vector3D first = sensor.getLOS(midDate, 0);
+        final Vector3D last  = sensor.getLOS(midDate, sensor.getNbPixels() - 1);
+        if (Vector3D.dotProduct(singularVector, Vector3D.crossProduct(first, last)) >= 0) {
+            return singularVector;
+        } else {
+            return singularVector.negate();
+        }
+
+    }
+
+    /** Get the underlying sensor.
+     * @return underlying sensor
+     */
+    public LineSensor getSensor() {
+        return sensor;
+    }
+
+    /** Get converter between spacecraft and body.
+     * @return converter between spacecraft and body
+     */
+    public SpacecraftToObservedBody getScToBody() {
+        return scToBody;
+    }
+
+    /** Get the minimum line number in the search interval.
+     * @return minimum line number in the search interval
+     */
+    public int getMinLine() {
+        return minLine;
+    }
+
+    /** Get the maximum line number in the search interval.
+     * @return maximum line number in the search interval
+     */
+    public int getMaxLine() {
+        return maxLine;
+    }
+
+    /** Get the maximum number of evaluations.
+     * @return maximum number of evaluations
+     */
+    public int getMaxEval() {
+        return maxEval;
+    }
+
+    /** Get the accuracy to use for finding crossing line number.
+     * @return accuracy to use for finding crossing line number
+     */
+    public double getAccuracy() {
+        return accuracy;
+    }
+
+    /** Get the mean plane normal.
+     * <p>
+     * The normal is oriented such traversing pixels in increasing indices
+     * order corresponds is consistent with trigonometric order (i.e.
+     * counterclockwise).
+     * </p>
+     * @return mean plane normal
+     */
+    public Vector3D getMeanPlaneNormal() {
+        return meanPlaneNormal;
+    }
+
+    /** Get cached previous results.
+     * @return cached previous results
+     */
+    public Stream<CrossingResult> getCachedResults() {
+        return cachedResults.stream();
+    }
+
+    /** Container for mean plane crossing result. */
+    public static class CrossingResult {
+
+        /** Crossing date. */
+        private final AbsoluteDate crossingDate;
+
+        /** Crossing line. */
+        private final double crossingLine;
+
+        /** Target ground point. */
+        private final Vector3D target;
+
+        /** Target direction in spacecraft frame. */
+        private final Vector3D targetDirection;
+
+        /** Derivative of the target direction in spacecraft frame with respect to line number. */
+        private final Vector3D targetDirectionDerivative;
+
+        /** Simple constructor.
+         * @param crossingDate crossing date
+         * @param crossingLine crossing line
+         * @param target target ground point
+         * @param targetDirection target direction in spacecraft frame
+         * @param targetDirectionDerivative derivative of the target direction in spacecraft frame with respect to line number
+         */
+        public CrossingResult(final AbsoluteDate crossingDate, final double crossingLine,
+                              final Vector3D target,
+                              final Vector3D targetDirection,
+                              final Vector3D targetDirectionDerivative) {
+            this.crossingDate              = crossingDate;
+            this.crossingLine              = crossingLine;
+            this.target                    = target;
+            this.targetDirection           = targetDirection;
+            this.targetDirectionDerivative = targetDirectionDerivative;
+        }
+
+        /** Get the crossing date.
+         * @return crossing date
+         */
+        public AbsoluteDate getDate() {
+            return crossingDate;
+        }
+
+        /** Get the crossing line.
+         * @return crossing line
+         */
+        public double getLine() {
+            return crossingLine;
+        }
+
+        /** Get the target ground point.
+         * @return target ground point
+         */
+        public Vector3D getTarget() {
+            return target;
+        }
+
+        /** Get the normalized target direction in spacecraft frame at crossing.
+         * @return normalized target direction in spacecraft frame at crossing
+         * with respect to line number
+         */
+        public Vector3D getTargetDirection() {
+            return targetDirection;
+        }
+
+        /** Get the derivative of the normalized target direction in spacecraft frame at crossing.
+         * @return derivative of the normalized target direction in spacecraft frame at crossing
+         * with respect to line number
+         * @since 2.0
+         */
+        public Vector3D getTargetDirectionDerivative() {
+            return targetDirectionDerivative;
+        }
+
+    }
+
+    /** Find mean plane crossing.
+     * @param target target ground point
+     * @return line number and target direction at mean plane crossing,
+     * or null if search interval does not bracket a solution
+     */
+    public CrossingResult find(final Vector3D target) {
+
+        double crossingLine     = midLine;
+        Transform bodyToInert   = midBodyToInert;
+        Transform scToInert     = midScToInert;
+
+        // count the number of available results
+        if (cachedResults.size() >= 4) {
+            // we already have computed at lest 4 values, we attempt to build a linear
+            // model to guess a better start line
+            final double guessedCrossingLine = guessStartLine(target);
+            if (guessedCrossingLine >= minLine && guessedCrossingLine <= maxLine) {
+                crossingLine = guessedCrossingLine;
+                final AbsoluteDate date = sensor.getDate(crossingLine);
+                bodyToInert = scToBody.getBodyToInertial(date);
+                scToInert   = scToBody.getScToInertial(date);
+            }
+        }
+
+        final PVCoordinates targetPV = new PVCoordinates(target, Vector3D.ZERO);
+
+        // we don't use an Hipparchus solver here because we are more
+        // interested in reducing the number of evaluations than being accurate,
+        // as we know the solution is improved in the second stage of inverse location.
+        // We expect two or three evaluations only. Each new evaluation shows up quickly in
+        // the performances as it involves frames conversions
+        final double[]  crossingLineHistory = new double[maxEval];
+        final double[]  betaHistory         = new double[maxEval];
+        final double[]  betaDerHistory      = new double[maxEval];
+        boolean         atMin               = false;
+        boolean         atMax               = false;
+        for (int i = 0; i < maxEval; ++i) {
+
+            crossingLineHistory[i] = crossingLine;
+            final Vector3D[] targetDirection =
+                    evaluateLine(crossingLine, targetPV, bodyToInert, scToInert);
+            betaHistory[i] = FastMath.acos(Vector3D.dotProduct(targetDirection[0], meanPlaneNormal));
+            final double s = Vector3D.dotProduct(targetDirection[1], meanPlaneNormal);
+            betaDerHistory[i] = -s / FastMath.sqrt(1 - s * s);
+
+            final double deltaL;
+            if (i == 0) {
+                // simple Newton iteration
+                deltaL        = (0.5 * FastMath.PI - betaHistory[i]) / betaDerHistory[i];
+                crossingLine += deltaL;
+            } else {
+                // inverse cubic iteration
+                final double a0    = betaHistory[i - 1] - 0.5 * FastMath.PI;
+                final double l0    = crossingLineHistory[i - 1];
+                final double d0    = betaDerHistory[i - 1];
+                final double a1    = betaHistory[i]     - 0.5 * FastMath.PI;
+                final double l1    = crossingLineHistory[i];
+                final double d1    = betaDerHistory[i];
+                final double a1Ma0 = a1 - a0;
+                crossingLine = ((l0 * (a1 - 3 * a0) - a0 * a1Ma0 / d0) * a1 * a1 +
+                                (l1 * (3 * a1 - a0) - a1 * a1Ma0 / d1) * a0 * a0
+                               ) / (a1Ma0 * a1Ma0 * a1Ma0);
+                deltaL = crossingLine - l1;
+            }
+            if (FastMath.abs(deltaL) <= accuracy) {
+                // return immediately, without doing any additional evaluation!
+                final CrossingResult crossingResult =
+                    new CrossingResult(sensor.getDate(crossingLine), crossingLine, target,
+                                       targetDirection[0], targetDirection[1]);
+                boolean isNew = true;
+                for (final CrossingResult existing : cachedResults) {
+                    isNew = isNew && FastMath.abs(crossingLine - existing.crossingLine) > accuracy;
+                }
+                if (isNew) {
+                    // this result is different from the existing ones,
+                    // it brings new sampling data to the cache
+                    if (cachedResults.size() >= CACHED_RESULTS) {
+                        cachedResults.remove(cachedResults.size() - 1);
+                    }
+                    cachedResults.add(0, crossingResult);
+                }
+                return crossingResult;
+            }
+            for (int j = 0; j < i; ++j) {
+                if (FastMath.abs(crossingLine - crossingLineHistory[j]) <= 1.0) {
+                    // rare case: we are stuck in a loop!
+                    // switch to a more robust (but slower) algorithm in this case
+                    final CrossingResult slowResult = slowFind(targetPV, crossingLine);
+                    if (slowResult == null) {
+                        return null;
+                    }
+                    if (cachedResults.size() >= CACHED_RESULTS) {
+                        cachedResults.remove(cachedResults.size() - 1);
+                    }
+                    cachedResults.add(0, slowResult);
+                    return cachedResults.get(0);
+                }
+            }
+
+            if (crossingLine < minLine) {
+                if (atMin) {
+                    // we were already trying at minLine and we need to go below that
+                    // give up as the solution is out of search interval
+                    return null;
+                }
+                atMin        = true;
+                crossingLine = minLine;
+            } else if (crossingLine > maxLine) {
+                if (atMax) {
+                    // we were already trying at maxLine and we need to go above that
+                    // give up as the solution is out of search interval
+                    return null;
+                }
+                atMax        = true;
+                crossingLine = maxLine;
+            } else {
+                // the next evaluation will be a regular point
+                atMin = false;
+                atMax = false;
+            }
+
+            final AbsoluteDate date = sensor.getDate(crossingLine);
+            bodyToInert = scToBody.getBodyToInertial(date);
+            scToInert   = scToBody.getScToInertial(date);
+        }
+
+        return null;
+
+    }
+
+    /** Guess a start line using the last four results.
+     * @param target target ground point
+     * @return guessed start line
+     */
+    private double guessStartLine(final Vector3D target) {
+        try {
+
+            // assume a linear model of the form: l = ax + by + cz + d
+            final int        n = cachedResults.size();
+            final RealMatrix m = new Array2DRowRealMatrix(n, 4);
+            final RealVector v = new ArrayRealVector(n);
+            for (int i = 0; i < n; ++i) {
+                final CrossingResult crossingResult = cachedResults.get(i);
+                m.setEntry(i, 0, crossingResult.getTarget().getX());
+                m.setEntry(i, 1, crossingResult.getTarget().getY());
+                m.setEntry(i, 2, crossingResult.getTarget().getZ());
+                m.setEntry(i, 3, 1.0);
+                v.setEntry(i, crossingResult.getLine());
+            }
+
+            final DecompositionSolver solver = new QRDecomposition(m, Precision.SAFE_MIN).getSolver();
+            final RealVector coefficients = solver.solve(v);
+
+            // apply the linear model
+            return target.getX() * coefficients.getEntry(0) +
+                   target.getY() * coefficients.getEntry(1) +
+                   target.getZ() * coefficients.getEntry(2) +
+                   coefficients.getEntry(3);
+
+        } catch (MathIllegalStateException sme) {
+            // the points are not independent
+            return Double.NaN;
+        }
+    }
+
+    /** Find mean plane crossing using a slow but robust method.
+     * @param targetPV target ground point
+     * @param initialGuess initial guess for the crossing line
+     * @return line number and target direction at mean plane crossing,
+     * or null if search interval does not bracket a solution
+     */
+    private CrossingResult slowFind(final PVCoordinates targetPV, final double initialGuess) {
+
+        try {
+            // safety check
+            final double startValue;
+            if (initialGuess < minLine || initialGuess > maxLine) {
+                startValue = midLine;
+            } else {
+                startValue = initialGuess;
+            }
+
+            final UnivariateSolver solver = new BracketingNthOrderBrentSolver(accuracy, 5);
+            final double crossingLine = solver.solve(maxEval, new UnivariateFunction() {
+                /** {@inheritDoc} */
+                @Override
+                public double value(final double x) {
+                    try {
+                        final AbsoluteDate date = sensor.getDate(x);
+                        final Vector3D[] targetDirection =
+                                evaluateLine(x, targetPV, scToBody.getBodyToInertial(date), scToBody.getScToInertial(date));
+                        return 0.5 * FastMath.PI - FastMath.acos(Vector3D.dotProduct(targetDirection[0], meanPlaneNormal));
+                    } catch (RuggedException re) {
+                        throw new RuggedInternalError(re);
+                    }
+                }
+            }, minLine, maxLine, startValue);
+
+            final AbsoluteDate date = sensor.getDate(crossingLine);
+            final Vector3D[] targetDirection =
+                    evaluateLine(crossingLine, targetPV, scToBody.getBodyToInertial(date), scToBody.getScToInertial(date));
+            return new CrossingResult(sensor.getDate(crossingLine), crossingLine, targetPV.getPosition(),
+                                      targetDirection[0], targetDirection[1]);
+
+        } catch (MathIllegalArgumentException nbe) {
+            return null;
+        }
+    }
+
+    /** Evaluate geometry for a given line number.
+     * @param lineNumber current line number
+     * @param targetPV target ground point
+     * @param bodyToInert transform from observed body to inertial frame, for current line
+     * @param scToInert transform from inertial frame to spacecraft frame, for current line
+     * @return target direction in spacecraft frame, with its first derivative
+     * with respect to line number
+     */
+    private Vector3D[] evaluateLine(final double lineNumber, final PVCoordinates targetPV,
+                                    final Transform bodyToInert, final Transform scToInert) {
+
+        // compute the transform between spacecraft and observed body
+        final PVCoordinates refInert =
+                scToInert.transformPVCoordinates(new PVCoordinates(sensor.getPosition(), Vector3D.ZERO));
+
+        final PVCoordinates targetInert;
+        if (lightTimeCorrection) {
+            // apply light time correction
+            final Vector3D iT     = bodyToInert.transformPosition(targetPV.getPosition());
+            final double   deltaT = refInert.getPosition().distance(iT) / Constants.SPEED_OF_LIGHT;
+            targetInert           = bodyToInert.shiftedBy(-deltaT).transformPVCoordinates(targetPV);
+        } else {
+            // don't apply light time correction
+            targetInert = bodyToInert.transformPVCoordinates(targetPV);
+        }
+
+        final PVCoordinates lInert    = new PVCoordinates(refInert, targetInert);
+        final Transform     inertToSc = scToInert.getInverse();
+        final Vector3D obsLInert;
+        final Vector3D obsLInertDot;
+        if (aberrationOfLightCorrection) {
+
+            // apply aberration of light correction
+            // as the spacecraft velocity is small with respect to speed of light,
+            // we use classical velocity addition and not relativistic velocity addition
+            // we have: c * lInert + vsat = k * obsLInert
+            final PVCoordinates spacecraftPV = scToInert.transformPVCoordinates(PVCoordinates.ZERO);
+            final Vector3D  l        = lInert.getPosition().normalize();
+            final Vector3D  lDot     = normalizedDot(lInert.getPosition(), lInert.getVelocity());
+            final Vector3D  kObs     = new Vector3D(Constants.SPEED_OF_LIGHT, l, +1.0, spacecraftPV.getVelocity());
+            obsLInert = kObs.normalize();
+
+            // the following derivative is computed under the assumption the spacecraft velocity
+            // is constant in inertial frame ... It is obviously not true, but as this velocity
+            // is very small with respect to speed of light, the error is expected to remain small
+            obsLInertDot = normalizedDot(kObs, new Vector3D(Constants.SPEED_OF_LIGHT, lDot));
+
+        } else {
+
+            // don't apply aberration of light correction
+            obsLInert    = lInert.getPosition().normalize();
+            obsLInertDot = normalizedDot(lInert.getPosition(), lInert.getVelocity());
+
+        }
+        final Vector3D dir    = inertToSc.transformVector(obsLInert);
+        final Vector3D dirDot = new Vector3D(+1.0, inertToSc.transformVector(obsLInertDot),
+                                             -1.0, Vector3D.crossProduct(inertToSc.getRotationRate(), dir));
+
+        // combine vector value and derivative
+        final double rate = sensor.getRate(lineNumber);
+        return new Vector3D[] {
+            dir, dirDot.scalarMultiply(1.0 / rate)
+        };
+
+    }
+
+    /** Compute the derivative of normalized vector.
+     * @param u base vector
+     * @param uDot derivative of the base vector
+     * @return vDot, where v = u / ||u||
+     */
+    private Vector3D normalizedDot(final Vector3D u, final Vector3D uDot) {
+        final double n = u.getNorm();
+        return new Vector3D(1.0 / n, uDot, -Vector3D.dotProduct(u, uDot) / (n * n * n), u);
+    }
+
+}

src/main/java/org/orekit/rugged/api/SensorPixel.java → src/main/java/org/orekit/rugged/linesensor/SensorPixel.java

-/* Copyright 2013-2014 CS Systèmes d'Information
- * Licensed to CS Systèmes d'Information (CS) under one or more
+/* Copyright 2013-2022 CS GROUP
+ * Licensed to CS GROUP (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
@@ -14,7 +14,7 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
-package org.orekit.rugged.api;
+package org.orekit.rugged.linesensor;
 
 import java.io.Serializable;
 

src/main/java/org/orekit/rugged/linesensor/SensorPixelCrossing.java

+/* Copyright 2013-2022 CS GROUP
+ * Licensed to CS GROUP (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.linesensor;
+
+import org.hipparchus.analysis.UnivariateFunction;
+import org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver;
+import org.hipparchus.analysis.solvers.UnivariateSolver;
+import org.hipparchus.exception.MathIllegalArgumentException;
+import org.hipparchus.geometry.euclidean.threed.Vector3D;
+import org.hipparchus.util.FastMath;
+import org.orekit.rugged.errors.RuggedException;
+import org.orekit.rugged.errors.RuggedInternalError;
+import org.orekit.time.AbsoluteDate;
+
+/** Class devoted to locate where ground point crosses a sensor line.
+ * <p>
+ * This class is used in the first stage of inverse location.
+ * </p>
+ * @author Luc Maisonobe
+ */
+public class SensorPixelCrossing {
+
+    /** Margin before and after end pixels, in order to avoid search failures near boundaries. */
+    private static final double MARGIN = 10.0;
+
+    /** Line sensor. */
+    private final LineSensor sensor;
+
+    /** Cross direction in spacecraft frame. */
+    private final Vector3D cross;
+
+    /** Maximum number of evaluations. */
+    private final int maxEval;
+
+    /** Accuracy to use for finding crossing line number. */
+    private final double accuracy;
+
+    /** Simple constructor.
+     * @param sensor sensor to consider
+     * @param meanNormal mean plane normal of the line sensor
+     * @param targetDirection target direction in spacecraft frame
+     * @param maxEval maximum number of evaluations
+     * @param accuracy accuracy to use for finding crossing line number
+     */
+    public SensorPixelCrossing(final LineSensor sensor, final Vector3D meanNormal,
+                               final Vector3D targetDirection,
+                               final int maxEval, final double accuracy) {
+        this.sensor   = sensor;
+        this.cross    = Vector3D.crossProduct(meanNormal, targetDirection).normalize();
+        this.maxEval  = maxEval;
+        this.accuracy = accuracy;
+    }
+
+    /** Locate pixel along sensor line.
+     * @param date current date
+     * @return pixel location ({@code Double.NaN} if the first and last
+     * pixels of the line do not bracket a location)
+     */
+    public double locatePixel(final AbsoluteDate date) {
+        try {
+
+            // set up function evaluating to 0.0 where target matches pixel
+            final UnivariateFunction f = new UnivariateFunction() {
+                /** {@inheritDoc} */
+                @Override
+                public double value(final double x) {
+                    try {
+                        return Vector3D.angle(cross, getLOS(date, x)) - 0.5 * FastMath.PI;
+                    } catch (RuggedException re) {
+                        throw new RuggedInternalError(re);
+                    }
+                }
+            };
+
+            // find the root
+            final UnivariateSolver solver =
+                    new BracketingNthOrderBrentSolver(0.0, accuracy, 5);
+            return solver.solve(maxEval, f, -MARGIN, sensor.getNbPixels() - 1 + MARGIN);
+
+        } catch (MathIllegalArgumentException nbe) {
+            // there are no solutions in the search interval
+            return Double.NaN;
+        }
+    }
+
+    /** Interpolate sensor pixels at some pixel index.
+     * @param date current date
+     * @param x pixel index
+     * @return interpolated direction for specified index
+     */
+    private Vector3D getLOS(final AbsoluteDate date, final double x) {
+
+        // find surrounding pixels
+        final int iInf = FastMath.max(0, FastMath.min(sensor.getNbPixels() - 2, (int) FastMath.floor(x)));
+        final int iSup = iInf + 1;
+
+        // interpolate
+        return new Vector3D(iSup - x, sensor.getLOS(date, iInf),
+                            x - iInf, sensor.getLOS(date, iSup)).normalize();
+
+    }
+
+}

src/main/java/org/orekit/rugged/linesensor/package-info.java

+/* Copyright 2013-2022 CS GROUP
+ * Licensed to CS GROUP (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.
+ */
+/**
+ *
+ * This package provides all the tools needed to deal with line sensor model.
+ *
+ * @author Luc Maisonobe
+ * @author Guylaine Prat
+ *
+ */
+package org.orekit.rugged.linesensor;

src/main/java/org/orekit/rugged/los/FixedRotation.java

+/* Copyright 2013-2022 CS GROUP
+ * Licensed to CS GROUP (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.los;
+
+import java.util.stream.Stream;
+
+import org.hipparchus.analysis.differentiation.Derivative;
+import org.hipparchus.geometry.euclidean.threed.FieldRotation;
+import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
+import org.hipparchus.geometry.euclidean.threed.Rotation;
+import org.hipparchus.geometry.euclidean.threed.RotationConvention;
+import org.hipparchus.geometry.euclidean.threed.Vector3D;
+import org.hipparchus.util.FastMath;
+import org.orekit.rugged.utils.DerivativeGenerator;
+import org.orekit.utils.ParameterDriver;
+import org.orekit.utils.ParameterObserver;
+
+/** {@link TimeIndependentLOSTransform LOS transform} based on a fixed rotation.
+ * @author Luc Maisonobe
+ * @see LOSBuilder
+ */
+public class FixedRotation implements TimeIndependentLOSTransform {
+
+    /** Parameters scaling factor.
+     * <p>
+     * We use a power of 2 to avoid numeric noise introduction
+     * in the multiplications/divisions sequences.
+     * </p>
+     */
+    private final double SCALE = FastMath.scalb(1.0, -20);
+
+    /** Rotation axis. */
+    private final Vector3D axis;
+
+    /** Underlying rotation. */
+    private Rotation rotation;
+
+    /** Underlying rotation with derivatives. */
+    private FieldRotation<?> rDS;
+
+    /** Driver for rotation angle. */
+    private final ParameterDriver angleDriver;
+
+    /** Simple constructor.
+     * <p>
+     * The single parameter is the rotation angle.
+     * </p>
+     * @param name name of the rotation (used for estimated parameters identification)
+     * @param axis rotation axis
+     * @param angle rotation angle
+     */
+    public FixedRotation(final String name, final Vector3D axis, final double angle) {
+
+        this.axis     = axis;
+        this.rotation = null;
+        this.rDS      = null;
+        this.angleDriver = new ParameterDriver(name, angle, SCALE, -2 * FastMath.PI, 2 * FastMath.PI);
+        angleDriver.addObserver(new ParameterObserver() {
+            @Override
+            public void valueChanged(final double previousValue, final ParameterDriver driver) {
+                // reset rotations to null, they will be evaluated lazily if needed
+                rotation = null;
+                rDS      = null;
+            }
+        });
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public Stream<ParameterDriver> getParametersDrivers() {
+        return Stream.of(angleDriver);
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public Vector3D transformLOS(final int i, final Vector3D los) {
+        if (rotation == null) {
+            // lazy evaluation of the rotation
+            rotation = new Rotation(axis, angleDriver.getValue(), RotationConvention.VECTOR_OPERATOR);
+        }
+        return rotation.applyTo(los);
+    }
+
+    /** {@inheritDoc} */
+    @SuppressWarnings("unchecked")
+    @Override
+    public <T extends Derivative<T>> FieldVector3D<T> transformLOS(final int i, final FieldVector3D<T> los,
+                                                                   final DerivativeGenerator<T> generator) {
+        final FieldRotation<T> rD;
+        if (rDS == null || !rDS.getQ0().getField().equals(generator.getField())) {
+
+            // lazy evaluation of the rotation
+            final FieldVector3D<T> axisDS =
+                            new FieldVector3D<>(generator.constant(axis.getX()),
+                                                generator.constant(axis.getY()),
+                                                generator.constant(axis.getZ()));
+            final T angleDS = generator.variable(angleDriver);
+            rD = new FieldRotation<>(axisDS, angleDS, RotationConvention.VECTOR_OPERATOR);
+
+            // cache evaluated rotation
+            rDS = rD;
+
+        } else {
+            // reuse cached value
+            rD  = (FieldRotation<T>) rDS;
+        }
+
+        return rD.applyTo(los);
+
+    }
+
+}

src/main/java/org/orekit/rugged/los/FixedZHomothety.java

+/* Copyright 2013-2022 CS GROUP
+ * Licensed to CS GROUP (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.los;
+
+import java.util.stream.Stream;
+
+import org.hipparchus.analysis.differentiation.Derivative;
+import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
+import org.hipparchus.geometry.euclidean.threed.Vector3D;
+import org.hipparchus.util.FastMath;
+import org.orekit.rugged.utils.DerivativeGenerator;
+import org.orekit.utils.ParameterDriver;
+import org.orekit.utils.ParameterObserver;
+
+/** {@link TimeIndependentLOSTransform LOS transform} based on a homothety along the Z axis.
+ * @author Lucie Labatallee
+ * @author Guylaine Prat
+ * @see LOSBuilder
+ * @since 2.0
+ */
+public class FixedZHomothety implements TimeIndependentLOSTransform {
+
+    /** Parameters scaling factor.
+     * <p>
+     * We use a power of 2 to avoid numeric noise introduction
+     * in the multiplications/divisions sequences.
+     * </p>
+     */
+    private final double SCALE = FastMath.scalb(1.0, 0);
+
+    /** Homothety factor. */
+    private double factor;
+
+    /** Underlying homothety with derivatives. */
+    private Derivative<?> factorDS;
+
+    /** Driver for homothety factor. */
+    private final ParameterDriver factorDriver;
+
+    /** Simple constructor.
+     * <p>
+     * The single parameter is the homothety factor.
+     * </p>
+     * @param name name of the homothety (used for estimated parameters identification)
+     * @param factorvalue homothety factor
+     */
+    public FixedZHomothety(final String name, final double factorvalue) {
+
+        this.factor   = factorvalue;
+        this.factorDS = null;
+        this.factorDriver = new ParameterDriver(name, factorvalue, SCALE, 0, Double.POSITIVE_INFINITY);
+        factorDriver.addObserver(new ParameterObserver() {
+            @Override
+            public void valueChanged(final double previousValue, final ParameterDriver driver) {
+                // reset factor to zero, they will be evaluated lazily if needed
+                factor = 0.0;
+                factorDS = null;
+            }
+        });
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public Stream<ParameterDriver> getParametersDrivers() {
+        return Stream.of(factorDriver);
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public Vector3D transformLOS(final int i, final Vector3D los) {
+
+        if (factor == 0.0) {
+            // lazy evaluation of the homothety
+            factor = factorDriver.getValue();
+        }
+        return new Vector3D(los.getX(), los.getY(), factor * los.getZ());
+    }
+
+    /** {@inheritDoc} */
+    @SuppressWarnings("unchecked")
+    @Override
+    public <T extends Derivative<T>> FieldVector3D<T> transformLOS(final int i, final FieldVector3D<T> los,
+                                                                   final DerivativeGenerator<T> generator) {
+        final T factorD;
+        if (factorDS == null || !factorDS.getField().equals(generator.getField())) {
+
+            // lazy evaluation of the homothety
+            factorD = generator.variable(factorDriver);
+
+            // cache evaluated homothety
+            factorDS = factorD;
+
+        } else {
+            // reuse cached value
+            factorD  = (T) factorDS;
+        }
+
+        return new FieldVector3D<>(los.getX(), los.getY(), factorD.multiply(los.getZ()));
+
+    }
+
+}

src/main/java/org/orekit/rugged/los/LOSBuilder.java

+/* Copyright 2013-2022 CS GROUP
+ * Licensed to CS GROUP (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.los;
+
+import java.util.ArrayList;
+import java.util.Arrays;
+import java.util.List;
+import java.util.stream.Stream;
+
+import org.hipparchus.analysis.differentiation.Derivative;
+import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
+import org.hipparchus.geometry.euclidean.threed.Vector3D;
+import org.orekit.rugged.utils.DerivativeGenerator;
+import org.orekit.time.AbsoluteDate;
+import org.orekit.utils.ParameterDriver;
+import org.orekit.utils.ParameterObserver;
+
+/** Builder for lines-of-sight list.
+ * <p>
+ * This class implements the <em>builder pattern</em> to create {@link TimeDependentLOS} instances.
+ * It does so by using a <em>fluent API</em> in order to clarify reading and allow
+ * later extensions with new configuration parameters.
+ * </p>
+ * <p>
+ * This builder aims at creating lines-of-sight directions which are
+ * the result of several transforms applied to an initial list of raw
+ * directions. It therefore allows to take into account the optical
+ * path due to mirrors and the alignments of sensors frames with respect
+ * to a spacecraft.
+ * </p>
+ * @see TimeDependentLOS
+ * @see <a href="https://en.wikipedia.org/wiki/Builder_pattern">Builder pattern (wikipedia)</a>
+ * @see <a href="https://en.wikipedia.org/wiki/Fluent_interface">Fluent interface (wikipedia)</a>
+ * @author Luc Maisonobe
+ */
+public class LOSBuilder {
+
+    /** Raw fixed line-of-sights. */
+    private final List<Vector3D> rawLOS;
+
+    /** Transforms to be applied. */
+    private final List<LOSTransform> transforms;
+
+    /** Flag for time-independent only transforms. */
+    private boolean timeIndependent;
+
+    /** Create builder.
+     * @param rawLOS raw fixed lines-of-sight
+     */
+    public LOSBuilder(final List<Vector3D> rawLOS) {
+        this.rawLOS          = rawLOS;
+        this.transforms      = new ArrayList<>();
+        this.timeIndependent = true;
+    }
+
+    /** Add a transform to be applied after the already registered transforms.
+     * @param transform transform to be applied to the lines-of-sight
+     * @return the builder instance
+     */
+    public LOSBuilder addTransform(final TimeIndependentLOSTransform transform) {
+        transforms.add(new TransformAdapter(transform));
+        return this;
+    }
+
+    /** Add a transform to be applied after the already registered transforms.
+     * @param transform transform to be applied to the lines-of-sight
+     * @return the builder instance
+     */
+    public LOSBuilder addTransform(final LOSTransform transform) {
+        transforms.add(transform);
+        timeIndependent = false;
+        return this;
+    }
+
+    /** Build a lines-of-sight provider.
+     * @return lines-of-sight provider
+     */
+    public TimeDependentLOS build() {
+
+        if (timeIndependent) {
+            // fast implementation for time-independent lines-of-sight
+            return new FixedLOS(rawLOS, transforms);
+        } else {
+            // regular implementation, for time-dependent lines-of-sight
+            return new TransformsSequenceLOS(rawLOS, transforms);
+        }
+
+    }
+
+    /** Adapter from time-independent transform to time-dependent transform. */
+    private static class TransformAdapter implements LOSTransform {
+
+        /** Underlying transform. */
+        private final TimeIndependentLOSTransform transform;
+
+        /** Simple constructor.
+         * @param transform underlying time-independent transform
+         */
+        TransformAdapter(final TimeIndependentLOSTransform transform) {
+            this.transform = transform;
+        }
+
+        /** {@inheritDoc} */
+        @Override
+        public Vector3D transformLOS(final int i, final Vector3D los, final AbsoluteDate date) {
+            return transform.transformLOS(i, los);
+        }
+
+        /** {@inheritDoc} */
+        @Override
+        public <T extends Derivative<T>> FieldVector3D<T> transformLOS(final int i, final FieldVector3D<T> los,
+                                                                       final AbsoluteDate date, final DerivativeGenerator<T> generator) {
+            return transform.transformLOS(i, los, generator);
+        }
+
+        /** {@inheritDoc} */
+        @Override
+        public Stream<ParameterDriver> getParametersDrivers() {
+            return transform.getParametersDrivers();
+        }
+
+    }
+
+    /** Implement time-independent LOS by recomputing directions by applying all transforms each time. */
+    private static class TransformsSequenceLOS implements TimeDependentLOS {
+
+        /** Raw direction. */
+        private final Vector3D[] raw;
+
+        /** Transforms to be applied. */
+        private final List<LOSTransform> transforms;
+
+        /** Simple constructor.
+         * @param raw raw directions
+         * @param transforms transforms to apply
+         */
+        TransformsSequenceLOS(final List<Vector3D> raw, final List<LOSTransform> transforms) {
+
+            // copy the lists, to ensure immutability of the built object,
+            // in case addTransform is called again after build
+            // or the raw LOS list is changed by caller
+            this.raw = new Vector3D[raw.size()];
+            for (int i = 0; i < raw.size(); ++i) {
+                this.raw[i] = raw.get(i);
+            }
+
+            this.transforms = new ArrayList<>(transforms);
+
+        }
+
+        /** {@inheritDoc} */
+        public int getNbPixels() {
+            return raw.length;
+        }
+
+        /** {@inheritDoc} */
+        @Override
+        public Vector3D getLOS(final int index, final AbsoluteDate date) {
+            Vector3D los = raw[index];
+            for (final LOSTransform transform : transforms) {
+                los = transform.transformLOS(index, los, date);
+            }
+            return los.normalize();
+        }
+
+        /** {@inheritDoc} */
+        @Override
+        public <T extends Derivative<T>> FieldVector3D<T> getLOSDerivatives(final int index, final AbsoluteDate date,
+                                                                            final DerivativeGenerator<T> generator) {
+
+            // the raw line of sights are considered to be constant
+            FieldVector3D<T> los = new FieldVector3D<>(generator.constant(raw[index].getX()),
+                                                       generator.constant(raw[index].getY()),
+                                                       generator.constant(raw[index].getZ()));
+
+            // apply the transforms, which depend on parameters and hence may introduce non-zero derivatives
+            for (final LOSTransform transform : transforms) {
+                los = transform.transformLOS(index, los, date, generator);
+            }
+
+            return los.normalize();
+
+        }
+
+        @Override
+        public Stream<ParameterDriver> getParametersDrivers() {
+            Stream<ParameterDriver> drivers = Stream.<ParameterDriver>empty();
+            for (final LOSTransform transform : transforms) {
+                drivers = Stream.concat(drivers, transform.getParametersDrivers());
+            }
+            return drivers;
+        }
+
+    }
+
+    /** Implement time-independent LOS by computing directions only when parameters are changed. */
+    private static class FixedLOS extends TransformsSequenceLOS {
+
+        /** transformed direction for los. */
+        private final Vector3D[] transformed;
+
+        /** Simple constructor.
+         * @param raw raw directions
+         * @param transforms transforms to apply (must be time-independent!)
+         */
+        FixedLOS(final List<Vector3D> raw, final List<LOSTransform> transforms) {
+
+            super(raw, transforms);
+            transformed = new Vector3D[raw.size()];
+
+            // we will reset the transforms to null when parameters are changed
+            final ParameterObserver resettingObserver = new ParameterObserver() {
+                /** {@inheritDoc} */
+                @Override
+                public void valueChanged(final double previousValue, final ParameterDriver driver) {
+                    Arrays.fill(transformed, null);
+                }
+            };
+            getParametersDrivers().forEach(driver -> {
+                driver.addObserver(resettingObserver);
+            });
+        }
+
+        /** {@inheritDoc} */
+        @Override
+        public Vector3D getLOS(final int index, final AbsoluteDate date) {
+            if (transformed[index] == null) {
+                // recompute the transformed los direction only if needed
+                transformed[index] = super.getLOS(index, date);
+            }
+            return transformed[index];
+        }
+
+    }
+
+}

src/main/java/org/orekit/rugged/los/LOSTransform.java

+/* Copyright 2013-2022 CS GROUP
+ * Licensed to CS GROUP (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.los;
+
+import java.util.stream.Stream;
+
+import org.hipparchus.analysis.differentiation.Derivative;
+import org.hipparchus.analysis.differentiation.DerivativeStructure;
+import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
+import org.hipparchus.geometry.euclidean.threed.Vector3D;
+import org.orekit.rugged.utils.DerivativeGenerator;
+import org.orekit.time.AbsoluteDate;
+import org.orekit.utils.ParameterDriver;
+
+/** Interface for lines-of-sight transforms.
+ * @author Luc Maisonobe
+ * @see LOSBuilder
+ */
+public interface LOSTransform {
+
+    /** Transform a line-of-sight.
+     * @param i los pixel index
+     * @param los line-of-sight to transform
+     * @param date current date
+     * @return transformed line-of-sight
+     */
+    Vector3D transformLOS(int i, Vector3D los, AbsoluteDate date);
+
+    /** Transform a line-of-sight and its partial derivatives.
+     * <p>
+     * This method is used for LOS calibration purposes. It allows to compute
+     * the Jacobian matrix of the LOS with respect to the parameters, which
+     * are typically polynomials coefficients representing rotation angles.
+     * These polynomials can be used for example to model thermo-elastic effects.
+     * </p>
+     * @param <T> derivative type
+     * @param index los pixel index
+     * @param los line-of-sight to transform
+     * @param date date
+     * @param generator generator to use for building {@link DerivativeStructure} instances
+     * @return line of sight, and its first partial derivatives with respect to the parameters
+     */
+    <T extends Derivative<T>> FieldVector3D<T> transformLOS(int index, FieldVector3D<T> los,
+                                                            AbsoluteDate date, DerivativeGenerator<T> generator);
+
+    /** Get the drivers for LOS parameters.
+     * @return drivers for LOS parameters
+     * @since 2.0
+     */
+    Stream<ParameterDriver> getParametersDrivers();
+
+}

src/main/java/org/orekit/rugged/los/PolynomialRotation.java

+/* Copyright 2013-2022 CS GROUP
+ * Licensed to CS GROUP (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.los;
+
+import java.util.stream.Stream;
+
+import org.hipparchus.Field;
+import org.hipparchus.analysis.differentiation.Derivative;
+import org.hipparchus.analysis.polynomials.PolynomialFunction;
+import org.hipparchus.geometry.euclidean.threed.FieldRotation;
+import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
+import org.hipparchus.geometry.euclidean.threed.Rotation;
+import org.hipparchus.geometry.euclidean.threed.RotationConvention;
+import org.hipparchus.geometry.euclidean.threed.Vector3D;
+import org.hipparchus.util.FastMath;
+import org.hipparchus.util.MathArrays;
+import org.orekit.rugged.utils.DerivativeGenerator;
+import org.orekit.time.AbsoluteDate;
+import org.orekit.utils.ParameterDriver;
+import org.orekit.utils.ParameterObserver;
+
+/** {@link LOSTransform LOS transform} based on a rotation with polynomial angle.
+ * @author Luc Maisonobe
+ * @see LOSBuilder
+ */
+public class PolynomialRotation implements LOSTransform {
+
+    /** Parameters scaling factor.
+     * <p>
+     * We use a power of 2 to avoid numeric noise introduction
+     * in the multiplications/divisions sequences.
+     * </p>
+     */
+    private final double SCALE = FastMath.scalb(1.0, -20);
+
+    /** Rotation axis. */
+    private final Vector3D axis;
+
+    /** Rotation angle polynomial. */
+    private PolynomialFunction angle;
+
+    /** Rotation axis and derivatives. */
+    private FieldVector3D<?> axisDS;
+
+    /** Rotation angle polynomial and derivatives. */
+    private Derivative<?>[] angleDS;
+
+    /** Reference date for polynomial evaluation. */
+    private final AbsoluteDate referenceDate;
+
+    /** Drivers for rotation angle polynomial coefficients. */
+    private final ParameterDriver[] coefficientsDrivers;
+
+    /** Simple constructor.
+     * <p>
+     * The angle of the rotation is evaluated as a polynomial in t,
+     * where t is the duration in seconds between evaluation date and
+     * reference date. The parameters are the polynomial coefficients,
+     * with the constant term at index 0.
+     * </p>
+     * @param name name of the rotation (used for estimated parameters identification)
+     * @param axis rotation axis
+     * @param referenceDate reference date for the polynomial angle
+     * @param angleCoeffs polynomial coefficients of the polynomial angle,
+     * with the constant term at index 0
+     */
+    public PolynomialRotation(final String name,
+                              final Vector3D axis,
+                              final AbsoluteDate referenceDate,
+                              final double... angleCoeffs) {
+        this.axis                = axis;
+        this.referenceDate       = referenceDate;
+        this.coefficientsDrivers = new ParameterDriver[angleCoeffs.length];
+        final ParameterObserver resettingObserver = new ParameterObserver() {
+            @Override
+            public void valueChanged(final double previousValue, final ParameterDriver driver) {
+                // reset rotations to null, they will be evaluated lazily if needed
+                angle   = null;
+                axisDS  = null;
+                angleDS = null;
+            }
+        };
+        for (int i = 0; i < angleCoeffs.length; ++i) {
+            coefficientsDrivers[i] = new ParameterDriver(name + "[" + i + "]", angleCoeffs[i], SCALE,
+                    Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY);
+            coefficientsDrivers[i].addObserver(resettingObserver);
+        }
+    }
+
+    /** Simple constructor.
+     * <p>
+     * The angle of the rotation is evaluated as a polynomial in t,
+     * where t is the duration in seconds between evaluation date and
+     * reference date. The parameters are the polynomial coefficients,
+     * with the constant term at index 0.
+     * </p>
+     * @param name name of the rotation (used for estimated parameters identification)
+     * @param axis rotation axis
+     * @param referenceDate reference date for the polynomial angle
+     * @param angle polynomial angle
+     */
+    public PolynomialRotation(final String name,
+                              final Vector3D axis,
+                              final AbsoluteDate referenceDate,
+                              final PolynomialFunction angle) {
+        this(name, axis, referenceDate, angle.getCoefficients());
+    }
+
+    /** {@inheritDoc}
+     * @since 2.0
+     */
+    @Override
+    public Stream<ParameterDriver> getParametersDrivers() {
+        return Stream.of(coefficientsDrivers);
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public Vector3D transformLOS(final int i, final Vector3D los, final AbsoluteDate date) {
+        if (angle == null) {
+            // lazy evaluation of the rotation
+            final double[] coefficients = new double[coefficientsDrivers.length];
+            for (int k = 0; k < coefficients.length; ++k) {
+                coefficients[k] = coefficientsDrivers[k].getValue();
+            }
+            angle = new PolynomialFunction(coefficients);
+        }
+        return new Rotation(axis,
+                            angle.value(date.durationFrom(referenceDate)),
+                            RotationConvention.VECTOR_OPERATOR).applyTo(los);
+    }
+
+    /** {@inheritDoc} */
+    @SuppressWarnings("unchecked")
+    @Override
+    public <T extends Derivative<T>> FieldVector3D<T> transformLOS(final int i, final FieldVector3D<T> los,
+                                                                   final AbsoluteDate date,
+                                                                   final DerivativeGenerator<T> generator) {
+
+        final Field<T> field = generator.getField();
+        final FieldVector3D<T> axisD;
+        final T[] angleD;
+        if (axisDS == null || !axisDS.getX().getField().equals(field)) {
+
+            // lazy evaluation of the rotation
+            axisD = new FieldVector3D<>(generator.constant(axis.getX()),
+                                        generator.constant(axis.getY()),
+                                        generator.constant(axis.getZ()));
+            angleD = MathArrays.buildArray(field, coefficientsDrivers.length);
+            for (int k = 0; k < angleD.length; ++k) {
+                angleD[k] = generator.variable(coefficientsDrivers[k]);
+            }
+
+            // cache evaluated rotation parameters
+            axisDS  = axisD;
+            angleDS = angleD;
+
+        } else {
+            // reuse cached values
+            axisD  = (FieldVector3D<T>) axisDS;
+            angleD = (T[]) angleDS;
+        }
+
+        // evaluate polynomial, with all its partial derivatives
+        final double t = date.durationFrom(referenceDate);
+        T alpha = field.getZero();
+        for (int k = angleDS.length - 1; k >= 0; --k) {
+            alpha = alpha.multiply(t).add(angleD[k]);
+        }
+
+        return new FieldRotation<>(axisD, alpha, RotationConvention.VECTOR_OPERATOR).applyTo(los);
+
+    }
+
+}

src/main/java/org/orekit/rugged/los/TimeDependentLOS.java

+/* Copyright 2013-2022 CS GROUP
+ * Licensed to CS GROUP (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.los;
+
+import java.util.stream.Stream;
+
+import org.hipparchus.analysis.differentiation.Derivative;
+import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
+import org.hipparchus.geometry.euclidean.threed.Vector3D;
+import org.orekit.rugged.utils.DerivativeGenerator;
+import org.orekit.time.AbsoluteDate;
+import org.orekit.utils.ParameterDriver;
+
+/** Interface representing a line-of-sight which depends on time.
+ * @see org.orekit.rugged.linesensor.LineSensor
+ * @author Luc Maisonobe
+ */
+public interface TimeDependentLOS {
+
+    /** Get the number of pixels.
+     * @return number of pixels
+     */
+    int getNbPixels();
+
+    /** Get the line of sight for a given date.
+     * @param index los pixel index
+     * @param date date
+     * @return line of sight
+     */
+    Vector3D getLOS(int index, AbsoluteDate date);
+
+    /** Get the line of sight and its partial derivatives for a given date.
+     * <p>
+     * This method is used for LOS calibration purposes. It allows to compute
+     * the Jacobian matrix of the LOS with respect to the estimated parameters, which
+     * are typically polynomials coefficients representing rotation angles.
+     * These polynomials can be used for example to model thermo-elastic effects.
+     * </p>
+     * <p>
+     * Note that in order for the partial derivatives to be properly set up, the
+     * {@link org.orekit.utils.ParameterDriver#setSelected(boolean) setSelected}
+     * method must have been set to {@code true} for the various parameters returned
+     * by {@link #getParametersDrivers()} that should be estimated.
+     * </p>
+     * @param <T> derivative type
+     * @param index los pixel index
+     * @param date date
+     * @param generator generator to use for building {@link Derivative} instances
+     * @return line of sight, and its first partial derivatives with respect to the parameters
+     * @since 2.0
+     */
+    <T extends Derivative<T>> FieldVector3D<T> getLOSDerivatives(int index, AbsoluteDate date,
+                                                                 DerivativeGenerator<T> generator);
+
+    /** Get the drivers for LOS parameters.
+     * @return drivers for LOS parameters
+     * @since 2.0
+     */
+    Stream<ParameterDriver> getParametersDrivers();
+
+}

src/main/java/org/orekit/rugged/los/TimeIndependentLOSTransform.java

+/* Copyright 2013-2022 CS GROUP
+ * Licensed to CS GROUP (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.los;
+
+import java.util.stream.Stream;
+
+import org.hipparchus.analysis.differentiation.Derivative;
+import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
+import org.hipparchus.geometry.euclidean.threed.Vector3D;
+import org.orekit.rugged.utils.DerivativeGenerator;
+import org.orekit.utils.ParameterDriver;
+
+/** Interface for lines-of-sight tranforms that do not depend on time.
+ * @author Luc Maisonobe
+ * @see LOSBuilder
+ */
+public interface TimeIndependentLOSTransform {
+
+    /** Transform a line-of-sight.
+     * @param i los pixel index
+     * @param los line-of-sight to transform
+     * @return transformed line-of-sight
+     */
+    Vector3D transformLOS(int i, Vector3D los);
+
+    /** Transform a line-of-sight and its partial derivatives.
+     * <p>
+     * This method is used for LOS calibration purposes. It allows to compute
+     * the Jacobian matrix of the LOS with respect to the parameters, which
+     * are typically polynomials coefficients representing rotation angles.
+     * These polynomials can be used for example to model thermo-elastic effects.
+     * </p>
+     * <p>
+     * Note that in order for the partial derivatives to be properly set up, the
+     * {@link org.orekit.utils.ParameterDriver#setSelected(boolean) setSelected}
+     * method must have been set to {@code true} for the various parameters returned
+     * by {@link #getParametersDrivers()} that should be estimated.
+     * </p>
+     * @param <T> derivative type
+     * @param index los pixel index
+     * @param los line-of-sight to transform
+     * @param generator generator to use for building {@link Derivative} instances
+     * @return line of sight, and its first partial derivatives with respect to the parameters
+     */
+    <T extends Derivative<T>> FieldVector3D<T> transformLOS(int index, FieldVector3D<T> los,
+                                                            DerivativeGenerator<T> generator);
+
+    /** Get the drivers for LOS parameters.
+     * @return drivers for LOS parameters
+     * @since 2.0
+     */
+    Stream<ParameterDriver> getParametersDrivers();
+
+}

src/main/java/org/orekit/rugged/los/package-info.java

+/* Copyright 2013-2022 CS GROUP
+ * Licensed to CS GROUP (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.
+ */
+/**
+ *
+ * This package provides all the tools to build lines-of-sight (LOS).
+ *
+ * @author Luc Maisonobe
+ * @author Guylaine Prat
+ *
+ */
+package org.orekit.rugged.los;

src/main/java/org/orekit/rugged/core/raster/SimpleTile.java → src/main/java/org/orekit/rugged/raster/SimpleTile.java

-/* Copyright 2013-2014 CS Systèmes d'Information
- * Licensed to CS Systèmes d'Information (CS) under one or more
+/* Copyright 2013-2022 CS GROUP
+ * Licensed to CS GROUP (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
@@ -14,23 +14,29 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
-package org.orekit.rugged.core.raster;
+package org.orekit.rugged.raster;
 
-import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
-import org.apache.commons.math3.util.FastMath;
-import org.apache.commons.math3.util.Precision;
-import org.orekit.bodies.GeodeticPoint;
-import org.orekit.rugged.api.RuggedException;
-import org.orekit.rugged.api.RuggedMessages;
+import java.util.Arrays;
+
+import org.hipparchus.geometry.euclidean.threed.Vector3D;
+import org.hipparchus.util.FastMath;
+import org.hipparchus.util.Precision;
+import org.orekit.rugged.errors.DumpManager;
+import org.orekit.rugged.errors.RuggedException;
+import org.orekit.rugged.errors.RuggedMessages;
+import org.orekit.rugged.utils.MaxSelector;
+import org.orekit.rugged.utils.MinSelector;
+import org.orekit.rugged.utils.NormalizedGeodeticPoint;
 
 
 /** Simple implementation of a {@link Tile}.
  * @see SimpleTileFactory
  * @author Luc Maisonobe
+ * @author Guylaine Prat
  */
 public class SimpleTile implements Tile {
 
-    /** Tolerance used to interpolate points slightly out of tile (in pixels). */
+    /** Tolerance used to interpolate points slightly out of tile (in cells). */
     private static final double TOLERANCE = 1.0 / 8.0;
 
     /** Minimum latitude. */
@@ -54,9 +60,21 @@ public class SimpleTile implements Tile {
     /** Minimum elevation. */
     private double minElevation;
 
+    /** Latitude index of min elevation. */
+    private int minElevationLatitudeIndex;
+
+    /** Longitude index of min elevation. */
+    private int minElevationLongitudeIndex;
+
     /** Maximum elevation. */
     private double maxElevation;
 
+    /** Latitude index of max elevation. */
+    private int maxElevationLatitudeIndex;
+
+    /** Longitude index of max elevation. */
+    private int maxElevationLongitudeIndex;
+
     /** Elevation array. */
     private double[] elevations;
 
@@ -72,27 +90,31 @@ public class SimpleTile implements Tile {
     @Override
     public void setGeometry(final double newMinLatitude, final double newMinLongitude,
                             final double newLatitudeStep, final double newLongitudeStep,
-                            final int newLatitudeRows, final int newLongitudeColumns)
-        throws RuggedException {
-        this.minLatitude      = newMinLatitude;
-        this.minLongitude     = newMinLongitude;
-        this.latitudeStep     = newLatitudeStep;
-        this.longitudeStep    = newLongitudeStep;
-        this.latitudeRows     = newLatitudeRows;
-        this.longitudeColumns = newLongitudeColumns;
-        this.minElevation     = Double.POSITIVE_INFINITY;
-        this.maxElevation     = Double.NEGATIVE_INFINITY;
+                            final int newLatitudeRows, final int newLongitudeColumns) {
+        this.minLatitude                = newMinLatitude;
+        this.minLongitude               = newMinLongitude;
+        this.latitudeStep               = newLatitudeStep;
+        this.longitudeStep              = newLongitudeStep;
+        this.latitudeRows               = newLatitudeRows;
+        this.longitudeColumns           = newLongitudeColumns;
+        this.minElevation               = Double.POSITIVE_INFINITY;
+        this.minElevationLatitudeIndex  = -1;
+        this.minElevationLongitudeIndex = -1;
+        this.maxElevation               = Double.NEGATIVE_INFINITY;
+        this.maxElevationLatitudeIndex  = -1;
+        this.maxElevationLongitudeIndex = -1;
 
         if (newLatitudeRows < 1 || newLongitudeColumns < 1) {
             throw new RuggedException(RuggedMessages.EMPTY_TILE, newLatitudeRows, newLongitudeColumns);
         }
         this.elevations = new double[newLatitudeRows * newLongitudeColumns];
+        Arrays.fill(elevations, Double.NaN);
 
     }
 
     /** {@inheritDoc} */
     @Override
-    public void tileUpdateCompleted() throws RuggedException {
+    public void tileUpdateCompleted() {
         processUpdatedElevation(elevations);
     }
 
@@ -174,6 +196,18 @@ public class SimpleTile implements Tile {
         return minElevation;
     }
 
+    /** {@inheritDoc} */
+    @Override
+    public int getMinElevationLatitudeIndex() {
+        return minElevationLatitudeIndex;
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public int getMinElevationLongitudeIndex() {
+        return minElevationLongitudeIndex;
+    }
+
     /** {@inheritDoc} */
     @Override
     public double getMaxElevation() {
@@ -182,34 +216,55 @@ public class SimpleTile implements Tile {
 
     /** {@inheritDoc} */
     @Override
-    public void setElevation(final int latitudeIndex, final int longitudeIndex,
-                             final double elevation) throws RuggedException {
+    public int getMaxElevationLatitudeIndex() {
+        return maxElevationLatitudeIndex;
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public int getMaxElevationLongitudeIndex() {
+        return maxElevationLongitudeIndex;
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public void setElevation(final int latitudeIndex, final int longitudeIndex, final double elevation) {
+
         if (latitudeIndex  < 0 || latitudeIndex  > (latitudeRows - 1) ||
             longitudeIndex < 0 || longitudeIndex > (longitudeColumns - 1)) {
             throw new RuggedException(RuggedMessages.OUT_OF_TILE_INDICES,
                                       latitudeIndex, longitudeIndex,
                                       latitudeRows - 1, longitudeColumns - 1);
         }
-        minElevation = FastMath.min(minElevation, elevation);
-        maxElevation = FastMath.max(maxElevation, elevation);
+        if (MinSelector.getInstance().selectFirst(elevation, minElevation)) {
+            minElevation               = elevation;
+            minElevationLatitudeIndex  = latitudeIndex;
+            minElevationLongitudeIndex = longitudeIndex;
+        }
+        if (MaxSelector.getInstance().selectFirst(elevation, maxElevation)) {
+            maxElevation               = elevation;
+            maxElevationLatitudeIndex  = latitudeIndex;
+            maxElevationLongitudeIndex = longitudeIndex;
+        }
         elevations[latitudeIndex * getLongitudeColumns() + longitudeIndex] = elevation;
     }
 
     /** {@inheritDoc} */
     @Override
     public double getElevationAtIndices(final int latitudeIndex, final int longitudeIndex) {
-        return elevations[latitudeIndex * getLongitudeColumns() + longitudeIndex];
+        final double elevation = elevations[latitudeIndex * getLongitudeColumns() + longitudeIndex];
+        DumpManager.dumpTileCell(this, latitudeIndex, longitudeIndex, elevation);
+        return elevation;
     }
 
     /** {@inheritDoc}
      * <p>
-     * This classes uses an arbitrary 1/8 pixel tolerance for interpolating
+     * This classes uses an arbitrary 1/8 cell tolerance for interpolating
      * slightly out of tile points.
      * </p>
      */
     @Override
-    public double interpolateElevation(final double latitude, final double longitude)
-        throws RuggedException {
+    public double interpolateElevation(final double latitude, final double longitude) {
 
         final double doubleLatitudeIndex  = getDoubleLatitudeIndex(latitude);
         final double doubleLongitudeIndex = getDoubleLontitudeIndex(longitude);
@@ -238,6 +293,7 @@ public class SimpleTile implements Tile {
         final double e10  = getElevationAtIndices(latitudeIndex,     longitudeIndex + 1);
         final double e01  = getElevationAtIndices(latitudeIndex + 1, longitudeIndex);
         final double e11  = getElevationAtIndices(latitudeIndex + 1, longitudeIndex + 1);
+
         return (e00 * (1.0 - dLon) + dLon * e10) * (1.0 - dLat) +
                (e01 * (1.0 - dLon) + dLon * e11) * dLat;
 
@@ -245,22 +301,31 @@ public class SimpleTile implements Tile {
 
     /** {@inheritDoc} */
     @Override
-    public GeodeticPoint pixelIntersection(final GeodeticPoint p, final Vector3D los,
-                                           final int latitudeIndex, final int longitudeIndex)
-        throws RuggedException {
-
-        // Digital Elevation Mode coordinates at pixel vertices
-        final double x00 = getLongitudeAtIndex(longitudeIndex);
-        final double y00 = getLatitudeAtIndex(latitudeIndex);
-        final double z00 = getElevationAtIndices(latitudeIndex,     longitudeIndex);
-        final double z01 = getElevationAtIndices(latitudeIndex + 1, longitudeIndex);
-        final double z10 = getElevationAtIndices(latitudeIndex,     longitudeIndex + 1);
-        final double z11 = getElevationAtIndices(latitudeIndex + 1, longitudeIndex + 1);
+    public NormalizedGeodeticPoint cellIntersection(final NormalizedGeodeticPoint p, final Vector3D los,
+                                                    final int latitudeIndex, final int longitudeIndex) {
+
+        // ensure neighboring cells to not fall out of tile
+        final int iLat  = FastMath.max(0, FastMath.min(latitudeRows     - 2, latitudeIndex));
+        final int jLong = FastMath.max(0, FastMath.min(longitudeColumns - 2, longitudeIndex));
+
+        // Digital Elevation Mode coordinates at cell vertices
+        final double x00 = getLongitudeAtIndex(jLong);
+        final double y00 = getLatitudeAtIndex(iLat);
+        final double z00 = getElevationAtIndices(iLat,     jLong);
+        final double z01 = getElevationAtIndices(iLat + 1, jLong);
+        final double z10 = getElevationAtIndices(iLat,     jLong + 1);
+        final double z11 = getElevationAtIndices(iLat + 1, jLong + 1);
+
+        // normalize back to tile coordinates
+        final NormalizedGeodeticPoint tileP = new NormalizedGeodeticPoint(p.getLatitude(),
+                                                                          p.getLongitude(),
+                                                                          p.getAltitude(),
+                                                                          x00);
 
         // line-of-sight coordinates at close points
-        final double dxA = (p.getLongitude() - x00) / longitudeStep;
-        final double dyA = (p.getLatitude()  - y00) / latitudeStep;
-        final double dzA = p.getAltitude();
+        final double dxA = (tileP.getLongitude() - x00) / longitudeStep;
+        final double dyA = (tileP.getLatitude()  - y00) / latitudeStep;
+        final double dzA = tileP.getAltitude();
         final double dxB = dxA + los.getX() / longitudeStep;
         final double dyB = dyA + los.getY() / latitudeStep;
         final double dzB = dzA + los.getZ();
@@ -308,8 +373,8 @@ public class SimpleTile implements Tile {
 
         }
 
-        final GeodeticPoint p1 = interpolate(t1, p, dxA, dyA, los);
-        final GeodeticPoint p2 = interpolate(t2, p, dxA, dyA, los);
+        final NormalizedGeodeticPoint p1 = interpolate(t1, tileP, dxA, dyA, los, x00);
+        final NormalizedGeodeticPoint p2 = interpolate(t2, tileP, dxA, dyA, los, x00);
 
         // select the first point along line-of-sight
         if (p1 == null) {
@@ -324,15 +389,17 @@ public class SimpleTile implements Tile {
 
     /** Interpolate point along a line.
      * @param t abscissa along the line
-     * @param p start point
-     * @param dxP relative coordinate of the start point with respect to current pixel
-     * @param dyP relative coordinate of the start point with respect to current pixel
+     * @param tileP start point, normalized to tile area
+     * @param dxP relative coordinate of the start point with respect to current cell
+     * @param dyP relative coordinate of the start point with respect to current cell
      * @param los direction of the line-of-sight, in geodetic space
+     * @param centralLongitude reference longitude lc such that the point longitude will
+     * be normalized between lc-π and lc+π
      * @return interpolated point along the line
      */
-    private GeodeticPoint interpolate(final double t, final GeodeticPoint p,
-                                      final double dxP, final double dyP,
-                                      final Vector3D los) {
+    private NormalizedGeodeticPoint interpolate(final double t, final NormalizedGeodeticPoint tileP,
+                                                final double dxP, final double dyP,
+                                                final Vector3D los, final double centralLongitude) {
 
         if (Double.isInfinite(t)) {
             return null;
@@ -341,9 +408,10 @@ public class SimpleTile implements Tile {
         final double dx = dxP + t * los.getX() / longitudeStep;
         final double dy = dyP + t * los.getY() / latitudeStep;
         if (dx >= -TOLERANCE && dx <= 1 + TOLERANCE && dy >= -TOLERANCE && dy <= 1 + TOLERANCE) {
-            return new GeodeticPoint(p.getLatitude()  + t * los.getY(),
-                                     p.getLongitude() + t * los.getX(),
-                                     p.getAltitude()  + t * los.getZ());
+            return new NormalizedGeodeticPoint(tileP.getLatitude()  + t * los.getY(),
+                                               tileP.getLongitude() + t * los.getX(),
+                                               tileP.getAltitude()  + t * los.getZ(),
+                                               centralLongitude);
         } else {
             return null;
         }
@@ -352,26 +420,26 @@ public class SimpleTile implements Tile {
 
     /** {@inheritDoc} */
     @Override
-    public int getLatitudeIndex(final double latitude) {
+    public int getFloorLatitudeIndex(final double latitude) {
         return (int) FastMath.floor(getDoubleLatitudeIndex(latitude));
     }
 
     /** {@inheritDoc} */
     @Override
-    public int getLongitudeIndex(final double longitude) {
+    public int getFloorLongitudeIndex(final double longitude) {
         return (int) FastMath.floor(getDoubleLontitudeIndex(longitude));
     }
 
     /** Get the latitude index of a point.
-     * @param latitude geodetic latitude
-     * @return latitute index (it may lie outside of the tile!)
+     * @param latitude geodetic latitude (rad)
+     * @return latitude index (it may lie outside of the tile!)
      */
     private double getDoubleLatitudeIndex(final double latitude) {
         return (latitude  - minLatitude)  / latitudeStep;
     }
 
     /** Get the longitude index of a point.
-     * @param longitude geodetic latitude
+     * @param longitude geodetic longitude (rad)
      * @return longitude index (it may lie outside of the tile!)
      */
     private double getDoubleLontitudeIndex(final double longitude) {
@@ -381,8 +449,8 @@ public class SimpleTile implements Tile {
     /** {@inheritDoc} */
     @Override
     public Location getLocation(final double latitude, final double longitude) {
-        final int latitudeIndex  = getLatitudeIndex(latitude);
-        final int longitudeIndex = getLongitudeIndex(longitude);
+        final int latitudeIndex  = getFloorLatitudeIndex(latitude);
+        final int longitudeIndex = getFloorLongitudeIndex(longitude);
         if (longitudeIndex < 0) {
             if (latitudeIndex < 0) {
                 return Location.SOUTH_WEST;
@@ -395,7 +463,7 @@ public class SimpleTile implements Tile {
             if (latitudeIndex < 0) {
                 return Location.SOUTH;
             } else if (latitudeIndex <= (latitudeRows - 2)) {
-                return Location.IN_TILE;
+                return Location.HAS_INTERPOLATION_NEIGHBORS;
             } else {
                 return Location.NORTH;
             }

src/main/java/org/orekit/rugged/core/raster/SimpleTileFactory.java → src/main/java/org/orekit/rugged/raster/SimpleTileFactory.java

-/* Copyright 2013-2014 CS Systèmes d'Information
- * Licensed to CS Systèmes d'Information (CS) under one or more
+/* Copyright 2013-2022 CS GROUP
+ * Licensed to CS GROUP (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
@@ -14,9 +14,7 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
-package org.orekit.rugged.core.raster;
-
-import org.orekit.rugged.core.raster.TileFactory;
+package org.orekit.rugged.raster;
 
 /** Simple implementation of a {@link TileFactory} for {@link SimpleTile}.
  * @author Luc Maisonobe

src/main/java/org/orekit/rugged/core/raster/Tile.java → src/main/java/org/orekit/rugged/raster/Tile.java

-/* Copyright 2013-2014 CS Systèmes d'Information
- * Licensed to CS Systèmes d'Information (CS) under one or more
+/* Copyright 2013-2022 CS GROUP
+ * Licensed to CS GROUP (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
@@ -14,95 +14,168 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
-package org.orekit.rugged.core.raster;
+package org.orekit.rugged.raster;
 
-import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
-import org.orekit.bodies.GeodeticPoint;
-import org.orekit.rugged.api.RuggedException;
-import org.orekit.rugged.api.UpdatableTile;
+import org.hipparchus.geometry.euclidean.threed.Vector3D;
+import org.orekit.rugged.utils.NormalizedGeodeticPoint;
 
 /** Interface representing a raster tile.
+ * <p>
+ * The elevations are considered to be at the <em>center</em> of each cells.
+ * The minimum latitude and longitude hence correspond to the <em>center</em>
+ * of the most South-West cell, and the maximum latitude and longitude
+ * correspond to the <em>center</em> of the most North-East cell.
+ * </p>
  * @author Luc Maisonobe
+ * @author Guylaine Prat
  */
 public interface Tile extends UpdatableTile {
 
-    /** Enumerate for point location with respect to tile. */
+    /** Enumerate for point location with respect to the interpolation grid of a tile.
+     * <p>
+     * Elevations in a tile are interpolated using the four neighboring points
+     * in a grid: (i, j), (i+1, j), (i, j+1), (i+1), (j+1). This implies that a point
+     * can be interpolated only if the elevation for these four points is available
+     * in the tile. A consequence is that a point in the northernmost row (resp.
+     * easternmost column) miss neighboring points at row j+1 (resp. neighboring points
+     * at column i+1) and therefore cannot be interpolated.
+     * </p>
+     * <p>
+     * This enumerate represent the position of a point taking this off-by-one property
+     * into account, the value {@link #HAS_INTERPOLATION_NEIGHBORS} correspond to points that
+     * do have the necessary four neightbors, whereas the other values correspond to points
+     * that are either completely outside of the tile or within the tile but in either the
+     * northernmost row or easternmost column.
+     * </p>
+     */
     enum Location {
 
-        /** Location for points out of tile, past the South-West corner. */
+        /** Location for points out of tile interpolation grid, in the South-West corner direction. */
         SOUTH_WEST,
 
-        /** Location for points out of tile, past the West edge. */
+        /** Location for points out of tile interpolation grid, in the West edge direction. */
         WEST,
 
-        /** Location for points out of tile, past the North-West corner. */
+        /** Location for points out of tile interpolation grid, in the North-West corner direction.
+         * <p>
+         * The point may still be in the tile, but in the northernmost row thus missing required
+         * interpolation points.
+         * </p>
+         */
         NORTH_WEST,
 
-        /** Location for points out of tile, past the North edge. */
+        /** Location for points out of tile interpolation grid, in the North edge direction.
+         * <p>
+         * The point may still be in the tile, but in the northernmost row thus missing required
+         * interpolation points.
+         * </p>
+         */
         NORTH,
 
-        /** Location for points out of tile, past the North-East corner. */
+        /** Location for points out of tile interpolation grid, in the North-East corner direction.
+         * <p>
+         * The point may still be in the tile, but either in the northernmost row or in the
+         * easternmost column thus missing required interpolation points.
+         * </p>
+         */
         NORTH_EAST,
 
-        /** Location for points out of tile, past the East edge. */
+        /** Location for points out of tile interpolation grid, in the East edge direction.
+         * <p>
+         * The point may still be in the tile, but in the easternmost column thus missing required
+         * interpolation points.
+         * </p>
+         */
         EAST,
 
-        /** Location for points out of tile, past the South-East corner. */
+        /** Location for points out of tile interpolation grid, in the South-East corner direction.
+         * <p>
+         * The point may still be in the tile, but in the easternmost column thus missing required
+         * interpolation points.
+         * </p>
+         */
         SOUTH_EAST,
 
-        /** Location for points out of tile, past the South edge. */
+        /** Location for points out of tile interpolation grid, in the South edge direction. */
         SOUTH,
 
-        /** Location for points within tile. */
-        IN_TILE
+        /** Location for points that do have interpolation neighbors.
+         * <p>
+         * The value corresponds to points that can be interpolated using their four
+         * neighboring points in the grid at indices (i, j), (i+1, j), (i, j+1), (i+1),
+         * (j+1). This implies that these points are neither in the northernmost latitude
+         * row nor in the easternmost longitude column.
+         * </p>
+         */
+        HAS_INTERPOLATION_NEIGHBORS
 
     }
 
     /** Hook called at the end of tile update completion.
-     * @exception RuggedException if something wrong occurs
-     * (missing data ...)
      */
-    void tileUpdateCompleted() throws RuggedException;
+    void tileUpdateCompleted();
 
-    /** Get minimum latitude.
-     * @return minimum latitude
+    /** Get minimum latitude of grid interpolation points.
+     * @return minimum latitude of grid interpolation points (rad)
+     * (latitude of the center of the cells of South row)
      */
     double getMinimumLatitude();
 
     /** Get the latitude at some index.
      * @param latitudeIndex latitude index
-     * @return latitude at the specified index
+     * @return latitude at the specified index (rad)
+     * (latitude of the center of the cells of specified row)
      */
     double getLatitudeAtIndex(int latitudeIndex);
 
     /** Get maximum latitude.
-     * @return maximum latitude
+     * <p>
+     * Beware that as a point at maximum latitude is the northernmost
+     * one of the grid, it doesn't have a northwards neighbor and
+     * therefore calling {@link #getLocation(double, double) getLocation}
+     * on such a latitude will return either {@link Location#NORTH_WEST},
+     * {@link Location#NORTH} or {@link Location#NORTH_EAST}, but can
+     * <em>never</em> return {@link Location#HAS_INTERPOLATION_NEIGHBORS}!
+     * </p>
+     * @return maximum latitude (rad)
+     * (latitude of the center of the cells of North row)
      */
     double getMaximumLatitude();
 
     /** Get minimum longitude.
-     * @return minimum longitude
+     * @return minimum longitude (rad)
+     * (longitude of the center of the cells of West column)
      */
     double getMinimumLongitude();
 
     /** Get the longitude at some index.
      * @param longitudeIndex longitude index
-     * @return longitude at the specified index
+     * @return longitude at the specified index (rad)
+     * (longitude of the center of the cells of specified column)
      */
     double getLongitudeAtIndex(int longitudeIndex);
 
     /** Get maximum longitude.
-     * @return maximum longitude
+     * <p>
+     * Beware that as a point at maximum longitude is the easternmost
+     * one of the grid, it doesn't have an eastwards neighbor and
+     * therefore calling {@link #getLocation(double, double) getLocation}
+     * on such a longitude will return either {@link Location#SOUTH_EAST},
+     * {@link Location#EAST} or {@link Location#NORTH_EAST}, but can
+     * <em>never</em> return {@link Location#HAS_INTERPOLATION_NEIGHBORS}!
+     * </p>
+     * @return maximum longitude (rad)
+     * (longitude of the center of the cells of East column)
      */
     double getMaximumLongitude();
 
     /** Get step in latitude (size of one raster element).
-     * @return step in latitude
+     * @return step in latitude (rad)
      */
     double getLatitudeStep();
 
     /** Get step in longitude (size of one raster element).
-     * @return step in longitude
+     * @return step in longitude (rad)
      */
     double getLongitudeStep();
 
@@ -116,67 +189,95 @@ public interface Tile extends UpdatableTile {
      */
     int getLongitudeColumns();
 
-    /** Get the latitude index of a point.
+    /** Get the floor latitude index of a point.
+     * <p>
+     * The specified latitude is always between index and index+1.
+     * </p>
      * @param latitude geodetic latitude
-     * @return latirute index (it may lie outside of the tile!)
+     * @return floor latitude index (it may lie outside of the tile!)
      */
-    int getLatitudeIndex(double latitude);
+    int getFloorLatitudeIndex(double latitude);
 
-    /** Get the longitude index of a point.
-     * @param longitude geodetic latitude
-     * @return longitude index (it may lie outside of the tile!)
+    /** Get the floor longitude index of a point.
+     * <p>
+     * The specified longitude is always between index and index+1.
+     * </p>
+     * @param longitude geodetic longitude
+     * @return floor longitude index (it may lie outside of the tile!)
      */
-    int getLongitudeIndex(double longitude);
+    int getFloorLongitudeIndex(double longitude);
 
     /** Get the minimum elevation in the tile.
-     * @return minimum elevation in the tile
+     * @return minimum elevation in the tile (m)
      */
     double getMinElevation();
 
-    /** Get the maximum elevation in the tile.
-     * @return maximum elevation in the tile
+    /** Get the latitude index of min elevation.
+     * @return latitude index of min elevation*/
+    int getMinElevationLatitudeIndex();
+
+    /** Get the longitude index of min elevation.
+     * @return longitude index of min elevation*/
+    int getMinElevationLongitudeIndex();
+
+   /** Get the maximum elevation in the tile.
+     * @return maximum elevation in the tile (m)
      */
     double getMaxElevation();
 
+    /** Get the latitude index of max elevation.
+     * @return latitude index of max elevation*/
+    int getMaxElevationLatitudeIndex();
+
+    /** Get the longitude index of max elevation.
+     * @return longitude index of max elevation*/
+    int getMaxElevationLongitudeIndex();
+
     /** Get the elevation of an exact grid point.
      * @param latitudeIndex grid point index along latitude
      * @param longitudeIndex grid point index along longitude
-     * @return elevation at grid point
-     * @exception RuggedException if indices are out of bound
+     * @return elevation at grid point (m)
      */
-    double getElevationAtIndices(int latitudeIndex, int longitudeIndex)
-        throws RuggedException;
+    double getElevationAtIndices(int latitudeIndex, int longitudeIndex);
 
     /** Interpolate elevation.
      * <p>
      * In order to cope with numerical accuracy issues when computing
-     * points at tile boundary, a slight tolerance (typically 1/8 pixel)
+     * points at tile boundary, a slight tolerance (typically 1/8 cell)
      * around the tile is allowed. Elevation can therefore be interpolated
      * (really extrapolated in this case) even for points slightly overshooting
-     * tile boundaries, using the closest tile pixel. Attempting to interpolate
+     * tile boundaries, using the closest tile cell. Attempting to interpolate
      * too far from the tile will trigger an exception.
      * </p>
      * @param latitude ground point latitude
      * @param longitude ground point longitude
-     * @return interpolated elevation
-     * @exception RuggedException if point is farthest from the tile than the tolerance
+     * @return interpolated elevation (m)
      */
-    double interpolateElevation(double latitude, double longitude)
-        throws RuggedException;
+    double interpolateElevation(double latitude, double longitude);
 
-    /** Find the intersection of a line-of-sight and a Digital Elevation Model pixel.
-     * @param p point on the line
+    /** Find the intersection of a line-of-sight and a Digital Elevation Model cell.
+     * <p>
+     * Beware that for continuity reasons, the point argument in {@code cellIntersection} is normalized
+     * with respect to other points used by the caller. This implies that the longitude may be
+     * outside of the [-π ; +π] interval (or the [0 ; 2π] interval, depending on the DEM). In particular,
+     * when a Line Of Sight crosses the antimeridian at  ±π longitude, the library may call the
+     * {@code cellIntersection} method with a point having a longitude of -π-ε to ensure this continuity.
+     * As tiles are stored with longitude clipped to a some DEM specific interval (either  [-π ; +π] or [0 ; 2π]),
+     * implementations MUST take care to clip the input point back to the tile interval using
+     * {@link org.hipparchus.util.MathUtils#normalizeAngle(double, double) MathUtils.normalizeAngle(p.getLongitude(),
+     * someLongitudeWithinTheTile)}. The output point normalization should also be made consistent with
+     * the current tile.
+     * </p>
+     * @param p point on the line (beware its longitude is <em>not</em> normalized with respect to tile)
      * @param los line-of-sight, in the topocentric frame (East, North, Zenith) of the point,
      * scaled to match radians in the horizontal plane and meters along the vertical axis
-     * @param latitudeIndex latitude index of the Digital Elevation Model pixel
-     * @param longitudeIndex longitude index of the Digital Elevation Model pixel
+     * @param latitudeIndex latitude index of the Digital Elevation Model cell
+     * @param longitudeIndex longitude index of the Digital Elevation Model cell
      * @return point corresponding to line-of-sight crossing the Digital Elevation Model surface
-     * if it lies within the pixel, null otherwise
-     * @exception RuggedException if intersection point cannot be computed
+     * if it lies within the cell, null otherwise
      */
-    GeodeticPoint pixelIntersection(GeodeticPoint p, Vector3D los,
-                                    int latitudeIndex, int longitudeIndex)
-        throws RuggedException;
+    NormalizedGeodeticPoint cellIntersection(NormalizedGeodeticPoint p, Vector3D los,
+                                             int latitudeIndex, int longitudeIndex);
 
     /** Check if a tile covers a ground point.
      * @param latitude ground point latitude

src/main/java/org/orekit/rugged/core/raster/TileFactory.java → src/main/java/org/orekit/rugged/raster/TileFactory.java

-/* Copyright 2013-2014 CS Systèmes d'Information
- * Licensed to CS Systèmes d'Information (CS) under one or more
+/* Copyright 2013-2022 CS GROUP
+ * Licensed to CS GROUP (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
@@ -14,7 +14,7 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
-package org.orekit.rugged.core.raster;
+package org.orekit.rugged.raster;
 
 
 /** Interface representing a factory for raster tile.