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

-/* Copyright 2013-2016 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

src/main/java/org/orekit/rugged/raster/TileUpdater.java

-/* Copyright 2013-2016 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
@@ -16,15 +16,14 @@
  */
 package org.orekit.rugged.raster;
 
-import org.orekit.rugged.errors.RuggedException;
-
 /** Interface used to update Digital Elevation Model tiles.
  * <p>
- * Implementations of this interface are must be provided by
+ * Implementations of this interface must be provided by
  * the image processing mission-specific layer, thus allowing
  * the Rugged library to access the Digital Elevation Model data.
  * </p>
- *  * @author Luc Maisonobe
+ * @author Luc Maisonobe
+ * @author Guylaine Prat
  */
 public interface TileUpdater {
 
@@ -76,12 +75,10 @@ public interface TileUpdater {
            latitudeStep, longitudeStep, latitudeRows, longitudeColumns)} call.
      *   </li>
      * </ul>
-     * @param latitude latitude that must be covered by the tile
-     * @param longitude longitude that must be covered by the tile
+     * @param latitude latitude that must be covered by the tile (rad)
+     * @param longitude longitude that must be covered by the tile (rad)
      * @param tile to update
-     * @exception RuggedException if tile cannot be updated
      */
-    void updateTile(double latitude, double longitude, UpdatableTile tile)
-        throws RuggedException;
+    void updateTile(double latitude, double longitude, UpdatableTile tile);
 
 }

src/main/java/org/orekit/rugged/raster/TilesCache.java

-/* Copyright 2013-2016 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
@@ -19,6 +19,7 @@ package org.orekit.rugged.raster;
 import org.hipparchus.util.FastMath;
 import java.lang.reflect.Array;
 
+import org.orekit.rugged.errors.DumpManager;
 import org.orekit.rugged.errors.RuggedException;
 import org.orekit.rugged.errors.RuggedMessages;
 
@@ -57,10 +58,8 @@ public class TilesCache<T extends Tile> {
      * @param latitude ground point latitude
      * @param longitude ground point longitude
      * @return tile covering the ground point
-     * @exception RuggedException if newly created tile cannot be updated
      */
-    public T getTile(final double latitude, final double longitude)
-        throws RuggedException {
+    public T getTile(final double latitude, final double longitude) {
 
         for (int i = 0; i < tiles.length; ++i) {
             final T tile = tiles[i];
@@ -88,7 +87,16 @@ public class TilesCache<T extends Tile> {
 
         // create the tile and retrieve its data
         final T tile = factory.createTile();
+
+        // In case dump is asked for, suspend the dump manager as we don't need to dump anything here
+        // For instance for SRTM DEM, the user needs to read Geoid data that are not useful in the dump
+        final Boolean wasSuspended = DumpManager.suspend();
+
         updater.updateTile(latitude, longitude, tile);
+
+        // Resume the dump manager if necessary
+        DumpManager.resume(wasSuspended);
+
         tile.tileUpdateCompleted();
 
         if (tile.getLocation(latitude, longitude) != Tile.Location.HAS_INTERPOLATION_NEIGHBORS) {

src/main/java/org/orekit/rugged/raster/UpdatableTile.java

-/* Copyright 2013-2016 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
@@ -16,26 +16,23 @@
  */
 package org.orekit.rugged.raster;
 
-import org.orekit.rugged.errors.RuggedException;
-
 /** Interface representing one tile of a raster Digital Elevation Model.
  * @author Luc Maisonobe
+ * @author Guylaine Prat
  */
 public interface UpdatableTile {
 
     /** Set the tile global geometry.
-     * @param minLatitude minimum latitude
-     * @param minLongitude minimum longitude
-     * @param latitudeStep step in latitude (size of one raster element)
-     * @param longitudeStep step in longitude (size of one raster element)
+     * @param minLatitude minimum latitude (rad)
+     * @param minLongitude minimum longitude (rad)
+     * @param latitudeStep step in latitude (size of one raster element) (rad)
+     * @param longitudeStep step in longitude (size of one raster element) (rad)
      * @param latitudeRows number of latitude rows
      * @param longitudeColumns number of longitude columns
-     * @exception RuggedException if tile is empty (zero rows or columns)
      */
     void setGeometry(double minLatitude, double minLongitude,
                      double latitudeStep, double longitudeStep,
-                     int latitudeRows, int longitudeColumns)
-        throws RuggedException;
+                     int latitudeRows, int longitudeColumns);
 
     /** Set the elevation for one raster element.
      * <p>
@@ -51,9 +48,7 @@ public interface UpdatableTile {
      * @param latitudeIndex index of latitude (row index)
      * @param longitudeIndex index of longitude (column index)
      * @param elevation elevation (m)
-     * @exception RuggedException if indices are out of bound
      */
-    void setElevation(int latitudeIndex, int longitudeIndex, double elevation)
-        throws RuggedException;
+    void setElevation(int latitudeIndex, int longitudeIndex, double elevation);
 
 }

src/main/java/org/orekit/rugged/raster/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 used to update Digital Elevation
+ * Model tiles to be implemented by the user,
+ * the interface representing a raster tile, as well as a simple model.
+ * It provides also the interface representing a factory for raster tile,
+ * as well as a simple implementation.
+ *
+ * @author Luc Maisonobe
+ * @author Guylaine Prat
+ *
+ */
+package org.orekit.rugged.raster;

src/main/java/org/orekit/rugged/refraction/AtmosphericComputationParameters.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.refraction;
+
+import org.orekit.rugged.errors.RuggedException;
+import org.orekit.rugged.errors.RuggedMessages;
+import org.orekit.rugged.linesensor.LineSensor;
+import org.orekit.rugged.utils.GridCreation;
+
+/**
+ * Atmospheric refraction computation parameters.
+ * Defines for inverse location a set of parameters in order to be able to perform the computation.
+ * @author Guylaine Prat
+ * @since 2.1
+ */
+
+public class AtmosphericComputationParameters {
+
+    /** Margin for definition of the interpolation grid.
+     * To be inside the min line and max line range to avoid problem with inverse location grid computation. */
+    private static final int MARGIN_LINE = 10;
+
+    /** Default value for pixel step. */
+    private static final int DEFAULT_STEP_PIXEL = 100;
+    /** Default value for line step. */
+    private static final int DEFAULT_STEP_LINE = 100;
+
+    /** Default margin for computation of inverse location with atmospheric refraction correction.
+    * @since 3.0
+    */
+    private static final double DEFAULT_INVLOC_MARGIN = 0.8;
+
+    /** Actual values for pixel step in case default are overwritten. */
+    private int pixelStep;
+    /** Actual values for line step in case default are overwritten. */
+    private int lineStep;
+
+    /** Actual values for inverse location margin with atmospheric refraction  in case default are overwritten.
+    * @since 3.0
+    */
+    private double invlocMargin;
+
+    // Definition of grids for sensor (u = along pixel; v = along line)
+    /** Linear grid in pixel. */
+    private double[] uGrid;
+    /** Linear grid in line. */
+    private double[] vGrid;
+    /** Size of uGrid = nbPixelGrid. */
+    private int nbPixelGrid;
+    /** Size of vGrid = nbLineGrid. */
+    private int nbLineGrid;
+
+    // Definition of the associated sensor
+    /** Current min line. */
+    private double minLineSensor = Double.NaN;
+    /** Current max line. */
+    private double maxLineSensor = Double.NaN;
+    /** Current sensor name. */
+    private String sensorName = null;
+
+    /**
+     * Default constructor.
+     */
+    public AtmosphericComputationParameters() {
+        this.pixelStep = DEFAULT_STEP_PIXEL;
+        this.lineStep = DEFAULT_STEP_LINE;
+        this.invlocMargin = DEFAULT_INVLOC_MARGIN;
+    }
+
+    /** Configuration of the interpolation grid. This grid is associated to the given sensor,
+     * with the given min and max lines.
+     * @param sensor line sensor
+     * @param minLine min line defined for the inverse location
+     * @param maxLine max line defined for the inverse location
+     */
+    public void configureCorrectionGrid(final LineSensor sensor, final int minLine, final int maxLine) {
+
+        // Keep information about the sensor and the required search lines.
+        // Needed to test if the grid is initialized with this context.
+        this.minLineSensor = minLine;
+        this.maxLineSensor = maxLine;
+        this.sensorName = sensor.getName();
+
+        // Compute the number of pixels and lines for the grid (round value is sufficient)
+        final int sensorNbPxs = sensor.getNbPixels();
+        this.nbPixelGrid = sensorNbPxs / this.pixelStep;
+
+        // check the validity of the min and max lines
+        if ((maxLine - minLine + 1 - 2 * MARGIN_LINE) < 2 * this.lineStep) {
+            final String info = ": (maxLine - minLine + 1 - 2*" + MARGIN_LINE + ") < 2*" + this.lineStep;
+            throw new RuggedException(RuggedMessages.INVALID_RANGE_FOR_LINES, minLine, maxLine, info);
+        }
+        this.nbLineGrid = (maxLine - minLine + 1 - 2 * MARGIN_LINE) / this.lineStep;
+
+        // Compute the linear grids in pixel (u index) and line (v index)
+        this.uGrid = GridCreation.createLinearGrid(0, sensorNbPxs - 1, this.nbPixelGrid);
+        this.vGrid = GridCreation.createLinearGrid(minLine + MARGIN_LINE, maxLine - MARGIN_LINE, this.nbLineGrid);
+
+    }
+
+    /**
+     * Set the grid steps in pixel and line (used to compute inverse location).
+     * Overwrite the default values, for time optimization if necessary.
+     * @param gridPixelStep grid pixel step for the inverse location computation
+     * @param gridLineStep grid line step for the inverse location computation
+     */
+    public void setGridSteps(final int gridPixelStep, final int gridLineStep) {
+
+        if (gridPixelStep <= 0) {
+            final String reason = " pixelStep <= 0";
+            throw new RuggedException(RuggedMessages.INVALID_STEP, gridPixelStep, reason);
+        }
+        if (gridLineStep <= 0) {
+            final String reason = " lineStep <= 0";
+            throw new RuggedException(RuggedMessages.INVALID_STEP, gridLineStep, reason);
+        }
+        this.pixelStep = gridPixelStep;
+        this.lineStep = gridLineStep;
+    }
+
+    /**
+     * Set the margin for computation of inverse location with atmospheric refraction correction.
+     * Overwrite the default value DEFAULT_INVLOC_MARGIN.
+     * No check is done about this margin. A recommended value is around 1.
+     * @param inverseLocMargin margin in pixel size to compute inverse location with atmospheric refraction correction.
+     * @since 3.0
+     */
+    public void setInverseLocMargin(final double inverseLocMargin) {
+        this.invlocMargin = inverseLocMargin;
+    }
+
+    /**
+     * @return the inverse location margin for computation of inverse location with atmospheric refraction correction.
+    * @since 3.0
+    */
+    public double getInverseLocMargin () {
+        return this.invlocMargin;
+    }
+
+    /**
+    * @return the default inverse location margin for computation of inverse location with atmospheric refraction correction.
+    * @since 3.0
+    */
+    public double getDefaultInverseLocMargin () {
+        return DEFAULT_INVLOC_MARGIN;
+    }
+
+    /**
+     * @return the size of pixel grid
+     */
+    public int getNbPixelGrid() {
+        return nbPixelGrid;
+    }
+
+    /**
+     * @return the size of line grid
+     */
+    public int getNbLineGrid() {
+        return nbLineGrid;
+    }
+
+    /**
+     * @return the pixel grid
+     */
+    public double[] getUgrid() {
+        return uGrid.clone();
+    }
+
+    /**
+     * @return the line grid
+     */
+    public double[] getVgrid() {
+        return vGrid.clone();
+    }
+
+    /**
+     * @return the min line used to compute the current grids
+     */
+    public double getMinLineSensor() {
+        return minLineSensor;
+    }
+
+    /**
+     * @return the max line used to compute the current grids
+     */
+    public double getMaxLineSensor() {
+        return maxLineSensor;
+    }
+
+    /**
+     * @return the sensor name used to compute the current grids
+     */
+    public String getSensorName() {
+        return sensorName;
+    }
+}

src/main/java/org/orekit/rugged/refraction/AtmosphericRefraction.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.refraction;
+
+import org.hipparchus.analysis.interpolation.BilinearInterpolatingFunction;
+import org.hipparchus.geometry.euclidean.threed.Vector3D;
+import org.orekit.rugged.errors.RuggedException;
+import org.orekit.rugged.errors.RuggedMessages;
+import org.orekit.rugged.intersection.IntersectionAlgorithm;
+import org.orekit.rugged.linesensor.LineSensor;
+import org.orekit.rugged.linesensor.SensorPixel;
+import org.orekit.rugged.utils.NormalizedGeodeticPoint;
+
+/**
+ * Base class for atmospheric refraction model.
+ * @author Sergio Esteves
+ * @author Guylaine Prat
+ * @since 2.0
+ */
+public abstract class AtmosphericRefraction {
+
+    /** Flag to tell if we must compute the correction.
+     * By default: computation is set up.
+     * @since 2.1
+     */
+    private boolean mustBeComputed;
+
+    /** The current atmospheric parameters.
+     * @since 2.1
+     */
+    private AtmosphericComputationParameters atmosphericParams;
+
+    /** Bilinear interpolating function for pixel (used by inverse location).
+     * @since 2.1
+    */
+    private BilinearInterpolatingFunction bifPixel;
+
+    /** Bilinear interpolating function of line (used by inverse location).
+     * @since 2.1
+    */
+    private BilinearInterpolatingFunction bifLine;
+
+    /**
+     * Default constructor.
+     */
+    protected AtmosphericRefraction() {
+        // Set up the atmospheric parameters ... with lazy evaluation of the grid (done only if necessary)
+        this.atmosphericParams = new AtmosphericComputationParameters();
+        this.mustBeComputed    = true;
+        this.bifPixel          = null;
+        this.bifLine           = null;
+    }
+
+    /** Apply correction to the intersected point with an atmospheric refraction model.
+     * @param satPos satellite position, in <em>body frame</em>
+     * @param satLos satellite line of sight, in <em>body frame</em>
+     * @param rawIntersection intersection point before refraction correction
+     * @param algorithm intersection algorithm
+     * @return corrected point with the effect of atmospheric refraction
+     * {@link org.orekit.rugged.utils.ExtendedEllipsoid#pointAtAltitude(Vector3D, Vector3D, double)} or see
+     * {@link org.orekit.rugged.intersection.IntersectionAlgorithm#refineIntersection(org.orekit.rugged.utils.ExtendedEllipsoid, Vector3D, Vector3D, NormalizedGeodeticPoint)}
+     */
+    public abstract NormalizedGeodeticPoint applyCorrection(Vector3D satPos, Vector3D satLos, NormalizedGeodeticPoint rawIntersection,
+                                            IntersectionAlgorithm algorithm);
+
+    /** Deactivate computation (needed for the inverse location computation).
+     * @since 2.1
+     */
+    public void deactivateComputation() {
+        this.mustBeComputed = false;
+    }
+
+    /** Reactivate computation (needed for the inverse location computation).
+     * @since 2.1
+     */
+    public void reactivateComputation() {
+        this.mustBeComputed = true;
+    }
+
+    /** Tell if the computation must be performed.
+     * @return true if computation must be performed; false otherwise
+     * @since 2.1
+     */
+    public boolean mustBeComputed() {
+        return mustBeComputed;
+    }
+
+    /** Configuration of the interpolation grid. This grid is associated to the given sensor,
+     * with the given min and max lines.
+     * @param sensor line sensor
+     * @param minLine min line defined for the inverse location
+     * @param maxLine max line defined for the inverse location
+     * @since 2.1
+     */
+    public void configureCorrectionGrid(final LineSensor sensor, final int minLine, final int maxLine) {
+
+        atmosphericParams.configureCorrectionGrid(sensor, minLine, maxLine);
+    }
+
+   /** Check if the current atmospheric parameters are the same as the asked ones.
+    * @param sensorName the asked sensor name
+    * @param minLine the asked min line
+    * @param maxLine the asked max line
+    * @return true if same context; false otherwise
+    * @since 2.1
+    */
+    public Boolean isSameContext(final String sensorName, final int minLine, final int maxLine) {
+
+        return Double.compare(atmosphericParams.getMinLineSensor(), minLine) == 0 &&
+               Double.compare(atmosphericParams.getMaxLineSensor(), maxLine) == 0 &&
+               atmosphericParams.getSensorName().compareTo(sensorName) == 0;
+    }
+
+    /** Get the computation parameters.
+     * @return the AtmosphericComputationParameters
+     * @since 2.1
+     */
+    public AtmosphericComputationParameters getComputationParameters() {
+        return atmosphericParams;
+    }
+
+    /** Set the grid steps in pixel and line (used to compute inverse location).
+     * Overwrite the default values, for time optimization for instance.
+     * @param pixelStep pixel step for the inverse location computation
+     * @param lineStep line step for the inverse location computation
+     * @since 2.1
+     */
+    public void setGridSteps(final int pixelStep, final int lineStep) {
+        atmosphericParams.setGridSteps(pixelStep, lineStep);
+    }
+
+    /**
+     * Set the margin for computation of inverse location with atmospheric refraction correction.
+     * Overwrite the default value DEFAULT_INVLOC_MARGIN.
+     * No check is done about this margin. A recommended value is around 1.
+     * @param inverseLocMargin margin in pixel size to compute inverse location with atmospheric refraction correction.
+     * @since 3.0
+     */
+    public void setInverseLocMargin(final double inverseLocMargin) {
+        atmosphericParams.setInverseLocMargin(inverseLocMargin);
+    }
+
+    /** Compute the correction functions for pixel and lines.
+     * The corrections are computed for pixels and lines, on a regular grid at sensor level.
+     * The corrections are based on the difference on grid nodes (where direct loc is known with atmosphere refraction)
+     * and the sensor pixel found by inverse loc without atmosphere refraction.
+     * The bilinear interpolating functions are then computed for pixel and for line.
+     * Need to be computed only once for a given sensor with the same minLine and maxLine.
+     * @param sensorPixelGridInverseWithout inverse location grid WITHOUT atmospheric refraction
+     * @since 2.1
+     */
+    public void computeGridCorrectionFunctions(final SensorPixel[][] sensorPixelGridInverseWithout) {
+
+        final int nbPixelGrid = atmosphericParams.getNbPixelGrid();
+        final int nbLineGrid = atmosphericParams.getNbLineGrid();
+        final double[] pixelGrid = atmosphericParams.getUgrid();
+        final double[] lineGrid = atmosphericParams.getVgrid();
+
+        // Initialize the needed diff functions
+        final double[][] gridDiffPixel = new double[nbPixelGrid][nbLineGrid];
+        final double[][] gridDiffLine = new double[nbPixelGrid][nbLineGrid];
+
+        // Compute the difference between grids nodes WITH - without atmosphere
+        for (int lineIndex = 0; lineIndex < nbLineGrid; lineIndex++) {
+            for (int pixelIndex = 0; pixelIndex < nbPixelGrid; pixelIndex++) {
+
+                if (sensorPixelGridInverseWithout[pixelIndex][lineIndex] != null) {
+                    final double diffLine = lineGrid[lineIndex] - sensorPixelGridInverseWithout[pixelIndex][lineIndex].getLineNumber();
+                    final double diffPixel = pixelGrid[pixelIndex] - sensorPixelGridInverseWithout[pixelIndex][lineIndex].getPixelNumber();
+                    gridDiffPixel[pixelIndex][lineIndex] = diffPixel;
+                    gridDiffLine[pixelIndex][lineIndex] = diffLine;
+
+                } else {
+                    // Impossible to find the sensor pixel in the given range lines
+                    throw new RuggedException(RuggedMessages.SENSOR_PIXEL_NOT_FOUND_IN_RANGE_LINES,
+                                              atmosphericParams.getMinLineSensor(), atmosphericParams.getMaxLineSensor());
+                }
+            }
+        }
+        // Definition of the interpolating function for pixel and for line
+        this.bifPixel = new BilinearInterpolatingFunction(pixelGrid, lineGrid, gridDiffPixel);
+        this.bifLine = new BilinearInterpolatingFunction(pixelGrid, lineGrid, gridDiffLine);
+    }
+
+    /**
+     * @return the bilinear interpolating function for pixel correction
+     */
+    public BilinearInterpolatingFunction getBifPixel() {
+        return bifPixel;
+    }
+
+    /**
+     * @return the bilinear interpolating function for line correction
+     */
+    public BilinearInterpolatingFunction getBifLine() {
+        return bifLine;
+    }
+}

src/main/java/org/orekit/rugged/refraction/ConstantRefractionLayer.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.refraction;
+
+/**
+ * Class that represents a constant refraction layer to be used with {@link MultiLayerModel}.
+ * @author Sergio Esteves
+ * @author Guylaine Prat
+ * @since 2.0
+ */
+public class ConstantRefractionLayer {
+
+    /** lowest altitude of this layer (m). */
+    private final Double lowestAltitude;
+
+    /** refractive index of this layer. */
+    private final double refractiveIndex;
+
+    /** Simple constructor.
+     * @param lowestAltitude lowest altitude of the layer (m)
+     * @param refractiveIndex refractive index of the layer
+     */
+    public ConstantRefractionLayer(final double lowestAltitude, final double refractiveIndex) {
+        this.lowestAltitude = lowestAltitude;
+        this.refractiveIndex = refractiveIndex;
+    }
+
+    /**
+     * @return the lowest altitude of the layer (m)
+     */
+    public double getLowestAltitude() {
+        return lowestAltitude;
+    }
+
+    /**
+     * @return the refractive index of the layer
+     */
+    public double getRefractiveIndex() {
+        return refractiveIndex;
+    }
+}

src/main/java/org/orekit/rugged/refraction/MultiLayerModel.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.refraction;
+
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.List;
+
+import org.hipparchus.geometry.euclidean.threed.Vector3D;
+import org.hipparchus.util.FastMath;
+import org.orekit.bodies.GeodeticPoint;
+import org.orekit.rugged.errors.RuggedException;
+import org.orekit.rugged.errors.RuggedMessages;
+import org.orekit.rugged.intersection.IntersectionAlgorithm;
+import org.orekit.rugged.utils.ExtendedEllipsoid;
+import org.orekit.rugged.utils.NormalizedGeodeticPoint;
+
+/**
+ * Atmospheric refraction model based on multiple layers with associated refractive index.
+ * @author Sergio Esteves
+ * @author Luc Maisonobe
+ * @author Guylaine Prat
+ * @since 2.0
+ */
+public class MultiLayerModel extends AtmosphericRefraction {
+
+    /** Observed body ellipsoid. */
+    private final ExtendedEllipsoid ellipsoid;
+
+    /** Constant refraction layers. */
+    private final List<ConstantRefractionLayer> refractionLayers;
+
+    /** Atmosphere lowest altitude (m). */
+    private final double atmosphereLowestAltitude;
+
+    /** Simple constructor.
+     * <p>
+     * This model uses a built-in set of layers.
+     * </p>
+     * @param ellipsoid the ellipsoid to be used.
+     */
+    public MultiLayerModel(final ExtendedEllipsoid ellipsoid) {
+
+        super();
+
+        this.ellipsoid = ellipsoid;
+
+        this.refractionLayers = new ArrayList<>(15);
+        this.refractionLayers.add(new ConstantRefractionLayer(100000.00, 1.000000));
+        this.refractionLayers.add(new ConstantRefractionLayer( 50000.00, 1.000000));
+        this.refractionLayers.add(new ConstantRefractionLayer( 40000.00, 1.000001));
+        this.refractionLayers.add(new ConstantRefractionLayer( 30000.00, 1.000004));
+        this.refractionLayers.add(new ConstantRefractionLayer( 23000.00, 1.000012));
+        this.refractionLayers.add(new ConstantRefractionLayer( 18000.00, 1.000028));
+        this.refractionLayers.add(new ConstantRefractionLayer( 14000.00, 1.000052));
+        this.refractionLayers.add(new ConstantRefractionLayer( 11000.00, 1.000083));
+        this.refractionLayers.add(new ConstantRefractionLayer(  9000.00, 1.000106));
+        this.refractionLayers.add(new ConstantRefractionLayer(  7000.00, 1.000134));
+        this.refractionLayers.add(new ConstantRefractionLayer(  5000.00, 1.000167));
+        this.refractionLayers.add(new ConstantRefractionLayer(  3000.00, 1.000206));
+        this.refractionLayers.add(new ConstantRefractionLayer(  1000.00, 1.000252));
+        this.refractionLayers.add(new ConstantRefractionLayer(     0.00, 1.000278));
+        this.refractionLayers.add(new ConstantRefractionLayer( -1000.00, 1.000306));
+
+        // get the lowest altitude of the atmospheric model
+        this.atmosphereLowestAltitude = refractionLayers.get(refractionLayers.size() - 1).getLowestAltitude();
+    }
+
+    /** Simple constructor.
+     * @param ellipsoid the ellipsoid to be used.
+     * @param refractionLayers the refraction layers to be used with this model (layers can be in any order).
+     */
+    public MultiLayerModel(final ExtendedEllipsoid ellipsoid, final List<ConstantRefractionLayer> refractionLayers) {
+
+        super();
+        // at this stage no optimization is set up: no optimization grid is defined
+
+        this.ellipsoid = ellipsoid;
+        this.refractionLayers = new ArrayList<>(refractionLayers);
+
+        // sort the layers from the highest (index = 0) to the lowest (index = size - 1)
+        Collections.sort(this.refractionLayers,
+            (l1, l2) -> Double.compare(l2.getLowestAltitude(), l1.getLowestAltitude()));
+
+        // get the lowest altitude of the model
+        atmosphereLowestAltitude = this.refractionLayers.get(this.refractionLayers.size() - 1).getLowestAltitude();
+    }
+
+    /** Compute the (position, LOS) of the intersection with the lowest atmospheric layer.
+     * @param satPos satellite position, in body frame
+     * @param satLos satellite line of sight, in body frame
+     * @param rawIntersection intersection point without refraction correction
+     * @return the intersection position and LOS with the lowest atmospheric layer
+     * @since 2.1
+     */
+    private IntersectionLOS computeToLowestAtmosphereLayer(
+                            final Vector3D satPos, final Vector3D satLos,
+                            final NormalizedGeodeticPoint rawIntersection) {
+
+        if (rawIntersection.getAltitude() < atmosphereLowestAltitude) {
+            throw new RuggedException(RuggedMessages.NO_LAYER_DATA, rawIntersection.getAltitude(),
+                                      atmosphereLowestAltitude);
+        }
+
+        Vector3D pos = satPos;
+        Vector3D los = satLos.normalize();
+
+        // Compute the intersection point with the lowest layer of atmosphere
+        double n1 = -1;
+        GeodeticPoint gp = ellipsoid.transform(satPos, ellipsoid.getBodyFrame(), null);
+
+        // Perform the exact computation (no optimization)
+        // TBN: the refractionLayers is ordered from the highest to the lowest
+        for (ConstantRefractionLayer refractionLayer : refractionLayers) {
+
+            if (refractionLayer.getLowestAltitude() > gp.getAltitude()) {
+                continue;
+            }
+
+            final double n2 = refractionLayer.getRefractiveIndex();
+
+            if (n1 > 0) {
+
+                // when we get here, we have already performed one iteration in the loop
+                // so gp is the los intersection with the layers interface (it was a
+                // point on ground at loop initialization, but is overridden at each iteration)
+
+                // get new los by applying Snell's law at atmosphere layers interfaces
+                // we avoid computing sequences of inverse-trigo/trigo/inverse-trigo functions
+                // we just use linear algebra and square roots, it is faster and more accurate
+
+                // at interface crossing, the interface normal is z, the local zenith direction
+                // the ray direction (i.e. los) is u in the upper layer and v in the lower layer
+                // v is in the (u, zenith) plane, so we can say
+                //  (1) v = α u + β z
+                // with α>0 as u and v are roughly in the same direction as the ray is slightly bent
+
+                // let θ₁ be the los incidence angle at interface crossing
+                // θ₁ = π - angle(u, zenith) is between 0 and π/2 for a downwards observation
+                // let θ₂ be the exit angle at interface crossing
+                // from Snell's law, we have n₁ sin θ₁ = n₂ sin θ₂ and θ₂ is also between 0 and π/2
+                // we have:
+                //   (2) u·z = -cos θ₁
+                //   (3) v·z = -cos θ₂
+                // combining equations (1), (2) and (3) and remembering z·z = 1 as z is normalized , we get
+                //   (4) β = α cos θ₁ - cos θ₂
+                // with all the expressions above, we can rewrite the fact v is normalized:
+                //       1 = v·v
+                //         = α² u·u + 2αβ u·z + β² z·z
+                //         = α² - 2αβ cos θ₁ + β²
+                //         = α² - 2α² cos² θ₁ + 2 α cos θ₁ cos θ₂ + α² cos² θ₁ - 2 α cos θ₁ cos θ₂ + cos² θ₂
+                //         = α²(1 - cos² θ₁) + cos² θ₂
+                // hence α² = (1 - cos² θ₂)/(1 - cos² θ₁)
+                //          = sin² θ₂ / sin² θ₁
+                // as α is positive, and both θ₁ and θ₂ are between 0 and π/2, we finally get
+                //       α  = sin θ₂ / sin θ₁
+                //   (5) α  = n₁/n₂
+                // the α coefficient is independent from the incidence angle,
+                // it depends only on the ratio of refractive indices!
+                //
+                // back to equation (4) and using again the fact θ₂ is between 0 and π/2, we can now write
+                //       β = α cos θ₁ - cos θ₂
+                //         = n₁/n₂ cos θ₁ - cos θ₂
+                //         = n₁/n₂ cos θ₁ - √(1 - sin² θ₂)
+                //         = n₁/n₂ cos θ₁ - √(1 - (n₁/n₂)² sin² θ₁)
+                //         = n₁/n₂ cos θ₁ - √(1 - (n₁/n₂)² (1 - cos² θ₁))
+                //         = n₁/n₂ cos θ₁ - √(1 - (n₁/n₂)² + (n₁/n₂)² cos² θ₁)
+                //   (6) β = -k - √(k² - ζ)
+                // where ζ = (n₁/n₂)² - 1 and k = n₁/n₂ u·z, which is negative, and close to -1 for
+                // nadir observations. As we expect atmosphere models to have small transitions between
+                // layers, we have to handle accurately the case where n₁/n₂ ≈ 1 so ζ ≈ 0. In this case,
+                // a cancellation occurs inside the square root: √(k² - ζ) ≈ √k² ≈ -k (because k is negative).
+                // So β ≈ -k + k ≈ 0 and another cancellation occurs, outside of the square root.
+                // This means that despite equation (6) is mathematically correct, it is prone to numerical
+                // errors when consecutive layers have close refractive indices. A better equivalent
+                // expression is needed. The fact β is close to 0 in this case was expected because
+                // equation (1) reads v = α u + β z, and α = n₁/n₂, so when n₁/n₂ ≈ 1, we have
+                // α ≈ 1 and β ≈ 0, so v ≈ u: when two layers have similar refractive indices, the
+                // propagation direction is almost unchanged.
+                //
+                // The first step for the accurate computation of β is to compute ζ = (n₁/n₂)² - 1
+                // accurately and avoid a cancellation just after a division (which is the least accurate
+                // of the four operations) and a squaring. We will simply use:
+                //   ζ = (n₁/n₂)² - 1
+                //     = (n₁ - n₂) (n₁ + n₂) / n₂²
+                // The cancellation is still there, but it occurs in the subtraction n₁ - n₂, which are
+                // the most accurate values we can get.
+                // The second step for the accurate computation of β is to rewrite equation (6)
+                // by both multiplying and dividing by the dual factor -k + √(k² - ζ):
+                //     β = -k - √(k² - ζ)
+                //       = (-k - √(k² - ζ)) * (-k + √(k² - ζ)) / (-k + √(k² - ζ))
+                //       = (k² - (k² - ζ)) / (-k + √(k² - ζ))
+                // (7) β = ζ / (-k + √(k² - ζ))
+                // expression (7) is more stable numerically than expression (6), because when ζ ≈ 0
+                // its denominator is about -2k, there are no cancellation anymore after the square root.
+                // β is computed with the same accuracy as ζ
+                final double alpha = n1 / n2;
+                final double k     = alpha * Vector3D.dotProduct(los, gp.getZenith());
+                final double zeta  = (n1 - n2) * (n1 + n2) / (n2 * n2);
+                final double beta  = zeta / (FastMath.sqrt(k * k - zeta) - k);
+                los = new Vector3D(alpha, los, beta, gp.getZenith());
+            }
+
+            // In case the altitude of the intersection without atmospheric refraction
+            // is above the lowest altitude of the atmosphere: stop the search
+            if (rawIntersection.getAltitude() > refractionLayer.getLowestAltitude()) {
+                break;
+            }
+
+            // Get for the intersection point: the position and the LOS
+            pos = ellipsoid.pointAtAltitude(pos, los, refractionLayer.getLowestAltitude());
+            gp  = ellipsoid.transform(pos, ellipsoid.getBodyFrame(), null);
+
+            n1 = n2;
+        }
+        return new IntersectionLOS(pos, los);
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public NormalizedGeodeticPoint applyCorrection(final Vector3D satPos, final Vector3D satLos,
+                                                   final NormalizedGeodeticPoint rawIntersection,
+                                                   final IntersectionAlgorithm algorithm) {
+
+        final IntersectionLOS intersectionLOS = computeToLowestAtmosphereLayer(satPos, satLos, rawIntersection);
+        final Vector3D pos = intersectionLOS.getIntersectionPos();
+        final Vector3D los = intersectionLOS.getIntersectionLos();
+
+        // at this stage the pos belongs to the lowest atmospheric layer.
+        // We can compute the intersection of line of sight (los) with Digital Elevation Model
+        // as usual (without atmospheric refraction).
+        return algorithm.refineIntersection(ellipsoid, pos, los, rawIntersection);
+    }
+
+} // end of class MultiLayerModel
+
+/** Container for the (position, LOS) of the intersection with the lowest atmospheric layer.
+ * @author Guylaine Prat
+ * @since 2.1
+ */
+class IntersectionLOS {
+
+    /** Position of the intersection with the lowest atmospheric layer. */
+    private Vector3D intersectionPos;
+    /** LOS of the intersection with the lowest atmospheric layer. */
+    private Vector3D intersectionLos;
+
+    /** Default constructor.
+     * @param intersectionPos position of the intersection
+     * @param intersectionLos los of the intersection
+     */
+    IntersectionLOS(final Vector3D intersectionPos, final Vector3D intersectionLos) {
+        this.intersectionPos = intersectionPos;
+        this.intersectionLos = intersectionLos;
+    }
+    public Vector3D getIntersectionPos() {
+        return intersectionPos;
+    }
+
+    public Vector3D getIntersectionLos() {
+        return intersectionLos;
+    }
+}

src/main/java/org/orekit/rugged/refraction/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 atmospheric refraction model, as well as
+ * some classical models.
+ *
+ * @author Luc Maisonobe
+ * @author Guylaine Prat
+ * @author Sergio Esteves
+ *
+ */
+package org.orekit.rugged.refraction;

src/main/java/org/orekit/rugged/utils/AbsoluteDateArrayHandling.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.utils;
+
+import org.hipparchus.exception.LocalizedCoreFormats;
+import org.orekit.errors.OrekitException;
+import org.orekit.time.AbsoluteDate;
+
+/** AbsoluteDateArrayHandling consist of additions to AbsoluteDate to handle arrays.
+ * @author Melina Vanel
+ */
+public class AbsoluteDateArrayHandling {
+
+    /** Dates array on which we want to apply time shift or compute duration. */
+    private final AbsoluteDate[] dates;
+
+    /** Simple constructor.
+     * @param dates is an array of absolute dates on which we want to apply time shift or
+     * compute duration
+     */
+    public AbsoluteDateArrayHandling(final AbsoluteDate[] dates) {
+        this.dates = dates.clone();
+    }
+
+    /** Get instance dates array.
+     * @return dates array
+     */
+    public AbsoluteDate[] getDates() {
+        return this.dates.clone();
+    }
+
+    /** Get time-shifted dates for several dates or several time shifts.
+     * If instance dates = [date1, date2, ..., daten] and argument
+     * dts = [dts1, dts2, ..., dtsn] then this function will return a matrix
+     * [[date1 shiftedby dts1, date1 shiftedBy dts2, ..., date1 shiftedBy dtsn],
+     * [date2 shiftedby dts1, date2 shiftedBy dts2, ..., date2 shiftedBy dtsn],
+     * [...]
+     * [daten shiftedby dts1, daten shiftedBy dts2, ..., date1 shiftedBy dtsn]].
+     * If ones want to apply only 1 time shift corresponding to 1 date see
+     * {@link #shiftedBy(double[])}.
+     * @param dts time shifts array in seconds we want to apply to dates
+     * @return a matrix of new dates, shifted with respect to wanted time
+     * shifts. If instance dates = [date1, date2, ..., daten] each line
+     * correspond to one date (for example date1 shiftedBy all timeshifts
+     * (building the different columns))
+     */
+    public AbsoluteDate[][] multipleShiftedBy(final double[] dts) {
+
+        final AbsoluteDate[][] datesShifted = new AbsoluteDate[dates.length][dts.length];
+        int index_dates = 0;
+
+        for (AbsoluteDate date: this.dates) {
+            final AbsoluteDate[] dateShifted = new AbsoluteDate[dts.length];
+            int index_dts = 0;
+            for (double dt: dts) {
+                dateShifted[index_dts] = date.shiftedBy(dt);
+                index_dts += 1;
+            }
+            datesShifted[index_dates] = dateShifted;
+            index_dates += 1;
+
+        }
+        return (AbsoluteDate[][]) datesShifted;
+    }
+
+    /** Get time-shifted dates for several dates and corresponding time shifts.
+     * If instance dates = [date1, date2, ..., daten] and argument
+     * dts = [dts1, dts2, ..., dtsn] then this function will return
+     * [date1 shiftedby dts1, date2 shiftedBy dts2, ..., daten shiftedBy dtsn]. If
+     * several time shift want to be applied on each date see
+     * {@link #multipleShiftedBy(double[])}.
+     * @param dts time shifts array in seconds we want to apply to corresponding dates.
+     * Warning, must be same length as dates.
+     * @return an 1D array of new dates, shifted with respect to wanted corresponding time
+     * shifts.
+     */
+    public AbsoluteDate[] shiftedBy(final double[] dts) {
+
+        // Check same dimensions
+        if (dates.length != dts.length) {
+            throw new OrekitException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
+                                      dates.length, dts.length);
+        }
+
+        final AbsoluteDate[] datesShifted = new AbsoluteDate[dates.length];
+        int index_dates = 0;
+
+        for (AbsoluteDate date: this.dates) {
+            datesShifted[index_dates] = date.shiftedBy(dts[index_dates]);
+            index_dates += 1;
+
+        }
+        return datesShifted;
+    }
+
+    /** Get array with durations between instances dates and given dates
+     * If instance dates = [date1, date2, ..., daten] and argument
+     * datesForDuration = [d1, d2, ..., dn] then this function will return a matrix
+     * [[date1 durationFrom d1, date1 durationFrom d2, ..., date1 durationFrom dn],
+     * [date2 durationFrom d1, date2 durationFrom d2, ..., date2 durationFrom dn],
+     * [...]
+     * [daten durationFrom d1, daten durationFrom d2, ..., date1 durationFrom dn]].
+     * If ones want to compute duration from only 1 date corresponding to 1 instance date see
+     * {@link #durationFrom(AbsoluteDate[])}.
+     * @param datesForDuration dates for which we want to compute the duration form instances dates
+     * @return a matrix of double representing durations from instance dates
+     * If instance dates = [date1, date2, ..., daten] each line
+     * correspond to one date (for example date1 duration from all given dates in arguments
+     * (building the different columns))
+     */
+    public double[][] multipleDurationFrom(final AbsoluteDate[] datesForDuration) {
+
+        final double[][] durationsFromDates = new double[dates.length][datesForDuration.length];
+        int index_dates = 0;
+
+        for (AbsoluteDate date: this.dates) {
+            final double[] durationFromDate = new double[datesForDuration.length];
+            int index_datesForDuration = 0;
+            for (AbsoluteDate dateForDuration: datesForDuration) {
+                durationFromDate[index_datesForDuration] = date.durationFrom(dateForDuration);
+                index_datesForDuration += 1;
+            }
+            durationsFromDates[index_dates] = durationFromDate;
+            index_dates += 1;
+
+        }
+        return (double[][]) durationsFromDates;
+    }
+
+    /** Get array with durations between instances dates and corresponding given dates
+     * If instance dates = [date1, date2, ..., daten] and argument
+     * datesForDuration = [d1, d2, ..., dn] then this function will return
+     * [date1 durationFrom d1, date2 durationFrom d2, ..., daten durationFrom dn]. If
+     * duration from from all arguments dates wants to be compute on each date see
+     * {@link #multipleDurationFrom(AbsoluteDate[])}.
+     * @param datesForDuration dates for which we want to compute the duration form instances dates.
+     * Warning must have same length as instance dates.
+     * @return a array of double representing durations between instance dates and corresponding
+     * argument dates
+     */
+    public double[] durationFrom(final AbsoluteDate[] datesForDuration) {
+
+        // Check same dimensions
+        if (dates.length != datesForDuration.length) {
+            throw new OrekitException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
+                                      dates.length, datesForDuration.length);
+        }
+
+        final double[] durationsFromDates = new double[dates.length];
+        int index_dates = 0;
+
+        for (AbsoluteDate date: this.dates) {
+            durationsFromDates[index_dates] = date.durationFrom(datesForDuration[index_dates]);
+            index_dates += 1;
+
+        }
+        return durationsFromDates;
+    }
+
+}

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

-/* Copyright 2013-2016 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
@@ -16,8 +16,6 @@
  */
 package org.orekit.rugged.utils;
 
-import java.util.List;
-
 import org.hipparchus.analysis.differentiation.DerivativeStructure;
 import org.orekit.utils.ParameterDriver;
 
@@ -27,29 +25,8 @@ import org.orekit.utils.ParameterDriver;
  * </p>
  * @author Luc Maisonobe
  * @since 2.0
+ * @deprecated as of 2.2, replaced by {@link DerivativeGenerator}
  */
-public interface DSGenerator {
-
-    /** Get the parameters selected for estimation.
-     * @return parameters selected for estimation
-     */
-    List<ParameterDriver> getSelected();
-
-    /** Generate a constant {@link DerivativeStructure}.
-     * @param value value of the constant
-     * @return constant {@link DerivativeStructure}
-     */
-    DerivativeStructure constant(double value);
-
-    /** Generate a {@link DerivativeStructure} representing the
-     * parameter driver either as a canonical variable or a constant.
-     * <p>
-     * The instance created is a variable only if the parameter
-     * has been selected for estimation, otherwise it is a constant.
-     * </p>
-     * @param driver driver for the variable
-     * @return variable {@link DerivativeStructure}
-     */
-    DerivativeStructure variable(ParameterDriver driver);
-
+public interface DSGenerator extends DerivativeGenerator<DerivativeStructure> {
+    // nothing specialized here
 }

src/main/java/org/orekit/rugged/utils/DerivativeGenerator.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.utils;
+
+import java.util.List;
+
+import org.hipparchus.Field;
+import org.hipparchus.analysis.differentiation.Derivative;
+import org.orekit.utils.ParameterDriver;
+
+/** Generator for {@link Derivative} instances from {@link ParameterDriver}.
+ * <p>
+ * Note that this interface is for Rugged library internal use only.
+ * </p>
+ * @author Luc Maisonobe
+ * @since 2.0
+ */
+public interface DerivativeGenerator<T extends Derivative<T>> {
+
+    /** Get the parameters selected for estimation.
+     * @return parameters selected for estimation
+     */
+    List<ParameterDriver> getSelected();
+
+    /** Generate a constant {@link Derivative}.
+     * @param value value of the constant
+     * @return constant {@link Derivative}
+     */
+    T constant(double value);
+
+    /** Generate a {@link Derivative} representing the
+     * parameter driver either as a canonical variable or a constant.
+     * <p>
+     * The instance created is a variable only if the parameter
+     * has been selected for estimation, otherwise it is a constant.
+     * </p>
+     * @param driver driver for the variable
+     * @return variable {@link Derivative}
+     */
+    T variable(ParameterDriver driver);
+
+    /** Get the {@link Field} to which the generated derivatives belongs.
+     * @return {@link Field} to which the generated derivatives belongs
+     * @since 2.2
+     */
+    default Field<T> getField() {
+        return constant(0).getField();
+    }
+
+}

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

-/* Copyright 2013-2016 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
@@ -22,7 +22,6 @@ import org.hipparchus.util.FastMath;
 import org.hipparchus.util.MathArrays;
 import org.orekit.bodies.GeodeticPoint;
 import org.orekit.bodies.OneAxisEllipsoid;
-import org.orekit.errors.OrekitException;
 import org.orekit.frames.Frame;
 import org.orekit.rugged.errors.DumpManager;
 import org.orekit.rugged.errors.RuggedException;
@@ -47,7 +46,7 @@ public class ExtendedEllipsoid extends OneAxisEllipsoid {
     private final double b2;
 
     /** Simple constructor.
-     * @param ae equatorial radius
+     * @param ae equatorial radius (m)
      * @param f the flattening (f = (a-b)/a)
      * @param bodyFrame body frame related to body shape
      * @see org.orekit.frames.FramesFactory#getITRF(org.orekit.utils.IERSConventions, boolean)
@@ -68,25 +67,22 @@ public class ExtendedEllipsoid extends OneAxisEllipsoid {
 
     /** {@inheritDoc} */
     @Override
-    public GeodeticPoint transform(final Vector3D point, final Frame frame, final AbsoluteDate date)
-        throws OrekitException {
+    public GeodeticPoint transform(final Vector3D point, final Frame frame, final AbsoluteDate date) {
         DumpManager.dumpEllipsoid(this);
         return super.transform(point, frame, date);
     }
 
     /** Get point at some latitude along a pixel line of sight.
-     * @param position cell position (in body frame)
+     * @param position cell position (in body frame) (m)
      * @param los pixel line-of-sight, not necessarily normalized (in body frame)
-     * @param latitude latitude with respect to ellipsoid
+     * @param latitude latitude with respect to ellipsoid (rad)
      * @param closeReference reference point used to select the closest solution
      * when there are two points at the desired latitude along the line, it should
-     * be close to los surface intersection
-     * @return point at latitude
-     * @exception RuggedException if no such point exists
+     * be close to los surface intersection (m)
+     * @return point at latitude (m)
      */
     public Vector3D pointAtLatitude(final Vector3D position, final Vector3D los,
-                                    final double latitude, final Vector3D closeReference)
-        throws RuggedException {
+                                    final double latitude, final Vector3D closeReference) {
 
         DumpManager.dumpEllipsoid(this);
 
@@ -161,14 +157,12 @@ public class ExtendedEllipsoid extends OneAxisEllipsoid {
     }
 
     /** Get point at some longitude along a pixel line of sight.
-     * @param position cell position (in body frame)
+     * @param position cell position (in body frame) (m)
      * @param los pixel line-of-sight, not necessarily normalized (in body frame)
-     * @param longitude longitude with respect to ellipsoid
-     * @return point at longitude
-     * @exception RuggedException if no such point exists
+     * @param longitude longitude with respect to ellipsoid (rad)
+     * @return point at longitude (m)
      */
-    public Vector3D pointAtLongitude(final Vector3D position, final Vector3D los, final double longitude)
-        throws RuggedException {
+    public Vector3D pointAtLongitude(final Vector3D position, final Vector3D los, final double longitude) {
 
         DumpManager.dumpEllipsoid(this);
 
@@ -186,85 +180,72 @@ public class ExtendedEllipsoid extends OneAxisEllipsoid {
     }
 
     /** Get point on ground along a pixel line of sight.
-     * @param position cell position (in body frame)
+     * @param position cell position (in body frame) (m)
      * @param los pixel line-of-sight, not necessarily normalized (in body frame)
      * @param centralLongitude reference longitude lc such that the point longitude will
-     * be normalized between lc-π and lc+π
+     * be normalized between lc-π and lc+π (rad)
      * @return point on ground
-     * @exception RuggedException if no such point exists (typically line-of-sight missing body)
      */
     public NormalizedGeodeticPoint pointOnGround(final Vector3D position, final Vector3D los,
-                                                 final double centralLongitude)
-        throws RuggedException {
-        try {
-            DumpManager.dumpEllipsoid(this);
-            final GeodeticPoint gp =
-                    getIntersectionPoint(new Line(position, new Vector3D(1, position, 1e6, los), 1.0e-12),
-                                         position, getBodyFrame(), null);
-            if (gp == null) {
-                throw new RuggedException(RuggedMessages.LINE_OF_SIGHT_DOES_NOT_REACH_GROUND);
-            }
-            return new NormalizedGeodeticPoint(gp.getLatitude(), gp.getLongitude(), gp.getAltitude(),
-                                               centralLongitude);
-        } catch (OrekitException oe) {
-            throw new RuggedException(oe, oe.getSpecifier(), oe.getParts());
+                                                 final double centralLongitude) {
+
+        DumpManager.dumpEllipsoid(this);
+        final GeodeticPoint gp =
+                getIntersectionPoint(new Line(position, new Vector3D(1, position, 1e6, los), 1.0e-12),
+                        position, getBodyFrame(), null);
+        if (gp == null) {
+            throw new RuggedException(RuggedMessages.LINE_OF_SIGHT_DOES_NOT_REACH_GROUND);
         }
+        return new NormalizedGeodeticPoint(gp.getLatitude(), gp.getLongitude(), gp.getAltitude(),
+                centralLongitude);
     }
 
     /** Get point at some altitude along a pixel line of sight.
-     * @param position cell position (in body frame)
+     * @param position cell position (in body frame) (m)
      * @param los pixel line-of-sight, not necessarily normalized (in body frame)
-     * @param altitude altitude with respect to ellipsoid
-     * @return point at altitude
-     * @exception RuggedException if no such point exists (typically too negative altitude)
+     * @param altitude altitude with respect to ellipsoid (m)
+     * @return point at altitude (m)
      */
-    public Vector3D pointAtAltitude(final Vector3D position, final Vector3D los, final double altitude)
-        throws RuggedException {
-        try {
-
-            DumpManager.dumpEllipsoid(this);
-
-            // point on line closest to origin
-            final double   los2   = los.getNormSq();
-            final double   dot    = Vector3D.dotProduct(position, los);
-            final double   k0     = -dot / los2;
-            final Vector3D close0 = new Vector3D(1, position, k0, los);
-
-            // very rough guess: if body is spherical, the desired point on line
-            // is at distance ae + altitude from origin
-            final double r        = getEquatorialRadius() + altitude;
-            final double delta2   = r * r - close0.getNormSq();
-            if (delta2 < 0) {
-                throw new RuggedException(RuggedMessages.LINE_OF_SIGHT_NEVER_CROSSES_ALTITUDE, altitude);
-            }
-            final double deltaK   = FastMath.sqrt(delta2 / los2);
-            final double k1       = k0 + deltaK;
-            final double k2       = k0 - deltaK;
-            double k              = (FastMath.abs(k1) <= FastMath.abs(k2)) ? k1 : k2;
-
-            // this loop generally converges in 3 iterations
-            for (int i = 0; i < 100; ++i) {
-
-                final Vector3D      point   = new Vector3D(1, position, k, los);
-                final GeodeticPoint gpK     = transform(point, getBodyFrame(), null);
-                final double        deltaH  = altitude - gpK.getAltitude();
-                if (FastMath.abs(deltaH) <= ALTITUDE_CONVERGENCE) {
-                    return point;
-                }
+    public Vector3D pointAtAltitude(final Vector3D position, final Vector3D los, final double altitude) {
 
-                // improve the offset using linear ratio between
-                // altitude variation and displacement along line-of-sight
-                k += deltaH / Vector3D.dotProduct(gpK.getZenith(), los);
+        DumpManager.dumpEllipsoid(this);
 
+        // point on line closest to origin
+        final double   los2   = los.getNormSq();
+        final double   dot    = Vector3D.dotProduct(position, los);
+        final double   k0     = -dot / los2;
+        final Vector3D close0 = new Vector3D(1, position, k0, los);
+
+        // very rough guess: if body is spherical, the desired point on line
+        // is at distance ae + altitude from origin
+        final double r        = getEquatorialRadius() + altitude;
+        final double delta2   = r * r - close0.getNormSq();
+        if (delta2 < 0) {
+            throw new RuggedException(RuggedMessages.LINE_OF_SIGHT_NEVER_CROSSES_ALTITUDE, altitude);
+        }
+        final double deltaK   = FastMath.sqrt(delta2 / los2);
+        final double k1       = k0 + deltaK;
+        final double k2       = k0 - deltaK;
+        double k              = (FastMath.abs(k1) <= FastMath.abs(k2)) ? k1 : k2;
+
+        // this loop generally converges in 3 iterations
+        for (int i = 0; i < 100; ++i) {
+
+            final Vector3D      point   = new Vector3D(1, position, k, los);
+            final GeodeticPoint gpK     = transform(point, getBodyFrame(), null);
+            final double        deltaH  = altitude - gpK.getAltitude();
+            if (FastMath.abs(deltaH) <= ALTITUDE_CONVERGENCE) {
+                return point;
             }
 
-            // this should never happen
-            throw new RuggedException(RuggedMessages.LINE_OF_SIGHT_NEVER_CROSSES_ALTITUDE, altitude);
+            // improve the offset using linear ratio between
+            // altitude variation and displacement along line-of-sight
+            k += deltaH / Vector3D.dotProduct(gpK.getZenith(), los);
 
-        } catch (OrekitException oe) {
-            // this should never happen
-            throw new RuggedException(oe, oe.getSpecifier(), oe.getParts());
         }
+
+        // this should never happen
+        throw new RuggedException(RuggedMessages.LINE_OF_SIGHT_NEVER_CROSSES_ALTITUDE, altitude);
     }
 
     /** Convert a line-of-sight from Cartesian to topocentric.
@@ -302,36 +283,27 @@ public class ExtendedEllipsoid extends OneAxisEllipsoid {
      * with respect to the primary point (in body frame and Cartesian coordinates)
      * @return line-of-sight in topocentric frame (East, North, Zenith) of the point,
      * scaled to match radians in the horizontal plane and meters along the vertical axis
-     * @exception RuggedException if points cannot be converted to geodetic coordinates
      */
-    public Vector3D convertLos(final Vector3D primary, final Vector3D secondary)
-        throws RuggedException {
-        try {
-
-            // switch to geodetic coordinates using primary point as reference
-            final GeodeticPoint point = transform(primary, getBodyFrame(), null);
-            final Vector3D      los   = secondary.subtract(primary);
+    public Vector3D convertLos(final Vector3D primary, final Vector3D secondary) {
 
-            // convert line of sight
-            return convertLos(point, los);
+        // switch to geodetic coordinates using primary point as reference
+        final GeodeticPoint point = transform(primary, getBodyFrame(), null);
+        final Vector3D      los   = secondary.subtract(primary);
 
-        } catch (OrekitException oe) {
-            throw new RuggedException(oe, oe.getSpecifier(), oe.getParts());
-        }
+        // convert line of sight
+        return convertLos(point, los);
     }
 
     /** Transform a cartesian point to a surface-relative point.
-     * @param point cartesian point
+     * @param point cartesian point (m)
      * @param frame frame in which cartesian point is expressed
      * @param date date of the computation (used for frames conversions)
      * @param centralLongitude reference longitude lc such that the point longitude will
-     * be normalized between lc-π and lc+π
+     * be normalized between lc-π and lc+π (rad)
      * @return point at the same location but as a surface-relative point
-     * @exception OrekitException if point cannot be converted to body frame
      */
     public NormalizedGeodeticPoint transform(final Vector3D point, final Frame frame, final AbsoluteDate date,
-                                             final double centralLongitude)
-        throws OrekitException {
+                                             final double centralLongitude) {
         final GeodeticPoint gp = transform(point, frame, date);
         return new NormalizedGeodeticPoint(gp.getLatitude(), gp.getLongitude(), gp.getAltitude(),
                                            centralLongitude);

src/main/java/org/orekit/rugged/utils/GridCreation.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.utils;
+
+/** Utility class for grids creation.
+ * @author Guylaine Prat
+ * @since 2.1
+ */
+public final class GridCreation {
+
+    /** Private constructor for utility class.
+     * Suppress default constructor for non instantiability ...
+     */
+    private GridCreation() {
+        super();
+    }
+
+    /** Create a linear grid between min and max value for a number n of points.
+     * TBN: no checks are performed here. Must be done by the calling method.
+     * @param min value for grid[0]
+     * @param max value for grid[n-1]
+     * @param n number of points
+     * @return the linear grid
+     */
+    public static double[] createLinearGrid(final double min, final double max, final int n) {
+
+        final double[] grid = new double[n];
+        for (int i = 0; i < n; ++i) {
+            grid[i] = ((n - 1 - i) * min + i * max) / (n - 1);
+        }
+        return grid;
+    }
+}

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

-/* Copyright 2013-2016 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

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

-/* Copyright 2013-2016 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

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

-/* Copyright 2013-2016 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

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

-/* Copyright 2013-2016 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
@@ -16,14 +16,13 @@
  */
 package org.orekit.rugged.utils;
 
-import org.hipparchus.geometry.euclidean.threed.Vector3D;
-import org.hipparchus.util.FastMath;
 import java.util.ArrayList;
 import java.util.List;
 
+import org.hipparchus.geometry.euclidean.threed.Vector3D;
+import org.hipparchus.util.FastMath;
 import org.orekit.bodies.GeodeticPoint;
 import org.orekit.bodies.OneAxisEllipsoid;
-import org.orekit.errors.OrekitException;
 import org.orekit.frames.Frame;
 import org.orekit.frames.Transform;
 import org.orekit.time.AbsoluteDate;
@@ -64,18 +63,16 @@ public class RoughVisibilityEstimator {
      * @param ellipsoid ground ellipsoid
      * @param frame frame in which position and velocity are defined (may be inertial or body frame)
      * @param positionsVelocities satellite position and velocity (m and m/s in specified frame)
-     * @exception OrekitException if position-velocity cannot be converted to body frame
      */
     public RoughVisibilityEstimator(final OneAxisEllipsoid ellipsoid, final Frame frame,
-                                    final List<TimeStampedPVCoordinates> positionsVelocities)
-        throws OrekitException {
+                                    final List<TimeStampedPVCoordinates> positionsVelocities) {
 
         this.ellipsoid = ellipsoid;
 
         // project spacecraft position-velocity to ground
         final Frame bodyFrame = ellipsoid.getBodyFrame();
         final int n = positionsVelocities.size();
-        this.pvGround = new ArrayList<TimeStampedPVCoordinates>(n);
+        this.pvGround = new ArrayList<>(n);
         for (final TimeStampedPVCoordinates pv : positionsVelocities) {
             final Transform t = frame.getTransformTo(bodyFrame, pv.getDate());
             pvGround.add(ellipsoid.projectToGround(t.transformPVCoordinates(pv), bodyFrame));

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

-/* Copyright 2013-2016 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