Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • orekit/rugged
  • sdinot/rugged
  • yzokras/rugged
  • youngcle/rugged-mod
4 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 1248 additions and 204 deletions
src/main/java/org/orekit/rugged/raster/Tile.java
View file @ acf902ea
/* 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
......@@ -17,8 +17,6 @@
package org.orekit.rugged.raster;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.rugged.errors.RuggedException;
import org.orekit.rugged.utils.NormalizedGeodeticPoint;
/** Interface representing a raster tile.
......@@ -29,6 +27,7 @@ import org.orekit.rugged.utils.NormalizedGeodeticPoint;
* correspond to the <em>center</em> of the most North-East cell.
* </p>
* @author Luc Maisonobe
* @author Guylaine Prat
*/
public interface Tile extends UpdatableTile {
......@@ -113,20 +112,18 @@ public interface Tile extends UpdatableTile {
}
/** 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 of grid interpolation points.
* @return 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);
......@@ -140,20 +137,20 @@ public interface Tile extends UpdatableTile {
* {@link Location#NORTH} or {@link Location#NORTH_EAST}, but can
* <em>never</em> return {@link Location#HAS_INTERPOLATION_NEIGHBORS}!
* </p>
* @return maximum latitude
* @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);
......@@ -167,18 +164,18 @@ public interface Tile extends UpdatableTile {
* {@link Location#EAST} or {@link Location#NORTH_EAST}, but can
* <em>never</em> return {@link Location#HAS_INTERPOLATION_NEIGHBORS}!
* </p>
* @return maximum longitude
* @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();
......@@ -211,7 +208,7 @@ public interface Tile extends UpdatableTile {
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();
......@@ -224,7 +221,7 @@ public interface Tile extends UpdatableTile {
int getMinElevationLongitudeIndex();
/** Get the maximum elevation in the tile.
* @return maximum elevation in the tile
* @return maximum elevation in the tile (m)
*/
double getMaxElevation();
......@@ -239,11 +236,9 @@ public interface Tile extends UpdatableTile {
/** 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>
......@@ -256,25 +251,33 @@ public interface Tile extends UpdatableTile {
* </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 cell.
* @param p point on the line
* <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 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 cell, null otherwise
* @exception RuggedException if intersection point cannot be computed
*/
NormalizedGeodeticPoint cellIntersection(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/raster/TileFactory.java
View file @ acf902ea
/* 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
View file @ acf902ea
/* 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
View file @ acf902ea
/* 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
View file @ acf902ea
/* 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 0 → 100644
View file @ acf902ea
/* 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 0 → 100644
View file @ acf902ea
/* 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 0 → 100644
View file @ acf902ea
/* 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 0 → 100644
View file @ acf902ea
/* 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 0 → 100644
View file @ acf902ea
/* 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 0 → 100644
View file @ acf902ea
/* 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 0 → 100644
View file @ acf902ea
/* 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
View file @ acf902ea
/* 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 0 → 100644
View file @ acf902ea
/* 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
View file @ acf902ea
/* 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 0 → 100644
View file @ acf902ea
/* 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
View file @ acf902ea
/* 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
View file @ acf902ea
/* 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
View file @ acf902ea
/* 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
View file @ acf902ea
/* 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));
......