diff --git a/Jenkinsfile b/Jenkinsfile new file mode 100644 index 0000000000000000000000000000000000000000..9f68828414b5c24d339f76984d255219eb84a274 --- /dev/null +++ b/Jenkinsfile @@ -0,0 +1,48 @@ +pipeline { + + agent any + tools { + maven 'mvn-default' + jdk 'openjdk-8' + } + + options { + timeout(time: 60, unit: 'MINUTES') + } + + stages { + + stage('Cleaning') { + steps { + sh 'git clean -fdx' + } + } + + stage('Build') { + steps { + script { + if ( env.BRANCH_NAME ==~ /^release-[.0-9]+$/ ) { + sh 'mvn verify assembly:single' + } + else { + sh 'mvn verify site' + } + } + } + } + } + + post { + always { + archiveArtifacts artifacts: 'target/*.jar', fingerprint: true + script { + if ( env.BRANCH_NAME ==~ /^release-[.0-9]+$/ ) { + archiveArtifacts artifacts: 'target/*.zip', fingerprint: true + } + } + checkstyle pattern: 'target/checkstyle-result.xml' + junit 'target/surefire-reports/*.xml' + jacoco execPattern:'target/**.exec', classPattern: '**/classes', sourcePattern: '**/src/main/java' + } + } +} diff --git a/checkstyle.xml b/checkstyle.xml index 6314d6616c6d2ee4a33ba0e12d3ce56255a5eee1..7fc65ee528ab744db95b38dbad022ac812bc301b 100644 --- a/checkstyle.xml +++ b/checkstyle.xml @@ -87,9 +87,19 @@ <property name="checkFormat" value="NoWhitespaceAfter"/> </module> <module name="SuppressionCommentFilter"> - <property name="offCommentFormat" value="CHECKSTYLE\: stop FallThrough check"/> - <property name="onCommentFormat" value="CHECKSTYLE\: resume FallThrough check"/> - <property name="checkFormat" value="FallThrough"/> + <property name="offCommentFormat" value="CHECKSTYLE\: stop Indentation check"/> + <property name="onCommentFormat" value="CHECKSTYLE\: resume Indentation check"/> + <property name="checkFormat" value="IndentationCheck"/> + </module> + <module name="SuppressionCommentFilter"> + <property name="offCommentFormat" value="CHECKSTYLE\: stop MultipleStringLiterals check"/> + <property name="onCommentFormat" value="CHECKSTYLE\: resume MultipleStringLiterals check"/> + <property name="checkFormat" value="MultipleStringLiteralsCheck"/> + </module> + <module name="SuppressionCommentFilter"> + <property name="offCommentFormat" value="CHECKSTYLE\: stop UnnecessaryParentheses check"/> + <property name="onCommentFormat" value="CHECKSTYLE\: resume UnnecessaryParentheses check"/> + <property name="checkFormat" value="UnnecessaryParentheses"/> </module> </module> <module name="RegexpHeader"> diff --git a/pom.xml b/pom.xml index b024b65209560cb5ccf19793eff8f052aaf41ac8..31ea7a972e83a0613a943f3f6586f79757af539e 100644 --- a/pom.xml +++ b/pom.xml @@ -18,32 +18,33 @@ <properties> <project.build.sourceEncoding>UTF-8</project.build.sourceEncoding> <project.reporting.outputEncoding>UTF-8</project.reporting.outputEncoding> - <rugged.findbugs-maven-plugin.version>3.0.4</rugged.findbugs-maven-plugin.version> - <rugged.jacoco-maven-plugin.version>0.7.9</rugged.jacoco-maven-plugin.version> - <rugged.maven-assembly-plugin.version>3.0.0</rugged.maven-assembly-plugin.version> - <rugged.maven-bundle-plugin.version>3.3.0</rugged.maven-bundle-plugin.version> + <rugged.spotbugs-maven-plugin.version>3.1.7</rugged.spotbugs-maven-plugin.version> + <rugged.jacoco-maven-plugin.version>0.8.2</rugged.jacoco-maven-plugin.version> + <rugged.maven-assembly-plugin.version>3.1.0</rugged.maven-assembly-plugin.version> + <rugged.maven-bundle-plugin.version>4.1.0</rugged.maven-bundle-plugin.version> <rugged.maven-changes-plugin.version>2.12.1</rugged.maven-changes-plugin.version> - <rugged.maven-checkstyle-plugin.version>2.17</rugged.maven-checkstyle-plugin.version> - <rugged.checkstyle.version>8.10</rugged.checkstyle.version> - <rugged.maven-clean-plugin.version>3.0.0</rugged.maven-clean-plugin.version> - <rugged.maven-compiler-plugin.version>3.6.1</rugged.maven-compiler-plugin.version> - <rugged.maven-javadoc-plugin.version>2.10.4</rugged.maven-javadoc-plugin.version> - <rugged.maven-jar-plugin.version>3.0.2</rugged.maven-jar-plugin.version> - <rugged.maven-jxr-plugin.version>2.5</rugged.maven-jxr-plugin.version> + <rugged.maven-checkstyle-plugin.version>3.0.0</rugged.maven-checkstyle-plugin.version> + <rugged.checkstyle.version>8.14</rugged.checkstyle.version> + <rugged.maven-clean-plugin.version>3.1.0</rugged.maven-clean-plugin.version> + <rugged.maven-compiler-plugin.version>3.8.0</rugged.maven-compiler-plugin.version> + <rugged.maven-javadoc-plugin.version>3.0.1</rugged.maven-javadoc-plugin.version> + <rugged.maven-jar-plugin.version>3.1.0</rugged.maven-jar-plugin.version> + <rugged.maven-jxr-plugin.version>3.0.0</rugged.maven-jxr-plugin.version> <rugged.plantuml-maven-plugin.version>1.2</rugged.plantuml-maven-plugin.version> - <rugged.plantuml.version>1.2017.15</rugged.plantuml.version> - <rugged.maven-project-info-reports-plugin.version>2.9</rugged.maven-project-info-reports-plugin.version> - <rugged.maven-resources-plugin.version>3.0.2</rugged.maven-resources-plugin.version> - <rugged.maven-site-plugin.version>3.6</rugged.maven-site-plugin.version> - <rugged.maven-surefire-plugin.version>2.20</rugged.maven-surefire-plugin.version> - <rugged.maven-surefire-report-plugin.version>2.20</rugged.maven-surefire-report-plugin.version> + <rugged.plantuml.version>1.2018.12</rugged.plantuml.version> + <rugged.maven-project-info-reports-plugin.version>3.0.0</rugged.maven-project-info-reports-plugin.version> + <rugged.maven-resources-plugin.version>3.1.0</rugged.maven-resources-plugin.version> + <rugged.maven-site-plugin.version>3.7.1</rugged.maven-site-plugin.version> + <rugged.maven-source-plugin.version>3.0.1</rugged.maven-source-plugin.version> + <rugged.maven-surefire-plugin.version>2.22.1</rugged.maven-surefire-plugin.version> + <rugged.maven-surefire-report-plugin.version>2.22.1</rugged.maven-surefire-report-plugin.version> <rugged.jgit.buildnumber.version>1.2.10</rugged.jgit.buildnumber.version> <rugged.build-helper-maven-plugin.version>3.0.0</rugged.build-helper-maven-plugin.version> <rugged.nexus-staging-maven-plugin.version>1.6.8</rugged.nexus-staging-maven-plugin.version> <rugged.maven-gpg-plugin.version>1.6</rugged.maven-gpg-plugin.version> - <rugged.maven-install-plugin.version>2.5.2</rugged.maven-install-plugin.version> - <rugged.orekit.version>9.2</rugged.orekit.version> - <rugged.hipparchus.version>1.3</rugged.hipparchus.version> + <rugged.maven-install-plugin.version>3.0.0-M1</rugged.maven-install-plugin.version> + <rugged.orekit.version>9.3-SNAPSHOT</rugged.orekit.version> + <rugged.hipparchus.version>1.4</rugged.hipparchus.version> <rugged.junit.version>4.12</rugged.junit.version> <rugged.compiler.source>1.8</rugged.compiler.source> <rugged.compiler.target>1.8</rugged.compiler.target> @@ -129,10 +130,10 @@ </scm> <issueManagement> - <system>Redmine</system> - <url>https://www.orekit.org/forge/projects/rugged/issues</url> + <system>Gitlab</system> + <url>https://gitlab.orekit.org/orekit/rugged/issues</url> </issueManagement> - + <mailingLists> <mailingList> <name>Shared Orekit/Rugged announces mailing list</name> @@ -448,14 +449,14 @@ <version>${rugged.maven-project-info-reports-plugin.version}</version> </plugin> <plugin> - <groupId>org.codehaus.mojo</groupId> - <artifactId>findbugs-maven-plugin</artifactId> - <version>${rugged.findbugs-maven-plugin.version}</version> + <groupId>com.github.spotbugs</groupId> + <artifactId>spotbugs-maven-plugin</artifactId> + <version>${rugged.spotbugs-maven-plugin.version}</version> <configuration> <threshold>Normal</threshold> <effort>Default</effort> - <excludeFilterFile>${basedir}/findbugs-exclude-filter.xml</excludeFilterFile> - </configuration> + <excludeFilterFile>${basedir}/spotbugs-exclude-filter.xml</excludeFilterFile> + </configuration> </plugin> <plugin> <groupId>org.apache.maven.plugins</groupId> @@ -584,6 +585,21 @@ <id>release</id> <build> <plugins> + <!-- + <plugin> + <groupId>org.apache.maven.plugins</groupId> + <artifactId>maven-source-plugin</artifactId> + <version>${orekit.maven-source-plugin.version}</version> + <executions> + <execution> + <id>attach-sources</id> + <goals> + <goal>jar-no-fork</goal> + </goals> + </execution> + </executions> + </plugin> + --> <plugin> <groupId>org.apache.maven.plugins</groupId> <artifactId>maven-javadoc-plugin</artifactId> diff --git a/findbugs-exclude-filter.xml b/spotbugs-exclude-filter.xml similarity index 76% rename from findbugs-exclude-filter.xml rename to spotbugs-exclude-filter.xml index ea36a9e9c86089cdee6bb556df0ee49d9a2d691b..f06c28a6db3204c09a1c41ab9da95459eee03108 100644 --- a/findbugs-exclude-filter.xml +++ b/spotbugs-exclude-filter.xml @@ -1,6 +1,6 @@ <?xml version="1.0"?> <!-- - This file contains some false positive bugs detected by findbugs. Their + This file contains some false positive bugs detected by spotbugs. Their false positive nature has been analyzed individually and they have been put here to instruct findbugs it must ignore them. --> diff --git a/src/main/assembly/source-distribution-assembly.xml b/src/main/assembly/source-distribution-assembly.xml index 91abe41918453dd3a7acea9970c8745a7ede56c4..c7f81d8b88bc8dfa23c5c9b0e9ab14ef1d51de83 100644 --- a/src/main/assembly/source-distribution-assembly.xml +++ b/src/main/assembly/source-distribution-assembly.xml @@ -12,7 +12,7 @@ <include>BUILDING.txt</include> <include>pom.xml</include> <include>checkstyle.xml</include> - <include>findbugs-exclude-filter.xml</include> + <include>spotbugs-exclude-filter.xml</include> <include>license-header.txt</include> </includes> <useDefaultExcludes>true</useDefaultExcludes> @@ -22,4 +22,4 @@ <useDefaultExcludes>true</useDefaultExcludes> </fileSet> </fileSets> -</assembly> \ No newline at end of file +</assembly> diff --git a/src/main/java/org/orekit/rugged/api/Rugged.java b/src/main/java/org/orekit/rugged/api/Rugged.java index 95f733edc2de2235f628ec264ddbe35a0d529f01..c3b26db71a785847cd1ca49a64008ce548b27453 100644 --- a/src/main/java/org/orekit/rugged/api/Rugged.java +++ b/src/main/java/org/orekit/rugged/api/Rugged.java @@ -1,4 +1,4 @@ -/* Copyright 2013-2017 CS Systèmes d'Information +/* Copyright 2013-2018 CS Systèmes d'Information * Licensed to CS Systèmes d'Information (CS) under one or more * contributor license agreements. See the NOTICE file distributed with * this work for additional information regarding copyright ownership. @@ -34,6 +34,7 @@ import org.orekit.rugged.linesensor.LineSensor; import org.orekit.rugged.linesensor.SensorMeanPlaneCrossing; import org.orekit.rugged.linesensor.SensorPixel; import org.orekit.rugged.linesensor.SensorPixelCrossing; +import org.orekit.rugged.los.PixelLOS; import org.orekit.rugged.refraction.AtmosphericRefraction; import org.orekit.rugged.utils.DSGenerator; import org.orekit.rugged.utils.ExtendedEllipsoid; @@ -46,9 +47,9 @@ import org.orekit.utils.PVCoordinates; /** Main class of Rugged library API. * @see RuggedBuilder * @author Luc Maisonobe - * @author Lucie LabatAllee - * @author Jonathan Guinet * @author Guylaine Prat + * @author Jonathan Guinet + * @author Lucie LabatAllee */ public class Rugged { @@ -60,9 +61,20 @@ public class Rugged { */ private static final double COARSE_INVERSE_LOCATION_ACCURACY = 0.01; - /** Maximum number of evaluations. */ + /** Maximum number of evaluations for crossing algorithms. */ private static final int MAX_EVAL = 50; + /** Margin for computation of inverse location with atmospheric refraction correction. */ + private static final double INVLOC_MARGIN = 0.5; + + /** Threshold for pixel convergence in fixed point method + * (for inverse location with atmospheric refraction correction). */ + private static final double PIXEL_CV_THRESHOLD = 1.e-4; + + /** Threshold for line convergence in fixed point method + * (for inverse location with atmospheric refraction correction). */ + private static final double LINE_CV_THRESHOLD = 1.e-4; + /** Reference ellipsoid. */ private final ExtendedEllipsoid ellipsoid; @@ -231,178 +243,188 @@ public class Rugged { public GeodeticPoint[] directLocation(final String sensorName, final double lineNumber) throws RuggedException { - // compute the approximate transform between spacecraft and observed body - final LineSensor sensor = getLineSensor(sensorName); - final AbsoluteDate date = sensor.getDate(lineNumber); + final LineSensor sensor = getLineSensor(sensorName); + final Vector3D sensorPosition = sensor.getPosition(); + final AbsoluteDate date = sensor.getDate(lineNumber); + + // Compute the transform for the date + // from spacecraft to inertial final Transform scToInert = scToBody.getScToInertial(date); + // from inertial to body final Transform inertToBody = scToBody.getInertialToBody(date); - final Transform approximate = new Transform(date, scToInert, inertToBody); - final Vector3D spacecraftVelocity = - scToInert.transformPVCoordinates(PVCoordinates.ZERO).getVelocity(); + // Compute spacecraft velocity in inertial frame + final Vector3D spacecraftVelocity = scToInert.transformPVCoordinates(PVCoordinates.ZERO).getVelocity(); + // Compute sensor position in inertial frame + // TBN: for simplicity, due to the size of sensor, we consider each pixel to be at sensor position + final Vector3D pInert = scToInert.transformPosition(sensorPosition); - // compute location of each pixel - final Vector3D pInert = scToInert.transformPosition(sensor.getPosition()); + // Compute location of each pixel final GeodeticPoint[] gp = new GeodeticPoint[sensor.getNbPixels()]; for (int i = 0; i < sensor.getNbPixels(); ++i) { - DumpManager.dumpDirectLocation(date, sensor.getPosition(), sensor.getLOS(date, i), lightTimeCorrection, - aberrationOfLightCorrection, atmosphericRefraction != null); + final Vector3D los = sensor.getLOS(date, i); + DumpManager.dumpDirectLocation(date, sensorPosition, los, lightTimeCorrection, + aberrationOfLightCorrection, atmosphericRefraction != null); - final Vector3D obsLInert = scToInert.transformVector(sensor.getLOS(date, i)); + // compute the line of sight in inertial frame (without correction) + final Vector3D obsLInert = scToInert.transformVector(los); final Vector3D lInert; + if (aberrationOfLightCorrection) { - // apply aberration of light correction - // as the spacecraft velocity is small with respect to speed of light, - // we use classical velocity addition and not relativistic velocity addition - // we look for a positive k such that: c * lInert + vsat = k * obsLInert - // with lInert normalized - final double a = obsLInert.getNormSq(); - final double b = -Vector3D.dotProduct(obsLInert, spacecraftVelocity); - final double c = spacecraftVelocity.getNormSq() - Constants.SPEED_OF_LIGHT * Constants.SPEED_OF_LIGHT; - final double s = FastMath.sqrt(b * b - a * c); - final double k = (b > 0) ? -c / (s + b) : (s - b) / a; - lInert = new Vector3D( k / Constants.SPEED_OF_LIGHT, obsLInert, - -1.0 / Constants.SPEED_OF_LIGHT, spacecraftVelocity); + // apply aberration of light correction on LOS + lInert = applyAberrationOfLightCorrection(obsLInert, spacecraftVelocity); } else { - // don't apply aberration of light correction + // don't apply aberration of light correction on LOS lInert = obsLInert; } if (lightTimeCorrection) { // compute DEM intersection with light time correction - final Vector3D sP = approximate.transformPosition(sensor.getPosition()); - final Vector3D sL = approximate.transformVector(sensor.getLOS(date, i)); - final Vector3D eP1 = ellipsoid.transform(ellipsoid.pointOnGround(sP, sL, 0.0)); - final double deltaT1 = eP1.distance(sP) / Constants.SPEED_OF_LIGHT; - final Transform shifted1 = inertToBody.shiftedBy(-deltaT1); - final NormalizedGeodeticPoint gp1 = algorithm.intersection(ellipsoid, - shifted1.transformPosition(pInert), - shifted1.transformVector(lInert)); - - final Vector3D eP2 = ellipsoid.transform(gp1); - final double deltaT2 = eP2.distance(sP) / Constants.SPEED_OF_LIGHT; - final Transform shifted2 = inertToBody.shiftedBy(-deltaT2); - gp[i] = algorithm.refineIntersection(ellipsoid, - shifted2.transformPosition(pInert), - shifted2.transformVector(lInert), - gp1); + // TBN: for simplicity, due to the size of sensor, we consider each pixel to be at sensor position + gp[i] = computeWithLightTimeCorrection(date, sensorPosition, los, scToInert, inertToBody, pInert, lInert); } else { // compute DEM intersection without light time correction final Vector3D pBody = inertToBody.transformPosition(pInert); final Vector3D lBody = inertToBody.transformVector(lInert); gp[i] = algorithm.refineIntersection(ellipsoid, pBody, lBody, - algorithm.intersection(ellipsoid, pBody, lBody)); + algorithm.intersection(ellipsoid, pBody, lBody)); } - if (atmosphericRefraction != null) { - // apply atmospheric refraction correction - final Vector3D pBody = inertToBody.transformPosition(pInert); - final Vector3D lBody = inertToBody.transformVector(lInert); - gp[i] = atmosphericRefraction.applyCorrection(pBody, lBody, (NormalizedGeodeticPoint) gp[i], algorithm); - } + // compute with atmospheric refraction correction if necessary + if (atmosphericRefraction != null && atmosphericRefraction.mustBeComputed()) { - DumpManager.dumpDirectLocationResult(gp[i]); + // Test if optimization is not required + if (!atmosphericRefraction.isOptimized()) { - } + // apply atmospheric refraction correction + final Vector3D pBody = inertToBody.transformPosition(pInert); + final Vector3D lBody = inertToBody.transformVector(lInert); + gp[i] = atmosphericRefraction.applyCorrection(pBody, lBody, (NormalizedGeodeticPoint) gp[i], algorithm); - return gp; + } else { // Optimization is required + // TODO algo with optimization + throw new RuggedException(RuggedMessages.UNINITIALIZED_CONTEXT, "atmospheric optimization not defined"); + } + } + DumpManager.dumpDirectLocationResult(gp[i]); + } + return gp; } /** Direct location of a single line-of-sight. + * <br> + * NB: if (the optionnal) atmospheric refraction must be computed with the optimization algorithm, + * see {@link #directLocation(AbsoluteDate, Vector3D, PixelLOS)}. * @param date date of the location - * @param position pixel position in spacecraft frame + * @param sensorPosition sensor position in spacecraft frame. For simplicity, due to the size of sensor, + * we consider each pixel to be at sensor position * @param los normalized line-of-sight in spacecraft frame * @return ground position of intersection point between specified los and ground - * @exception RuggedException if line cannot be localized, or sensor is unknown + * @exception RuggedException if line cannot be localized, sensor is unknown or problem with atmospheric data */ - public GeodeticPoint directLocation(final AbsoluteDate date, final Vector3D position, final Vector3D los) + public GeodeticPoint directLocation(final AbsoluteDate date, final Vector3D sensorPosition, final Vector3D los) throws RuggedException { - DumpManager.dumpDirectLocation(date, position, los, lightTimeCorrection, aberrationOfLightCorrection, + // Set the pixel to null in order not to compute atmosphere with optimization + final PixelLOS pixelLOS = new PixelLOS(null, los); + + return directLocation(date, sensorPosition, pixelLOS); + } + + + /** Direct location of a single line-of-sight. + * <br> + * NB: if the atmospheric refraction must be computed without the optimization algorithm, + * see {@link #directLocation(AbsoluteDate, Vector3D, Vector3D)}. + * @param date date of the location + * @param sensorPosition sensor position in spacecraft frame. For simplicity, due to the size of sensor, + * we consider each pixel to be at sensor position + * @param pixelLOS pixel definition with normalized line-of-sight in spacecraft frame + * @return ground position of intersection point between specified los and ground + * @exception RuggedException if line cannot be localized, sensor is unknown or problem with atmospheric data + */ + public GeodeticPoint directLocation(final AbsoluteDate date, final Vector3D sensorPosition, final PixelLOS pixelLOS) + throws RuggedException { + + final Vector3D los = pixelLOS.getLOS(); + // TODO change dump to add sensorPixel + DumpManager.dumpDirectLocation(date, sensorPosition, los, lightTimeCorrection, aberrationOfLightCorrection, atmosphericRefraction != null); - // compute the approximate transform between spacecraft and observed body + // Compute the transforms for the date + // from spacecraft to inertial final Transform scToInert = scToBody.getScToInertial(date); + // from inertial to body final Transform inertToBody = scToBody.getInertialToBody(date); - final Transform approximate = new Transform(date, scToInert, inertToBody); - final Vector3D spacecraftVelocity = - scToInert.transformPVCoordinates(PVCoordinates.ZERO).getVelocity(); - - // compute location of specified pixel - final Vector3D pInert = scToInert.transformPosition(position); + // Compute spacecraft velocity in inertial frame + final Vector3D spacecraftVelocity = scToInert.transformPVCoordinates(PVCoordinates.ZERO).getVelocity(); + // Compute sensor position in inertial frame + // TBN: for simplicity, due to the size of sensor, we consider each pixel to be at sensor position + final Vector3D pInert = scToInert.transformPosition(sensorPosition); + // Compute the line of sight in inertial frame (without correction) final Vector3D obsLInert = scToInert.transformVector(los); + final Vector3D lInert; if (aberrationOfLightCorrection) { - // apply aberration of light correction - // as the spacecraft velocity is small with respect to speed of light, - // we use classical velocity addition and not relativistic velocity addition - // we look for a positive k such that: c * lInert + vsat = k * obsLInert - // with lInert normalized - final double a = obsLInert.getNormSq(); - final double b = -Vector3D.dotProduct(obsLInert, spacecraftVelocity); - final double c = spacecraftVelocity.getNormSq() - Constants.SPEED_OF_LIGHT * Constants.SPEED_OF_LIGHT; - final double s = FastMath.sqrt(b * b - a * c); - final double k = (b > 0) ? -c / (s + b) : (s - b) / a; - lInert = new Vector3D( k / Constants.SPEED_OF_LIGHT, obsLInert, - -1.0 / Constants.SPEED_OF_LIGHT, spacecraftVelocity); + // apply aberration of light correction on LOS + lInert = applyAberrationOfLightCorrection(obsLInert, spacecraftVelocity); } else { - // don't apply aberration of light correction + // don't apply aberration of light correction on LOS lInert = obsLInert; } + // Compute ground location of specified pixel final NormalizedGeodeticPoint gp; + if (lightTimeCorrection) { // compute DEM intersection with light time correction - final Vector3D sP = approximate.transformPosition(position); - final Vector3D sL = approximate.transformVector(los); - final Vector3D eP1 = ellipsoid.transform(ellipsoid.pointOnGround(sP, sL, 0.0)); - final double deltaT1 = eP1.distance(sP) / Constants.SPEED_OF_LIGHT; - final Transform shifted1 = inertToBody.shiftedBy(-deltaT1); - final NormalizedGeodeticPoint gp1 = algorithm.intersection(ellipsoid, - shifted1.transformPosition(pInert), - shifted1.transformVector(lInert)); - - final Vector3D eP2 = ellipsoid.transform(gp1); - final double deltaT2 = eP2.distance(sP) / Constants.SPEED_OF_LIGHT; - final Transform shifted2 = inertToBody.shiftedBy(-deltaT2); - gp = algorithm.refineIntersection(ellipsoid, - shifted2.transformPosition(pInert), - shifted2.transformVector(lInert), - gp1); + // TBN: for simplicity, due to the size of sensor, we consider each pixel to be at sensor position + gp = computeWithLightTimeCorrection(date, sensorPosition, los, scToInert, inertToBody, pInert, lInert); } else { // compute DEM intersection without light time correction final Vector3D pBody = inertToBody.transformPosition(pInert); final Vector3D lBody = inertToBody.transformVector(lInert); gp = algorithm.refineIntersection(ellipsoid, pBody, lBody, - algorithm.intersection(ellipsoid, pBody, lBody)); + algorithm.intersection(ellipsoid, pBody, lBody)); } - final NormalizedGeodeticPoint result; - if (atmosphericRefraction != null) { - // apply atmospheric refraction correction - final Vector3D pBody = inertToBody.transformPosition(pInert); - final Vector3D lBody = inertToBody.transformVector(lInert); - result = atmosphericRefraction.applyCorrection(pBody, lBody, gp, algorithm); - } else { - // don't apply atmospheric refraction correction - result = gp; - } + NormalizedGeodeticPoint result = gp; + + // compute the ground location with atmospheric correction if asked for + if (atmosphericRefraction != null && atmosphericRefraction.mustBeComputed()) { + // TBN: two methods exist for computation of the atmospheric correction + // * a full computation + // * a time optimized computation based on an interpolation grid + + // Test if optimization is not required or if sensor pixel is not defined (impossible to perform optimization) + if (!atmosphericRefraction.isOptimized() || pixelLOS.getSensorPixel() == null) { + + // apply atmospheric refraction correction (full computation) + final Vector3D pBody = inertToBody.transformPosition(pInert); + final Vector3D lBody = inertToBody.transformVector(lInert); + result = atmosphericRefraction.applyCorrection(pBody, lBody, gp, algorithm); + + } else { // Optimization is required and sensor pixel is defined + + // TODO to be done + throw new RuggedException(RuggedMessages.UNINITIALIZED_CONTEXT, "Atmospheric optimization for direct loc"); + + } // end test on optimization is required + } // end test on atmosphericRefraction != null DumpManager.dumpDirectLocationResult(result); return result; - } /** Find the date at which sensor sees a ground point. * <p> - * This method is a partial {@link #inverseLocation(String, - * GeodeticPoint, int, int) inverse location} focusing only on date. + * This method is a partial {@link #inverseLocation(String, GeodeticPoint, int, int) inverse location} focusing only on date. * </p> * <p> * The point is given only by its latitude and longitude, the elevation is @@ -483,7 +505,6 @@ public class Rugged { } else { return sensor.getDate(crossingResult.getLine()); } - } /** Inverse location of a ground point. @@ -518,8 +539,7 @@ public class Rugged { final double latitude, final double longitude, final int minLine, final int maxLine) throws RuggedException { - final GeodeticPoint groundPoint = - new GeodeticPoint(latitude, longitude, algorithm.getElevation(latitude, longitude)); + final GeodeticPoint groundPoint = new GeodeticPoint(latitude, longitude, algorithm.getElevation(latitude, longitude)); return inverseLocation(sensorName, groundPoint, minLine, maxLine); } @@ -538,9 +558,9 @@ public class Rugged { * are only an example and should be adjusted depending on mission needs. * </p> * @param sensorName name of the line sensor - * @param point point to localize - * @param minLine minimum line number - * @param maxLine maximum line number + * @param point geodetic point to localize + * @param minLine minimum line number where the search will be performed + * @param maxLine maximum line number where the search will be performed * @return sensor pixel seeing point, or null if point cannot be seen between the * prescribed line numbers * @exception RuggedException if line cannot be localized, or sensor is unknown @@ -552,15 +572,99 @@ public class Rugged { throws RuggedException { final LineSensor sensor = getLineSensor(sensorName); - DumpManager.dumpInverseLocation(sensor, point, minLine, maxLine, - lightTimeCorrection, aberrationOfLightCorrection); + DumpManager.dumpInverseLocation(sensor, point, minLine, maxLine, lightTimeCorrection, aberrationOfLightCorrection); final SensorMeanPlaneCrossing planeCrossing = getPlaneCrossing(sensorName, minLine, maxLine); - DumpManager.dumpSensorMeanPlane(planeCrossing); + if (atmosphericRefraction == null || !atmosphericRefraction.mustBeComputed()) { + // Compute inverse location WITHOUT atmospheric refraction + return findSensorPixelWithoutAtmosphere(point, sensor, planeCrossing); + } else { + // Compute inverse location WITH atmospheric refraction + return findSensorPixelWithAtmosphere(point, sensor, minLine, maxLine); + } + } + + /** Apply aberration of light correction (for direct location). + * @param spacecraftVelocity spacecraft velocity in inertial frame + * @param obsLInert line of sight in inertial frame + * @return line of sight with aberration of light correction + */ + private Vector3D applyAberrationOfLightCorrection(final Vector3D obsLInert, final Vector3D spacecraftVelocity) { + + // As the spacecraft velocity is small with respect to speed of light, + // we use classical velocity addition and not relativistic velocity addition + // we look for a positive k such that: c * lInert + vsat = k * obsLInert + // with lInert normalized + final double a = obsLInert.getNormSq(); + final double b = -Vector3D.dotProduct(obsLInert, spacecraftVelocity); + final double c = spacecraftVelocity.getNormSq() - Constants.SPEED_OF_LIGHT * Constants.SPEED_OF_LIGHT; + + // a > 0 and c < 0 + final double s = FastMath.sqrt(b * b - a * c); + + // Only the k > 0 are kept as solutions (the solutions: -(s+b)/a and c/(s-b) are useless) + final double k = (b > 0) ? -c / (s + b) : (s - b) / a; + + final Vector3D lInert = new Vector3D( k / Constants.SPEED_OF_LIGHT, obsLInert, -1.0 / Constants.SPEED_OF_LIGHT, spacecraftVelocity); + return lInert; + } + + /** Compute the DEM intersection with light time correction. + * @param date date of the los + * @param sensorPosition sensor position in spacecraft frame + * @param los los in spacecraft frame + * @param scToInert transform for the date from spacecraft to inertial + * @param inertToBody transform for the date from inertial to body + * @param pInert sensor position in inertial frame + * @param lInert line of sight in inertial frame + * @return geodetic point with light time correction + * @throws RuggedException if intersection cannot be found + */ + private NormalizedGeodeticPoint computeWithLightTimeCorrection(final AbsoluteDate date, + final Vector3D sensorPosition, final Vector3D los, + final Transform scToInert, final Transform inertToBody, + final Vector3D pInert, final Vector3D lInert) + throws RuggedException { + + // compute the approximate transform between spacecraft and observed body + final Transform approximate = new Transform(date, scToInert, inertToBody); + + final Vector3D sL = approximate.transformVector(los); + final Vector3D sP = approximate.transformPosition(sensorPosition); + + final Vector3D eP1 = ellipsoid.transform(ellipsoid.pointOnGround(sP, sL, 0.0)); + final double deltaT1 = eP1.distance(sP) / Constants.SPEED_OF_LIGHT; + final Transform shifted1 = inertToBody.shiftedBy(-deltaT1); + final NormalizedGeodeticPoint gp1 = algorithm.intersection(ellipsoid, + shifted1.transformPosition(pInert), + shifted1.transformVector(lInert)); + + final Vector3D eP2 = ellipsoid.transform(gp1); + final double deltaT2 = eP2.distance(sP) / Constants.SPEED_OF_LIGHT; + final Transform shifted2 = inertToBody.shiftedBy(-deltaT2); + return algorithm.refineIntersection(ellipsoid, + shifted2.transformPosition(pInert), + shifted2.transformVector(lInert), + gp1); + } + + /** + * Find the sensor pixel WITHOUT atmospheric refraction correction. + * @param point geodetic point to localize + * @param sensor the line sensor + * @param planeCrossing the sensor mean plane crossing + * @return the sensor pixel crossing or null if cannot be found + * @throws RuggedException if sensor cannot be found + * @since 3.0 + */ + private SensorPixel findSensorPixelWithoutAtmosphere(final GeodeticPoint point, + final LineSensor sensor, final SensorMeanPlaneCrossing planeCrossing) + throws RuggedException { + // find approximately the sensor line at which ground point crosses sensor mean plane - final Vector3D target = ellipsoid.transform(point); + final Vector3D target = ellipsoid.transform(point); final SensorMeanPlaneCrossing.CrossingResult crossingResult = planeCrossing.find(target); if (crossingResult == null) { // target is out of search interval @@ -584,10 +688,8 @@ public class Rugged { final Vector3D lowLOS = sensor.getLOS(crossingResult.getDate(), lowIndex); final Vector3D highLOS = sensor.getLOS(crossingResult.getDate(), lowIndex + 1); final Vector3D localZ = Vector3D.crossProduct(lowLOS, highLOS).normalize(); - final double beta = FastMath.acos(Vector3D.dotProduct(crossingResult.getTargetDirection(), - localZ)); - final double s = Vector3D.dotProduct(crossingResult.getTargetDirectionDerivative(), - localZ); + final double beta = FastMath.acos(Vector3D.dotProduct(crossingResult.getTargetDirection(), localZ)); + final double s = Vector3D.dotProduct(crossingResult.getTargetDirectionDerivative(), localZ); final double betaDer = -s / FastMath.sqrt(1 - s * s); final double deltaL = (0.5 * FastMath.PI - beta) / betaDer; final double fixedLine = crossingResult.getLine() + deltaL; @@ -609,8 +711,187 @@ public class Rugged { final SensorPixel result = new SensorPixel(fixedLine, fixedPixel); DumpManager.dumpInverseLocationResult(result); + return result; + } + /** + * Find the sensor pixel WITH atmospheric refraction correction. + * @param point geodetic point to localize + * @param sensor the line sensor + * @param minLine minimum line number where the search will be performed + * @param maxLine maximum line number where the search will be performed + * @return the sensor pixel crossing or null if cannot be found + * @throws RuggedException if problem while computing correction + * @since 3.0 + */ + private SensorPixel findSensorPixelWithAtmosphere(final GeodeticPoint point, + final LineSensor sensor, final int minLine, final int maxLine) + throws RuggedException { + + // TBN: there is no direct way to compute the inverse location. + // The method is based on an interpolation grid associated with the fixed point method + + final String sensorName = sensor.getName(); + + // Compute a correction grid (at sensor level) + // =========================================== + // Need to be computed only once for a given sensor (with the same minLine and maxLine) + if (atmosphericRefraction.getBifPixel() == null || atmosphericRefraction.getBifLine() == null || // lazy evaluation + (!atmosphericRefraction.isSameContext(sensorName, minLine, maxLine))) { // Must be recomputed if the context changed + + // Definition of a regular grid (at sensor level) + atmosphericRefraction.configureCorrectionGrid(sensor, minLine, maxLine); + + // Get the grid knots + final int nbPixelGrid = atmosphericRefraction.getComputationParameters().getNbPixelGrid(); + final int nbLineGrid = atmosphericRefraction.getComputationParameters().getNbLineGrid(); + final double[] pixelGrid = atmosphericRefraction.getComputationParameters().getUgrid(); + final double[] lineGrid = atmosphericRefraction.getComputationParameters().getVgrid(); + + // Computation, for the sensor grid, of the direct location WITH atmospheric refraction + // (full computation) + atmosphericRefraction.reactivateComputation(); + final GeodeticPoint[][] geodeticGridWithAtmosphere = computeDirectLocOnGridWithAtmosphere(pixelGrid, lineGrid, sensor); + // pixelGrid and lineGrid are the knots where the direct loc is computed WITH atmosphere + + // Computation of the inverse location WITHOUT atmospheric refraction for the grid knots + atmosphericRefraction.deactivateComputation(); + final SensorPixel[][] sensorPixelGridInverseWithout = computeInverseLocOnGridWithoutAtmosphere(geodeticGridWithAtmosphere, + nbPixelGrid, nbLineGrid, sensor, minLine, maxLine); + atmosphericRefraction.reactivateComputation(); + + // Compute the grid correction functions (for pixel and line) + atmosphericRefraction.computeGridCorrectionFunctions(sensorPixelGridInverseWithout); + } + + // Fixed point method + // ================== + // Initialization + // -------------- + // compute the sensor pixel on the desired ground point WITHOUT atmosphere + atmosphericRefraction.deactivateComputation(); + final SensorPixel sp0 = inverseLocation(sensorName, point, minLine, maxLine); + if (sp0 == null) { + // Impossible to find the point in the given min line and max line + throw new RuggedException(RuggedMessages.INVALID_RANGE_FOR_LINES, minLine, maxLine, ""); + } + atmosphericRefraction.reactivateComputation(); + + // set up the starting point of the fixed point method + final double pixel0 = sp0.getPixelNumber(); + final double line0 = sp0.getLineNumber(); + + // Apply fixed point method until convergence in pixel and line + // ------------------------------------------------------------ + // compute the first (pixel, line) value: + // initial sensor pixel value + correction due to atmosphere at this same sensor pixel + double corrPixelPrevious = pixel0 + atmosphericRefraction.getBifPixel().value(pixel0, line0); + double corrLinePrevious = line0 + atmosphericRefraction.getBifLine().value(pixel0, line0); + + double deltaCorrPixel = Double.POSITIVE_INFINITY; + double deltaCorrLine = Double.POSITIVE_INFINITY; + + while (deltaCorrPixel > PIXEL_CV_THRESHOLD && deltaCorrLine > LINE_CV_THRESHOLD) { + // Compute the current (pixel, line) value = + // initial sensor pixel value + correction due to atmosphere on the previous sensor pixel + final double corrPixelCurrent = pixel0 + atmosphericRefraction.getBifPixel().value(corrPixelPrevious, corrLinePrevious); + final double corrLineCurrent = line0 + atmosphericRefraction.getBifLine().value(corrPixelPrevious, corrLinePrevious); + + // Compute the delta in pixel and line to check the convergence + deltaCorrPixel = FastMath.abs(corrPixelCurrent - corrPixelPrevious); + deltaCorrLine = FastMath.abs(corrLineCurrent - corrLinePrevious); + + // Store the (pixel, line) for next loop + corrPixelPrevious = corrPixelCurrent; + corrLinePrevious = corrLineCurrent; + } + // The sensor pixel is found ! + final SensorPixel sensorPixelWithAtmosphere = new SensorPixel(corrLinePrevious, corrPixelPrevious); + return sensorPixelWithAtmosphere; + } + + /** Compute the inverse location WITHOUT atmospheric refraction for the geodetic points + * associated to the sensor grid knots. + * @param groundGridWithAtmosphere ground grid found for sensor grid knots with atmosphere + * @param nbPixelGrid size of the pixel grid + * @param nbLineGrid size of the line grid + * @param sensor the line sensor + * @param minLine minimum line number where the search will be performed + * @param maxLine maximum line number where the search will be performed + * @return the sensor pixel grid computed without atmosphere + * @throws RuggedException if invalid range for lines + */ + private SensorPixel[][] computeInverseLocOnGridWithoutAtmosphere(final GeodeticPoint[][] groundGridWithAtmosphere, + final int nbPixelGrid, final int nbLineGrid, + final LineSensor sensor, final int minLine, final int maxLine) + throws RuggedException { + + final SensorPixel[][] sensorPixelGrid = new SensorPixel[nbPixelGrid][nbLineGrid]; + final String sensorName = sensor.getName(); + + for (int uIndex = 0; uIndex < nbPixelGrid; uIndex++) { + for (int vIndex = 0; vIndex < nbLineGrid; vIndex++) { + // Check if the geodetic point exists + if (groundGridWithAtmosphere[uIndex][vIndex] != null) { + final GeodeticPoint groundPoint = groundGridWithAtmosphere[uIndex][vIndex]; + final double currentLat = groundPoint.getLatitude(); + final double currentLon = groundPoint.getLongitude(); + + try { + sensorPixelGrid[uIndex][vIndex] = inverseLocation(sensorName, currentLat, currentLon, minLine, maxLine); + + // Check if the pixel is inside the sensor + if (sensorPixelGrid[uIndex][vIndex] != null && + (sensorPixelGrid[uIndex][vIndex].getPixelNumber() < (-INVLOC_MARGIN) || + sensorPixelGrid[uIndex][vIndex].getPixelNumber() > (INVLOC_MARGIN + sensor.getNbPixels() - 1)) ) { + // Impossible to find the point in the given min line and max line + throw new RuggedException(RuggedMessages.INVALID_RANGE_FOR_LINES, minLine, maxLine, ""); + } else if (sensorPixelGrid[uIndex][vIndex] == null) { + throw new RuggedException(RuggedMessages.INVALID_RANGE_FOR_LINES, minLine, maxLine, ""); + } + } catch (RuggedException re) { + throw new RuggedException(RuggedMessages.INVALID_RANGE_FOR_LINES, minLine, maxLine, ""); + } + + } else { // groundGrid[uIndex][vIndex] == null: impossible to compute inverse loc because ground point not defined + sensorPixelGrid[uIndex][vIndex] = null; + } // groundGrid[uIndex][vIndex] != null + } // end loop vIndex + } // end loop uIndex + return sensorPixelGrid; + } + + /** Computation, for the sensor pixels grid, of the direct location WITH atmospheric refraction. + * (full computation) + * @param pixelGrid the pixel grid + * @param lineGrid the line grid + * @param sensor the line sensor + * @return the ground grid computed with atmosphere + * @throws RuggedException if a problem occurs while computing direct location on sensor grid with atmospheric refraction + */ + private GeodeticPoint[][] computeDirectLocOnGridWithAtmosphere(final double[] pixelGrid, final double[] lineGrid, final LineSensor sensor) + throws RuggedException { + + final int nbPixelGrid = pixelGrid.length; + final int nbLineGrid = lineGrid.length; + final GeodeticPoint[][] groundGridWithAtmosphere = new GeodeticPoint[nbPixelGrid][nbLineGrid]; + final Vector3D sensorPosition = sensor.getPosition(); + + for (int uIndex = 0; uIndex < nbPixelGrid; uIndex++) { + final double pixelNumber = pixelGrid[uIndex]; + for (int vIndex = 0; vIndex < nbLineGrid; vIndex++) { + final double lineNumber = lineGrid[vIndex]; + final AbsoluteDate date = sensor.getDate(lineNumber); + final Vector3D los = sensor.getLOS(date, pixelNumber); + try { + groundGridWithAtmosphere[uIndex][vIndex] = directLocation(date, sensorPosition, los); + } catch (RuggedException re) { // This should never happen + throw new RuggedException(RuggedMessages.INTERNAL_ERROR, re); + } + } // end loop vIndex + } // end loop uIndex + return groundGridWithAtmosphere; } /** Compute distances between two line sensors. @@ -783,7 +1064,7 @@ public class Rugged { private SensorMeanPlaneCrossing getPlaneCrossing(final String sensorName, final int minLine, final int maxLine) throws RuggedException { - + final LineSensor sensor = getLineSensor(sensorName); SensorMeanPlaneCrossing planeCrossing = finders.get(sensorName); if (planeCrossing == null || @@ -802,7 +1083,6 @@ public class Rugged { } return planeCrossing; - } /** Set the mean plane crossing finder for a sensor. @@ -891,7 +1171,6 @@ public class Rugged { return new DerivativeStructure[] { fixedLine, fixedPixel }; - } /** Get transform from spacecraft to inertial frame. @@ -944,5 +1223,4 @@ public class Rugged { public SpacecraftToObservedBody getScToBody() { return scToBody; } - } diff --git a/src/main/java/org/orekit/rugged/errors/Dump.java b/src/main/java/org/orekit/rugged/errors/Dump.java index 13a7b36827913854a7d24045224c9838650cd329..a7f29c7a70c56b9d91441e200a5f89e16b8e9d32 100644 --- a/src/main/java/org/orekit/rugged/errors/Dump.java +++ b/src/main/java/org/orekit/rugged/errors/Dump.java @@ -145,21 +145,21 @@ class Dump { /** Dump a direct location computation. * @param date date of the location - * @param position pixel position in spacecraft frame + * @param sensorPosition sensor position in spacecraft frame * @param los normalized line-of-sight in spacecraft frame * @param lightTimeCorrection flag for light time correction * @param aberrationOfLightCorrection flag for aberration of light correction * @param refractionCorrection flag for refraction correction * @exception RuggedException if date cannot be converted to UTC */ - public void dumpDirectLocation(final AbsoluteDate date, final Vector3D position, final Vector3D los, + public void dumpDirectLocation(final AbsoluteDate date, final Vector3D sensorPosition, final Vector3D los, final boolean lightTimeCorrection, final boolean aberrationOfLightCorrection, final boolean refractionCorrection) throws RuggedException { writer.format(Locale.US, "direct location: date %s position %22.15e %22.15e %22.15e los %22.15e %22.15e %22.15e lightTime %b aberration %b refraction %b %n", convertDate(date), - position.getX(), position.getY(), position.getZ(), + sensorPosition.getX(), sensorPosition.getY(), sensorPosition.getZ(), los.getX(), los.getY(), los.getZ(), lightTimeCorrection, aberrationOfLightCorrection, refractionCorrection); } diff --git a/src/main/java/org/orekit/rugged/errors/DumpManager.java b/src/main/java/org/orekit/rugged/errors/DumpManager.java index 8c466f590afcde7cc5406a2896a53b7539e1f293..eefb9fad3033dd2d0a9416bbfe9419b41913bd3c 100644 --- a/src/main/java/org/orekit/rugged/errors/DumpManager.java +++ b/src/main/java/org/orekit/rugged/errors/DumpManager.java @@ -135,19 +135,19 @@ public class DumpManager { /** Dump a direct location computation. * @param date date of the location - * @param position pixel position in spacecraft frame + * @param sensorPosition sensor position in spacecraft frame * @param los normalized line-of-sight in spacecraft frame * @param lightTimeCorrection flag for light time correction * @param aberrationOfLightCorrection flag for aberration of light correction * @param refractionCorrection flag for refraction correction * @exception RuggedException if date cannot be converted to UTC */ - public static void dumpDirectLocation(final AbsoluteDate date, final Vector3D position, final Vector3D los, + public static void dumpDirectLocation(final AbsoluteDate date, final Vector3D sensorPosition, final Vector3D los, final boolean lightTimeCorrection, final boolean aberrationOfLightCorrection, final boolean refractionCorrection) throws RuggedException { if (isActive()) { - DUMP.get().dumpDirectLocation(date, position, los, lightTimeCorrection, aberrationOfLightCorrection, + DUMP.get().dumpDirectLocation(date, sensorPosition, los, lightTimeCorrection, aberrationOfLightCorrection, refractionCorrection); } } diff --git a/src/main/java/org/orekit/rugged/errors/RuggedMessages.java b/src/main/java/org/orekit/rugged/errors/RuggedMessages.java index 31b763e9a5c09c7c37d510e46a69b70c80866f2b..da289541206c49a5f409cd864ef8fd965996c4ba 100644 --- a/src/main/java/org/orekit/rugged/errors/RuggedMessages.java +++ b/src/main/java/org/orekit/rugged/errors/RuggedMessages.java @@ -80,7 +80,9 @@ public enum RuggedMessages implements Localizable { DUPLICATED_PARAMETER_NAME("a different parameter with name {0} already exists"), INVALID_RUGGED_NAME("invalid rugged name"), UNSUPPORTED_REFINING_CONTEXT("refining using {0} rugged instance is not handled"), - NO_LAYER_DATA("no atmospheric layer data at altitude {0} (lowest altitude: {1})"); + NO_LAYER_DATA("no atmospheric layer data at altitude {0} (lowest altitude: {1})"), + INVALID_STEP("step {0} is not valid : {1}"), + INVALID_RANGE_FOR_LINES("range between min line {0} and max line {1} too small {2}"); // CHECKSTYLE: resume JavadocVariable check diff --git a/src/main/java/org/orekit/rugged/linesensor/SensorPixel.java b/src/main/java/org/orekit/rugged/linesensor/SensorPixel.java index 3f5248cd7705ab2838cb942bf2b7f2555f8772cf..748a71ac782078405e9d3d00b72db44d71b123e8 100644 --- a/src/main/java/org/orekit/rugged/linesensor/SensorPixel.java +++ b/src/main/java/org/orekit/rugged/linesensor/SensorPixel.java @@ -17,9 +17,6 @@ package org.orekit.rugged.linesensor; import java.io.Serializable; -import java.text.NumberFormat; - -import org.hipparchus.util.CompositeFormat; /** Container for sensor pixel. * <p> @@ -63,15 +60,4 @@ public class SensorPixel implements Serializable { public double getPixelNumber() { return pixelNumber; } - - @Override - public String toString() { - final NumberFormat format = CompositeFormat.getDefaultNumberFormat(); - format.setMaximumFractionDigits(6); - return "{Line: " + - format.format(getLineNumber()) + - " , Pixel: " + - format.format(getPixelNumber()) + - "}"; - } } diff --git a/src/main/java/org/orekit/rugged/los/PixelLOS.java b/src/main/java/org/orekit/rugged/los/PixelLOS.java new file mode 100644 index 0000000000000000000000000000000000000000..0bed82a88f7bc92f5ffd41f505a7e670cae6be1d --- /dev/null +++ b/src/main/java/org/orekit/rugged/los/PixelLOS.java @@ -0,0 +1,62 @@ +/* Copyright 2013-2018 CS Systèmes d'Information + * Licensed to CS Systèmes d'Information (CS) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * CS licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.orekit.rugged.los; + +import java.io.Serializable; + +import org.hipparchus.geometry.euclidean.threed.Vector3D; +import org.orekit.rugged.linesensor.SensorPixel; + +/** Container for pixel line-of-sight. + * @author Guylaine Prat + * @since 3.0 + */ +public class PixelLOS implements Serializable { + + /** Serializable UID. */ + private static final long serialVersionUID = -6674056279573271367L; + + /** Sensor pixel. */ + private final SensorPixel sensorPixel; + + /** Pixel line-of-sight in spacecraft frame. */ + private final Vector3D los; + + /** + * Build a new instance. + * @param sensorPixel the sensor pixel cell + * @param los the pixel line-of-sight in spacecraft frame + */ + public PixelLOS(final SensorPixel sensorPixel, final Vector3D los) { + this.sensorPixel = sensorPixel; + this.los = los; + } + + /** + * @return the sensorPixel + */ + public SensorPixel getSensorPixel() { + return sensorPixel; + } + + /** + * @return the lOS in spacecraft frame + */ + public Vector3D getLOS() { + return los; + } +} diff --git a/src/main/java/org/orekit/rugged/refraction/AtmosphericComputationParameters.java b/src/main/java/org/orekit/rugged/refraction/AtmosphericComputationParameters.java new file mode 100644 index 0000000000000000000000000000000000000000..0ca61c4d0aa39c6a6234672abbf39112c2b4deb6 --- /dev/null +++ b/src/main/java/org/orekit/rugged/refraction/AtmosphericComputationParameters.java @@ -0,0 +1,177 @@ +/* Copyright 2013-2018 CS Systèmes d'Information + * Licensed to CS Systèmes d'Information (CS) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * CS licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.orekit.rugged.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 3.0 + */ +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; + + /** 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; + + // 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; + } + + /** 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 + * @throws RuggedException if invalid range for lines + */ + public void configureCorrectionGrid(final LineSensor sensor, final int minLine, final int maxLine) throws RuggedException { + + // 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; + + // CHECKSTYLE: stop UnnecessaryParentheses check + + // 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); + + // CHECKSTYLE: resume UnnecessaryParentheses check + + } + + /** + * 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 + * @throws RuggedException if invalid steps + */ + public void setGridSteps(final int gridPixelStep, final int gridLineStep) throws RuggedException { + + 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; + } + + /** + * @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; + } +} diff --git a/src/main/java/org/orekit/rugged/refraction/AtmosphericRefraction.java b/src/main/java/org/orekit/rugged/refraction/AtmosphericRefraction.java index 8ff9e3e291b4eeb624fdf98802179b6d10b33fcc..961db64f65087e2ef0f5940f8e7c2713b4f67598 100644 --- a/src/main/java/org/orekit/rugged/refraction/AtmosphericRefraction.java +++ b/src/main/java/org/orekit/rugged/refraction/AtmosphericRefraction.java @@ -1,4 +1,4 @@ -/* Copyright 2013-2017 CS Systèmes d'Information +/* Copyright 2013-2018 CS Systèmes d'Information * Licensed to CS Systèmes d'Information (CS) under one or more * contributor license agreements. See the NOTICE file distributed with * this work for additional information regarding copyright ownership. @@ -16,18 +16,57 @@ */ 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; /** * Interface for atmospheric refraction model. * @author Sergio Esteves + * @author Guylaine Prat * @since 2.0 */ -public interface AtmosphericRefraction { +public abstract class AtmosphericRefraction { + + /** Flag to tell if we must compute the correction. + * By default: computation is set up. + * @since 3.0 + */ + private boolean mustBeComputed = true; + + /** Flag to tell if we must compute the correction (for direct location) with an optimization grid. + * By default: optimization is not set up. + * @since 3.0 + */ + private boolean isOptimized = false; + + /** The current atmospheric parameters. + * @since 3.0 + */ + private AtmosphericComputationParameters atmosphericParams; + + /** Bilinear interpolating function for pixel (used by inverse location). + * @since 3.0 + */ + private BilinearInterpolatingFunction bifPixel = null; + + /** Bilinear interpolating function of line (used by inverse location). + * @since 3.0 + */ + private BilinearInterpolatingFunction bifLine = null; + + /** + * Default constructor. + */ + protected AtmosphericRefraction() { + // Set up the atmospheric parameters ... with lazy evaluation of the grid (done only if necessary) + this.atmosphericParams = new AtmosphericComputationParameters(); + } /** Apply correction to the intersected point with an atmospheric refraction model. * @param satPos satellite position, in <em>body frame</em> @@ -37,10 +76,151 @@ public interface AtmosphericRefraction { * @return corrected point with the effect of atmospheric refraction * @throws RuggedException if there is no refraction data at altitude of rawIntersection or see * {@link org.orekit.rugged.utils.ExtendedEllipsoid#pointAtAltitude(Vector3D, Vector3D, double)} or see - * {@link IntersectionAlgorithm#refineIntersection(ExtendedEllipsoid, Vector3D, Vector3D, NormalizedGeodeticPoint)} + * {@link org.orekit.rugged.intersection.IntersectionAlgorithm#refineIntersection(org.orekit.rugged.utils.ExtendedEllipsoid, Vector3D, Vector3D, NormalizedGeodeticPoint)} */ - NormalizedGeodeticPoint applyCorrection(Vector3D satPos, Vector3D satLos, NormalizedGeodeticPoint rawIntersection, + public abstract NormalizedGeodeticPoint applyCorrection(Vector3D satPos, Vector3D satLos, NormalizedGeodeticPoint rawIntersection, IntersectionAlgorithm algorithm) throws RuggedException; + /** Apply correction to the intersected point with an atmospheric refraction model, + * using a time optimized algorithm. + * @param lineSensor the line sensor + * @param sensorPixel the sensor pixel (must be defined) + * @param satPos satellite position, in <em>body frame</em> + * @param satLos sensor pixel 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 + * @throws RuggedException if there is no refraction data at altitude of rawIntersection or see + * {@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)} + * @since 3.0 + */ + public abstract NormalizedGeodeticPoint applyCorrection(LineSensor lineSensor, SensorPixel sensorPixel, + Vector3D satPos, Vector3D satLos, NormalizedGeodeticPoint rawIntersection, + IntersectionAlgorithm algorithm) + throws RuggedException; + + /** Deactivate computation (needed for the inverse location computation). + * @since 3.0 + */ + public void deactivateComputation() { + this.mustBeComputed = false; + } + + /** Reactivate computation (needed for the inverse location computation). + * @since 3.0 + */ + public void reactivateComputation() { + this.mustBeComputed = true; + } + + /** Tell if the computation must be performed. + * @return true if computation must be performed; false otherwise + * @since 3.0 + */ + public boolean mustBeComputed() { + return mustBeComputed; + } + + /** Tell if the computation (for direct location) must be optimized. + * @return true if computation must be optimized; false otherwise + * @since 3.0 + */ + public boolean isOptimized() { + return isOptimized; + } + + /** 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 + * @throws RuggedException if invalid range for lines + */ + public void configureCorrectionGrid(final LineSensor sensor, final int minLine, final int maxLine) throws RuggedException { + + 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 + */ + 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 + */ + 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 + * @throws RuggedException if invalid steps + */ + public void setGridSteps(final int pixelStep, final int lineStep) throws RuggedException { + atmosphericParams.setGridSteps(pixelStep, lineStep); + } + + /** 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 knots (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 + * @throws RuggedException if invalid range for lines + */ + public void computeGridCorrectionFunctions(final SensorPixel[][] sensorPixelGridInverseWithout) throws RuggedException { + + // Compute for a sensor grid, the associated ground grid, WITH atmospheric effect + // (for interpolations we need a regular grid) + // ================================================================================== + final int nbPixelGrid = atmosphericParams.getNbPixelGrid(); + final int nbLineGrid = atmosphericParams.getNbLineGrid(); + final double[] pixelGrid = atmosphericParams.getUgrid(); + final double[] lineGrid = atmosphericParams.getVgrid(); + final double[][] gridDiffPixel = new double[nbPixelGrid][nbLineGrid]; + final double[][] gridDiffLine = new double[nbPixelGrid][nbLineGrid]; + + // Compute the difference between grids knots 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 point in the given min line and max line + throw new RuggedException(RuggedMessages.INVALID_RANGE_FOR_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); + } + + public BilinearInterpolatingFunction getBifPixel() { + return bifPixel; + } + + public BilinearInterpolatingFunction getBifLine() { + return bifLine; + } } diff --git a/src/main/java/org/orekit/rugged/refraction/ConstantRefractionLayer.java b/src/main/java/org/orekit/rugged/refraction/ConstantRefractionLayer.java index cadc13dd9f34acf422777970fdac5aa1f75b8882..5deac72e9b0668252d26bf06a08e212b1dd52470 100644 --- a/src/main/java/org/orekit/rugged/refraction/ConstantRefractionLayer.java +++ b/src/main/java/org/orekit/rugged/refraction/ConstantRefractionLayer.java @@ -1,4 +1,4 @@ -/* Copyright 2013-2017 CS Systèmes d'Information +/* Copyright 2013-2018 CS Systèmes d'Information * Licensed to CS Systèmes d'Information (CS) under one or more * contributor license agreements. See the NOTICE file distributed with * this work for additional information regarding copyright ownership. @@ -19,18 +19,19 @@ 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. */ + /** 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 + * @param lowestAltitude lowest altitude of the layer (m) * @param refractiveIndex refractive index of the layer */ public ConstantRefractionLayer(final double lowestAltitude, final double refractiveIndex) { @@ -38,12 +39,17 @@ public class ConstantRefractionLayer { 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; } - } diff --git a/src/main/java/org/orekit/rugged/refraction/MultiLayerModel.java b/src/main/java/org/orekit/rugged/refraction/MultiLayerModel.java index 1a26f7db50258851524b6d517a06f0ea3a5fa4a2..db8bd800f1b2905fa01c4ca217ab6ad1f5cbeedc 100644 --- a/src/main/java/org/orekit/rugged/refraction/MultiLayerModel.java +++ b/src/main/java/org/orekit/rugged/refraction/MultiLayerModel.java @@ -1,4 +1,4 @@ -/* Copyright 2013-2017 CS Systèmes d'Information +/* Copyright 2013-2018 CS Systèmes d'Information * Licensed to CS Systèmes d'Information (CS) under one or more * contributor license agreements. See the NOTICE file distributed with * this work for additional information regarding copyright ownership. @@ -23,19 +23,22 @@ import java.util.List; import org.hipparchus.geometry.euclidean.threed.Vector3D; import org.hipparchus.util.FastMath; import org.orekit.bodies.GeodeticPoint; -import org.orekit.errors.OrekitException; 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.ExtendedEllipsoid; import org.orekit.rugged.utils.NormalizedGeodeticPoint; /** * Atmospheric refraction model based on multiple layers with associated refractive index. - * @author Sergio Esteves, Luc Maisonobe + * @author Sergio Esteves + * @author Luc Maisonobe + * @author Guylaine Prat * @since 2.0 */ -public class MultiLayerModel implements AtmosphericRefraction { +public class MultiLayerModel extends AtmosphericRefraction { /** Observed body ellipsoid. */ private final ExtendedEllipsoid ellipsoid; @@ -43,7 +46,7 @@ public class MultiLayerModel implements AtmosphericRefraction { /** Constant refraction layers. */ private final List<ConstantRefractionLayer> refractionLayers; - /** Atmosphere lowest altitude. */ + /** Atmosphere lowest altitude (m). */ private final double atmosphereLowestAltitude; /** Simple constructor. @@ -51,166 +54,243 @@ public class MultiLayerModel implements AtmosphericRefraction { * This model uses a built-in set of layers. * </p> * @param ellipsoid the ellipsoid to be used. + * @throws RuggedException if problem at layers construction */ - public MultiLayerModel(final ExtendedEllipsoid ellipsoid) { + public MultiLayerModel(final ExtendedEllipsoid ellipsoid) throws RuggedException { + + super(); + this.ellipsoid = ellipsoid; - refractionLayers = new ArrayList<ConstantRefractionLayer>(15); - refractionLayers.add(new ConstantRefractionLayer(100000.00, 1.000000)); - refractionLayers.add(new ConstantRefractionLayer( 50000.00, 1.000000)); - refractionLayers.add(new ConstantRefractionLayer( 40000.00, 1.000001)); - refractionLayers.add(new ConstantRefractionLayer( 30000.00, 1.000004)); - refractionLayers.add(new ConstantRefractionLayer( 23000.00, 1.000012)); - refractionLayers.add(new ConstantRefractionLayer( 18000.00, 1.000028)); - refractionLayers.add(new ConstantRefractionLayer( 14000.00, 1.000052)); - refractionLayers.add(new ConstantRefractionLayer( 11000.00, 1.000083)); - refractionLayers.add(new ConstantRefractionLayer( 9000.00, 1.000106)); - refractionLayers.add(new ConstantRefractionLayer( 7000.00, 1.000134)); - refractionLayers.add(new ConstantRefractionLayer( 5000.00, 1.000167)); - refractionLayers.add(new ConstantRefractionLayer( 3000.00, 1.000206)); - refractionLayers.add(new ConstantRefractionLayer( 1000.00, 1.000252)); - refractionLayers.add(new ConstantRefractionLayer( 0.00, 1.000278)); - refractionLayers.add(new ConstantRefractionLayer( -1000.00, 1.000306)); - - atmosphereLowestAltitude = refractionLayers.get(refractionLayers.size() - 1).getLowestAltitude(); + this.refractionLayers = new ArrayList<ConstantRefractionLayer>(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). + * @throws RuggedException if problem at layers construction */ - public MultiLayerModel(final ExtendedEllipsoid ellipsoid, final List<ConstantRefractionLayer> refractionLayers) { + public MultiLayerModel(final ExtendedEllipsoid ellipsoid, final List<ConstantRefractionLayer> refractionLayers) + throws RuggedException { + + 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 + * @throws RuggedException if no layer data at asked altitude + * @since 3.0 + */ + private IntersectionLOS computeToLowestAtmosphereLayer( + final Vector3D satPos, final Vector3D satLos, + final NormalizedGeodeticPoint rawIntersection) throws RuggedException { + + 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) + throws RuggedException { + + 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); + } + + /** {@inheritDoc} */ + @Override + public NormalizedGeodeticPoint applyCorrection(final LineSensor lineSensor, final SensorPixel sensorPixel, + final Vector3D satPos, final Vector3D satLos, + final NormalizedGeodeticPoint rawIntersection, + final IntersectionAlgorithm algorithm) throws RuggedException { - try { - if (rawIntersection.getAltitude() < atmosphereLowestAltitude) { - throw new RuggedException(RuggedMessages.NO_LAYER_DATA, rawIntersection.getAltitude(), - atmosphereLowestAltitude); - } + // TODO to be done + throw new RuggedException(RuggedMessages.UNINITIALIZED_CONTEXT, "Atmospheric optimization for direct loc"); + } - Vector3D pos = satPos; - Vector3D los = satLos.normalize(); - double n1 = -1; - GeodeticPoint gp = ellipsoid.transform(satPos, ellipsoid.getBodyFrame(), null); - - 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()); - } - - if (rawIntersection.getAltitude() > refractionLayer.getLowestAltitude()) { - break; - } - - // get intersection point - pos = ellipsoid.pointAtAltitude(pos, los, refractionLayer.getLowestAltitude()); - gp = ellipsoid.transform(pos, ellipsoid.getBodyFrame(), null); - - n1 = n2; +} // end of class MultiLayerModel - } +/** Container for the (position, LOS) of the intersection with the lowest atmospheric layer. + * @author Guylaine Prat + * @since 3.0 + */ +class IntersectionLOS { - return algorithm.refineIntersection(ellipsoid, pos, los, rawIntersection); + /** Position of the intersection with the lowest atmospheric layer. */ + private Vector3D intersectionPos; + /** LOS of the intersection with the lowest atmospheric layer. */ + private Vector3D intersectionLos; - } catch (OrekitException oe) { - throw new RuggedException(oe, oe.getSpecifier(), oe.getParts()); - } + /** 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; } } diff --git a/src/main/java/org/orekit/rugged/utils/GridCreation.java b/src/main/java/org/orekit/rugged/utils/GridCreation.java new file mode 100644 index 0000000000000000000000000000000000000000..d3c687bb9fad7aad27d0fe94663d11c81155d5ad --- /dev/null +++ b/src/main/java/org/orekit/rugged/utils/GridCreation.java @@ -0,0 +1,47 @@ +/* Copyright 2013-2018 CS Systèmes d'Information + * Licensed to CS Systèmes d'Information (CS) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * CS licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.orekit.rugged.utils; + +/** Utility class for grids creation. + * @author Guylaine Prat + * @since 3.0 + */ +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; + } +} diff --git a/src/main/resources/assets/org/orekit/rugged/RuggedMessages_da.utf8 b/src/main/resources/assets/org/orekit/rugged/RuggedMessages_da.utf8 index 8d6ce5dd2039d68e442724e0ce1283cc06a8c934..46205dd2816ea93adc9622222c1f5728d0d8c8a2 100644 --- a/src/main/resources/assets/org/orekit/rugged/RuggedMessages_da.utf8 +++ b/src/main/resources/assets/org/orekit/rugged/RuggedMessages_da.utf8 @@ -88,3 +88,8 @@ UNSUPPORTED_REFINING_CONTEXT = <MISSING TRANSLATION> # no atmospheric layer data at altitude {0} (lowest altitude: {1}) NO_LAYER_DATA = ingen atsmosfæriske lagdata ved højden {0} (laveste højde: {1}) +# step {0} is not valid : {1} +INVALID_STEP = <MISSING TRANSLATION> + +# range between min line {0} and max line {1} too small {2} + INVALID_RANGE_FOR_LINES = <MISSING TRANSLATION> diff --git a/src/main/resources/assets/org/orekit/rugged/RuggedMessages_de.utf8 b/src/main/resources/assets/org/orekit/rugged/RuggedMessages_de.utf8 index 065de378968f59ff30c8930ca6282baadb6bf63b..4d35013672ea41b91d365f0f42a7e068a93e611a 100644 --- a/src/main/resources/assets/org/orekit/rugged/RuggedMessages_de.utf8 +++ b/src/main/resources/assets/org/orekit/rugged/RuggedMessages_de.utf8 @@ -87,3 +87,9 @@ UNSUPPORTED_REFINING_CONTEXT = <MISSING TRANSLATION> # no atmospheric layer data at altitude {0} (lowest altitude: {1}) NO_LAYER_DATA = <MISSING TRANSLATION> + +# step {0} is not valid : {1} +INVALID_STEP = <MISSING TRANSLATION> + +# range between min line {0} and max line {1} too small {2} + INVALID_RANGE_FOR_LINES = <MISSING TRANSLATION> diff --git a/src/main/resources/assets/org/orekit/rugged/RuggedMessages_en.utf8 b/src/main/resources/assets/org/orekit/rugged/RuggedMessages_en.utf8 index fdd3477562e3244a13324262ff8c93fff7f5e795..d40525d9800be390862bceb512f22a810b06b10e 100644 --- a/src/main/resources/assets/org/orekit/rugged/RuggedMessages_en.utf8 +++ b/src/main/resources/assets/org/orekit/rugged/RuggedMessages_en.utf8 @@ -87,3 +87,9 @@ UNSUPPORTED_REFINING_CONTEXT = refining using {0} rugged instance is not handled # no atmospheric layer data at altitude {0} (lowest altitude: {1}) NO_LAYER_DATA = no atmospheric layer data at altitude {0} (lowest altitude: {1}) + +# step {0} is not valid : {1} +INVALID_STEP = step {0} is not valid : {1} + +# range between min line {0} and max line {1} too small {2} + INVALID_RANGE_FOR_LINES = range between min line {0} and max line {1} too small {2} diff --git a/src/main/resources/assets/org/orekit/rugged/RuggedMessages_es.utf8 b/src/main/resources/assets/org/orekit/rugged/RuggedMessages_es.utf8 index 567424c59b99cb68eac1678d4486b7b457f657e9..c8d7a424930cab46f902a28618b8f5809ce27191 100644 --- a/src/main/resources/assets/org/orekit/rugged/RuggedMessages_es.utf8 +++ b/src/main/resources/assets/org/orekit/rugged/RuggedMessages_es.utf8 @@ -86,4 +86,10 @@ INVALID_RUGGED_NAME = el nombre no se considera válido para Rugged UNSUPPORTED_REFINING_CONTEXT = no se puede refinar usando {0} instancias rugged # no atmospheric layer data at altitude {0} (lowest altitude: {1}) -NO_LAYER_DATA = no hay datos de atmósfera para la altitud {0} (altitud mÃnima: {1}) \ No newline at end of file +NO_LAYER_DATA = no hay datos de atmósfera para la altitud {0} (altitud mÃnima: {1}) + +# step {0} is not valid : {1} +INVALID_STEP = <MISSING TRANSLATION> + +# range between min line {0} and max line {1} too small {2} + INVALID_RANGE_FOR_LINES = <MISSING TRANSLATION> diff --git a/src/main/resources/assets/org/orekit/rugged/RuggedMessages_fr.utf8 b/src/main/resources/assets/org/orekit/rugged/RuggedMessages_fr.utf8 index 322d9d94c20dd3879e9bc5fa449be87889b12ea6..f4175bab17953fc918394d39a466e3cd4d6ce4fa 100644 --- a/src/main/resources/assets/org/orekit/rugged/RuggedMessages_fr.utf8 +++ b/src/main/resources/assets/org/orekit/rugged/RuggedMessages_fr.utf8 @@ -87,3 +87,9 @@ UNSUPPORTED_REFINING_CONTEXT = le cas d''affinage utilisant {0} instances de rug # no atmospheric layer data at altitude {0} (lowest altitude: {1}) NO_LAYER_DATA = pas de couche atmosphérique définie pour l''altitude {0} (altitude la plus basse : {1}) + +# step {0} is not valid : {1} +INVALID_STEP = le pas {0} n''est pas valable : {1} + +# range between min line {0} and max line {1} too small {2} + INVALID_RANGE_FOR_LINES = l''écrat entre la ligne min {0} et la ligne max {1} est trop petite {2} diff --git a/src/main/resources/assets/org/orekit/rugged/RuggedMessages_gl.utf8 b/src/main/resources/assets/org/orekit/rugged/RuggedMessages_gl.utf8 index 2ef204f2abe3e887d270911d20aaa7261c5d807a..1773d1b05cf601751da8d1e4afc5876c5676f1c8 100644 --- a/src/main/resources/assets/org/orekit/rugged/RuggedMessages_gl.utf8 +++ b/src/main/resources/assets/org/orekit/rugged/RuggedMessages_gl.utf8 @@ -86,4 +86,10 @@ INVALID_RUGGED_NAME = <MISSING TRANSLATION> UNSUPPORTED_REFINING_CONTEXT = <MISSING TRANSLATION> # no atmospheric layer data at altitude {0} (lowest altitude: {1}) -NO_LAYER_DATA = <MISSING TRANSLATION> \ No newline at end of file +NO_LAYER_DATA = <MISSING TRANSLATION> + +# step {0} is not valid : {1} +INVALID_STEP = <MISSING TRANSLATION> + +# range between min line {0} and max line {1} too small {2} + INVALID_RANGE_FOR_LINES = <MISSING TRANSLATION> diff --git a/src/main/resources/assets/org/orekit/rugged/RuggedMessages_it.utf8 b/src/main/resources/assets/org/orekit/rugged/RuggedMessages_it.utf8 index 87ee2df18e6d2679088b339aa688bde22f6f7f28..ac44c66578ce3802f2c0c3c47f999fc7c5711d02 100644 --- a/src/main/resources/assets/org/orekit/rugged/RuggedMessages_it.utf8 +++ b/src/main/resources/assets/org/orekit/rugged/RuggedMessages_it.utf8 @@ -87,3 +87,9 @@ UNSUPPORTED_REFINING_CONTEXT = il refining che usa {0} istanze rugged non è ges # no atmospheric layer data at altitude {0} (lowest altitude: {1}) NO_LAYER_DATA = nessun dato atmosferico all''altitudine {0} (altitudine più bassa: {1}) + +# step {0} is not valid : {1} +INVALID_STEP = <MISSING TRANSLATION> + +# range between min line {0} and max line {1} too small {2} + INVALID_RANGE_FOR_LINES = <MISSING TRANSLATION> diff --git a/src/main/resources/assets/org/orekit/rugged/RuggedMessages_no.utf8 b/src/main/resources/assets/org/orekit/rugged/RuggedMessages_no.utf8 index c51adee3934b8de6c7551c66101385bcca910c9e..f97635af31d411556fe070a23fd03b82a65352ef 100644 --- a/src/main/resources/assets/org/orekit/rugged/RuggedMessages_no.utf8 +++ b/src/main/resources/assets/org/orekit/rugged/RuggedMessages_no.utf8 @@ -86,4 +86,10 @@ INVALID_RUGGED_NAME = <MISSING TRANSLATION> UNSUPPORTED_REFINING_CONTEXT = <MISSING TRANSLATION> # no atmospheric layer data at altitude {0} (lowest altitude: {1}) -NO_LAYER_DATA = <MISSING TRANSLATION> \ No newline at end of file +NO_LAYER_DATA = <MISSING TRANSLATION> + +# step {0} is not valid : {1} +INVALID_STEP = <MISSING TRANSLATION> + +# range between min line {0} and max line {1} too small {2} + INVALID_RANGE_FOR_LINES = <MISSING TRANSLATION> diff --git a/src/main/resources/assets/org/orekit/rugged/RuggedMessages_ro.utf8 b/src/main/resources/assets/org/orekit/rugged/RuggedMessages_ro.utf8 index 99cf4324d4ddbee53a976227df3cab25a458f637..c3121a6046fe142f0b49ae25996d82bed6092ef8 100644 --- a/src/main/resources/assets/org/orekit/rugged/RuggedMessages_ro.utf8 +++ b/src/main/resources/assets/org/orekit/rugged/RuggedMessages_ro.utf8 @@ -86,4 +86,10 @@ INVALID_RUGGED_NAME = nume invalid pentru Rugged UNSUPPORTED_REFINING_CONTEXT = Estimarea realizată folosind {0} instanÈ›e Rugged nu este suportată # no atmospheric layer data at altitude {0} (lowest altitude: {1}) -NO_LAYER_DATA = nu există informaÈ›ii referitoare la atmosferă pentru altitudinea {0} (altitudinea minimă: {1}) \ No newline at end of file +NO_LAYER_DATA = nu există informaÈ›ii referitoare la atmosferă pentru altitudinea {0} (altitudinea minimă: {1}) + +# step {0} is not valid : {1} +INVALID_STEP = <MISSING TRANSLATION> + +# range between min line {0} and max line {1} too small {2} + INVALID_RANGE_FOR_LINES = <MISSING TRANSLATION> diff --git a/src/site/markdown/building.md b/src/site/markdown/building.md index 074089261a17227bbdf21035c326a3643010850c..c40c95eb2e5d34091cac46b9fbe540d122abef18 100644 --- a/src/site/markdown/building.md +++ b/src/site/markdown/building.md @@ -56,7 +56,7 @@ with the following command: For other commands like generating the site, or generating the [checkstyle](http://checkstyle.sourceforge.net/), -[findbugs](http://findbugs.sourceforge.net/) or +[spotbugs](https://spotbugs.github.io/) or [jacoco](http://www.eclemma.org/jacoco/) reports, see the maven plugins documentation at [maven site](http://maven.apache.org/plugins/index.html). diff --git a/src/site/markdown/guidelines.md b/src/site/markdown/guidelines.md index e9a2ab9db7a1848b275416562891fe0cd4fab527..f5e87394b88e530adf53a38d373c087f09731188 100644 --- a/src/site/markdown/guidelines.md +++ b/src/site/markdown/guidelines.md @@ -41,7 +41,7 @@ The second goal, robustness, has some specific implications for a low level component like Rugged. In some sense, it can be considered an extension of the previous goal as it can also be improved by testing. It can also be improved by automatic checking tools that analyze either source code or binary code. The -[findbugs](http://findbugs.sourceforge.net/) tool is already configured for +[spotbugs](https://spotbugs.github.io/) tool is already configured for automatic checks of the library using a maven plugin. This is however not sufficient. A library is intended to be used by applications @@ -204,9 +204,9 @@ root directory. seek for a line test coverage of at least 80% (more is better) -* *findbugs* (robustness) +* *spotbugs* (robustness) - fix _all_ errors and warnings found by findbugs + fix _all_ errors and warnings found by spotbugs * *no runtime assumptions* (robustness) diff --git a/src/test/java/org/orekit/rugged/api/RuggedBuilderTest.java b/src/test/java/org/orekit/rugged/api/RuggedBuilderTest.java index 3329c05f491702268f1d60a0ea7c7166d27564dd..0a222ec0b68ad13f85b72114d7984fa9e670836c 100644 --- a/src/test/java/org/orekit/rugged/api/RuggedBuilderTest.java +++ b/src/test/java/org/orekit/rugged/api/RuggedBuilderTest.java @@ -22,6 +22,7 @@ import org.hipparchus.geometry.euclidean.threed.RotationConvention; import org.hipparchus.geometry.euclidean.threed.Vector3D; import org.hipparchus.ode.nonstiff.DormandPrince853Integrator; import org.hipparchus.util.FastMath; + import java.io.ByteArrayInputStream; import java.io.ByteArrayOutputStream; import java.io.EOFException; @@ -29,6 +30,7 @@ import java.io.File; import java.io.FileOutputStream; import java.io.IOException; import java.io.StreamCorruptedException; +import java.lang.reflect.Field; import java.net.URISyntaxException; import java.util.ArrayList; import java.util.Arrays; @@ -74,6 +76,10 @@ import org.orekit.rugged.los.TimeDependentLOS; import org.orekit.rugged.raster.RandomLandscapeUpdater; import org.orekit.rugged.raster.TileUpdater; import org.orekit.rugged.raster.VolcanicConeElevationUpdater; +import org.orekit.rugged.refraction.AtmosphericRefraction; +import org.orekit.rugged.refraction.ConstantRefractionLayer; +import org.orekit.rugged.refraction.MultiLayerModel; +import org.orekit.rugged.utils.ExtendedEllipsoid; import org.orekit.time.AbsoluteDate; import org.orekit.time.TimeScale; import org.orekit.time.TimeScalesFactory; @@ -92,7 +98,7 @@ public class RuggedBuilderTest { @Test public void testSetContextWithEphemerides() - throws RuggedException, OrekitException, URISyntaxException { + throws RuggedException, OrekitException, URISyntaxException, NoSuchFieldException, SecurityException, IllegalArgumentException, IllegalAccessException { String path = getClass().getClassLoader().getResource("orekit-data").toURI().getPath(); DataProvidersManager.getInstance().addProvider(new DirectoryCrawler(new File(path))); @@ -195,7 +201,36 @@ public class RuggedBuilderTest { Assert.assertFalse(rugged.isAberrationOfLightCorrected()); Assert.assertFalse(builder.getLightTimeCorrection()); Assert.assertFalse(builder.getAberrationOfLightCorrection()); - + + AtmosphericRefraction atmosphericRefraction = new MultiLayerModel(builder.getEllipsoid()); + builder.setRefractionCorrection(atmosphericRefraction); + rugged = builder.build(); + + MultiLayerModel atmosphericRefractionFromBuilder = (MultiLayerModel) builder.getRefractionCorrection(); + Field atmos = atmosphericRefractionFromBuilder.getClass().getDeclaredField("ellipsoid"); + atmos.setAccessible(true); + ExtendedEllipsoid ellipsoidAtmos = (ExtendedEllipsoid) atmos.get(atmosphericRefractionFromBuilder); + Assert.assertEquals(builder.getEllipsoid().getEquatorialRadius(), ellipsoidAtmos.getEquatorialRadius(), 1.0e-9); + Assert.assertEquals(builder.getEllipsoid().getFlattening(), ellipsoidAtmos.getFlattening(), 1.0e-10); + + Field layers = atmosphericRefractionFromBuilder.getClass().getDeclaredField("refractionLayers"); + layers.setAccessible(true); + @SuppressWarnings("unchecked") + List<ConstantRefractionLayer> layersAtmos = (List<ConstantRefractionLayer>) layers.get(atmosphericRefractionFromBuilder); + + Field layersExpected = atmosphericRefraction.getClass().getDeclaredField("refractionLayers"); + layersExpected.setAccessible(true); + @SuppressWarnings("unchecked") + List<ConstantRefractionLayer> layersAtmosExpected = (List<ConstantRefractionLayer>) layersExpected.get(atmosphericRefraction); + Assert.assertEquals(layersAtmosExpected.size(), layersAtmos.size()); + + List<ConstantRefractionLayer> copyAtmosExpected = new ArrayList<ConstantRefractionLayer>(layersAtmosExpected); + List<ConstantRefractionLayer> copyAtmos = new ArrayList<ConstantRefractionLayer>(layersAtmos); + + Assert.assertTrue(copyAtmosExpected.removeAll(layersAtmos) && copyAtmos.removeAll(layersAtmosExpected)); + Assert.assertTrue(copyAtmosExpected.isEmpty() && copyAtmos.isEmpty()); + + Assert.assertEquals(AlgorithmId.DUVENHAGE, builder.getAlgorithm()); Assert.assertEquals(6378137.0, builder.getEllipsoid().getEquatorialRadius(), 1.0e-9); Assert.assertEquals(1.0 / 298.257222101, builder.getEllipsoid().getFlattening(), 1.0e-10); diff --git a/src/test/java/org/orekit/rugged/api/RuggedTest.java b/src/test/java/org/orekit/rugged/api/RuggedTest.java index 67ded92a05850a76115328ad2a8165da3d736dda..60eb44093f5ac3892358d6313223479cbb575024 100644 --- a/src/test/java/org/orekit/rugged/api/RuggedTest.java +++ b/src/test/java/org/orekit/rugged/api/RuggedTest.java @@ -1,4 +1,4 @@ -/* Copyright 2013-2017 CS Systèmes d'Information +/* Copyright 2013-2018 CS Systèmes d'Information * Licensed to CS Systèmes d'Information (CS) under one or more * contributor license agreements. See the NOTICE file distributed with * this work for additional information regarding copyright ownership. @@ -17,6 +17,9 @@ package org.orekit.rugged.api; +import static org.junit.Assert.assertEquals; +import static org.junit.Assert.fail; + import java.io.File; import java.io.IOException; import java.io.RandomAccessFile; @@ -59,6 +62,8 @@ import org.orekit.orbits.Orbit; import org.orekit.propagation.Propagator; import org.orekit.rugged.TestUtils; import org.orekit.rugged.adjustment.GroundOptimizationProblemBuilder; +import org.orekit.rugged.adjustment.measurements.Observables; +import org.orekit.rugged.adjustment.util.InitInterRefiningTest; import org.orekit.rugged.errors.RuggedException; import org.orekit.rugged.errors.RuggedExceptionWrapper; import org.orekit.rugged.errors.RuggedMessages; @@ -72,9 +77,10 @@ import org.orekit.rugged.los.TimeDependentLOS; import org.orekit.rugged.raster.RandomLandscapeUpdater; import org.orekit.rugged.raster.TileUpdater; import org.orekit.rugged.raster.VolcanicConeElevationUpdater; -import org.orekit.rugged.adjustment.measurements.Observables; -import org.orekit.rugged.adjustment.util.InitInterRefiningTest; +import org.orekit.rugged.refraction.AtmosphericRefraction; +import org.orekit.rugged.refraction.MultiLayerModel; import org.orekit.rugged.utils.DSGenerator; +import org.orekit.rugged.utils.GeodeticUtilities; import org.orekit.time.AbsoluteDate; import org.orekit.time.TimeScale; import org.orekit.time.TimeScalesFactory; @@ -86,7 +92,6 @@ import org.orekit.utils.ParameterDriver; import org.orekit.utils.TimeStampedAngularCoordinates; import org.orekit.utils.TimeStampedPVCoordinates; - public class RuggedTest { @Rule @@ -122,7 +127,8 @@ public class RuggedTest { // los: swath in the (YZ) plane, centered at +Z, ±10° aperture, 960 pixels Vector3D position = new Vector3D(1.5, 0, -0.2); TimeDependentLOS los = TestUtils.createLOSPerfectLine(Vector3D.PLUS_K, Vector3D.PLUS_I, - FastMath.toRadians(10.0), dimension).build(); + FastMath.toRadians(10.0), dimension). + build(); // linear datation model: at reference time we get line 1000, and the rate is one line every 1.5ms LineDatation lineDatation = new LinearLineDatation(crossing, dimension / 2, 1.0 / 1.5e-3); @@ -318,9 +324,119 @@ public class RuggedTest { Assert.assertTrue(Vector3D.distance(pWith, pWithout) > 20.0); Assert.assertTrue(Vector3D.distance(pWith, pWithout) < 20.5); } + } + + @Test + public void testAtmosphericRefractionCorrection() throws RuggedException, OrekitException, URISyntaxException { + String sensorName = "line"; + int dimension = 4000; + + RuggedBuilder builder = initRuggedForAtmosphericTests(dimension, sensorName); + // Build Rugged without atmospheric refraction model + Rugged ruggedWithout = builder.build(); + + // Defines atmospheric refraction model (with the default multi layers model) + AtmosphericRefraction atmosphericRefraction = new MultiLayerModel(ruggedWithout.getEllipsoid()); + int pixelStep = 100; + int lineStep = 100; + atmosphericRefraction.setGridSteps(pixelStep, lineStep); + + // Build Rugged with atmospheric refraction model + builder.setRefractionCorrection(atmosphericRefraction); + Rugged ruggedWith = builder.build(); + + LineSensor lineSensor = ruggedWithout.getLineSensor(sensorName); + int minLine = (int) FastMath.floor(lineSensor.getLine(ruggedWithout.getMinDate())); + int maxLine = (int) FastMath.ceil(lineSensor.getLine(ruggedWithout.getMaxDate())); + + final double pixelThreshold = 1.e-3; + final double lineThreshold = 1.e-2; + + final double epsilonPixel = pixelThreshold; + final double epsilonLine = lineThreshold; + double earthRadius = ruggedWithout.getEllipsoid().getEquatorialRadius(); + + // Direct loc on a line WITHOUT and WITH atmospheric correction + // ============================================================ + double chosenLine = 200.; + GeodeticPoint[] gpWithoutAtmosphericRefractionCorrection = ruggedWithout.directLocation(sensorName, chosenLine); + GeodeticPoint[] gpWithAtmosphericRefractionCorrection = ruggedWith.directLocation(sensorName, chosenLine); + + // Check the shift on the ground due to atmospheric correction + for (int i = 0; i < gpWithAtmosphericRefractionCorrection.length; i++) { + double currentRadius = earthRadius + (gpWithAtmosphericRefractionCorrection[i].getAltitude()+ gpWithoutAtmosphericRefractionCorrection[i].getAltitude())/2.; + double distance = GeodeticUtilities.computeDistanceInMeter(currentRadius, gpWithAtmosphericRefractionCorrection[i], gpWithoutAtmosphericRefractionCorrection[i]); + // Check if the distance is not 0 and < 2m + Assert.assertTrue(distance > 0.0); + Assert.assertTrue(distance < 2.); + } + + // Inverse loc WITH atmospheric correction + // ========================================================================== + for (int i = 0; i < gpWithAtmosphericRefractionCorrection.length; i++) { + + // to check if we go back to the initial point when taking the geodetic point with atmospheric correction + GeodeticPoint gpWith = gpWithAtmosphericRefractionCorrection[i]; + SensorPixel sensorPixelReverseWith = ruggedWith.inverseLocation(sensorName, gpWith, minLine, maxLine); + + if (sensorPixelReverseWith != null) { + assertEquals(i, sensorPixelReverseWith.getPixelNumber(), epsilonPixel); + assertEquals(chosenLine, sensorPixelReverseWith.getLineNumber(), epsilonLine); + } else { + fail("Inverse location failed for pixel " + i + " with atmospheric refraction correction for geodetic point computed with" ); + } + } // end loop on pixel i } + private RuggedBuilder initRuggedForAtmosphericTests(final int dimension, final String sensorName) throws URISyntaxException, RuggedException { + + String path = getClass().getClassLoader().getResource("orekit-data").toURI().getPath(); + DataProvidersManager.getInstance().addProvider(new DirectoryCrawler(new File(path))); + final BodyShape earth = TestUtils.createEarth(); + final Orbit orbit = TestUtils.createOrbit(Constants.EIGEN5C_EARTH_MU); + + // one line sensor + // position: 1.5m in front (+X) and 20 cm above (-Z) of the S/C center of mass + // los: swath in the (YZ) plane, looking at 5° roll, 2.6" per pixel + Vector3D position = new Vector3D(1.5, 0, -0.2); + TimeDependentLOS los = TestUtils.createLOSPerfectLine(new Rotation(Vector3D.PLUS_I, + FastMath.toRadians(5.0), + RotationConvention.VECTOR_OPERATOR).applyTo(Vector3D.PLUS_K), + Vector3D.PLUS_I, + FastMath.toRadians((dimension/2.) * 2.6 / 3600.0), dimension).build(); + + // With the orbit (795km), the size of the pixel on the ground is around : 10m + + // linear datation model: at reference time we get the middle line, and the rate is one line every 1.5ms + AbsoluteDate crossing = new AbsoluteDate("2012-01-01T12:30:00.000", TimeScalesFactory.getUTC()); + LineDatation lineDatation = new LinearLineDatation(crossing, dimension / 2, 1.0 / 1.5e-3); + int firstLine = 0; + int lastLine = dimension; + LineSensor lineSensor = new LineSensor(sensorName, lineDatation, position, los); + AbsoluteDate minDate = lineSensor.getDate(firstLine).shiftedBy(-1.0); + AbsoluteDate maxDate = lineSensor.getDate(lastLine).shiftedBy(+1.0); + + TileUpdater updater = new RandomLandscapeUpdater(800.0, 9000.0, 0.1, 0xf0a401650191f9f6l, FastMath.toRadians(2.0), 257); + + + RuggedBuilder builder = new RuggedBuilder(). + setDigitalElevationModel(updater, 8). + setAlgorithm(AlgorithmId.DUVENHAGE). + setEllipsoid(EllipsoidId.WGS84, BodyRotatingFrameId.ITRF). + setTimeSpan(minDate, maxDate, 0.001, 5.0). + setTrajectory(InertialFrameId.EME2000, + TestUtils.orbitToPV(orbit, earth, minDate.shiftedBy(-1.0), maxDate.shiftedBy(+1.0), 0.25), + 8, CartesianDerivativesFilter.USE_PV, + TestUtils.orbitToQ(orbit, earth, minDate.shiftedBy(-1.0), maxDate.shiftedBy(+1.0), 0.25), + 2, AngularDerivativesFilter.USE_R). + setLightTimeCorrection(false). + setAberrationOfLightCorrection(false). + addLineSensor(lineSensor); + + return builder; + } + @Test public void testFlatBodyCorrection() throws RuggedException, OrekitException, URISyntaxException { @@ -384,7 +500,7 @@ public class RuggedTest { } @Test - public void testLocationsinglePoint() + public void testLocationSinglePoint() throws RuggedException, OrekitException, URISyntaxException { int dimension = 200; @@ -584,12 +700,12 @@ public class RuggedTest { for (int i = 0; i < nbSensors; ++i) { // one line sensor // position: 1.5m in front (+X) and 20 cm above (-Z) of the S/C center of mass - // los: swath in the (YZ) plane, looking roughly at 50° roll (sensor-dependent), 2.6" per pixel + // los: swath in the (YZ) plane, looking roughly at 50° roll (sensor-dependent), 5.2" per pixel Vector3D position = new Vector3D(1.5, 0, -0.2); TimeDependentLOS los = TestUtils.createLOSPerfectLine(new Rotation(Vector3D.PLUS_I, FastMath.toRadians(50.0 - 0.001 * i), RotationConvention.VECTOR_OPERATOR).applyTo(Vector3D.PLUS_K), - Vector3D.PLUS_I, FastMath.toRadians(dimension * 2.6 / 3600.0), dimension).build(); + Vector3D.PLUS_I, FastMath.toRadians((dimension/2.) * 5.2 / 3600.0), dimension).build(); // linear datation model: at reference time we get roughly middle line, and the rate is one line every 1.5ms LineDatation lineDatation = new LinearLineDatation(crossing, i + dimension / 2, 1.0 / 1.5e-3); @@ -917,10 +1033,10 @@ public class RuggedTest { // one line sensor // position: 1.5m in front (+X) and 20 cm above (-Z) of the S/C center of mass - // los: swath in the (YZ) plane, looking at nadir, 2.6" per pixel, 3" sagitta + // los: swath in the (YZ) plane, looking at nadir, 5.2" per pixel, 3" sagitta Vector3D position = new Vector3D(1.5, 0, -0.2); TimeDependentLOS los = TestUtils.createLOSCurvedLine(Vector3D.PLUS_K, Vector3D.PLUS_I, - FastMath.toRadians(dimension * 2.6 / 3600.0), + FastMath.toRadians((dimension/2.) * 5.2 / 3600.0), FastMath.toRadians(3.0 / 3600.0), dimension); @@ -972,13 +1088,14 @@ public class RuggedTest { // one line sensor // position: 1.5m in front (+X) and 20 cm above (-Z) of the S/C center of mass - // los: swath in the (YZ) plane, looking at 50° roll, 2.6" per pixel + // los: swath in the (YZ) plane, looking at 50° roll, 5.2" per pixel Vector3D position = new Vector3D(1.5, 0, -0.2); TimeDependentLOS los = TestUtils.createLOSPerfectLine(new Rotation(Vector3D.PLUS_I, - FastMath.toRadians(50.0), + FastMath.toRadians(5.0), RotationConvention.VECTOR_OPERATOR).applyTo(Vector3D.PLUS_K), Vector3D.PLUS_I, - FastMath.toRadians(dimension * 2.6 / 3600.0), dimension).build(); + FastMath.toRadians((dimension/2.) * 5.2 / 3600.0), dimension).build(); + // In fact the pixel size = 5.2" as we construct the LOS with the full line (dimension) instead of dimension/2 // linear datation model: at reference time we get the middle line, and the rate is one line every 1.5ms LineDatation lineDatation = new LinearLineDatation(crossing, dimension / 2, 1.0 / 1.5e-3); @@ -1017,6 +1134,7 @@ public class RuggedTest { (1 - d) * gp[i].getLatitude() + d * gp[i + 1].getLatitude(), (1 - d) * gp[i].getLongitude() + d * gp[i + 1].getLongitude(), 0, dimension); + Assert.assertEquals(referenceLine, sp.getLineNumber(), maxLineError); Assert.assertEquals(p, sp.getPixelNumber(), maxPixelError); } @@ -1109,17 +1227,18 @@ public class RuggedTest { // one line sensor // position: 1.5m in front (+X) and 20 cm above (-Z) of the S/C center of mass - // los: swath in the (YZ) plane, looking at 50° roll, 2.6" per pixel + // los: swath in the (YZ) plane, looking at 50° roll, 5.2" per pixel Vector3D position = new Vector3D(1.5, 0, -0.2); LOSBuilder losBuilder = TestUtils.createLOSPerfectLine(new Rotation(Vector3D.PLUS_I, FastMath.toRadians(50.0), RotationConvention.VECTOR_OPERATOR).applyTo(Vector3D.PLUS_K), Vector3D.PLUS_I, - FastMath.toRadians(dimension * 2.6 / 3600.0), dimension); + FastMath.toRadians((dimension/2.) * 5.2 / 3600.0), dimension); losBuilder.addTransform(new FixedRotation("roll", Vector3D.MINUS_I, 0.0)); losBuilder.addTransform(new FixedRotation("pitch", Vector3D.MINUS_J, 0.0)); TimeDependentLOS los = losBuilder.build(); + // In fact the pixel size = 5.2" as we construct the LOS with the full line (dimension) instead of dimension/2 // linear datation model: at reference time we get the middle line, and the rate is one line every 1.5ms LineDatation lineDatation = new LinearLineDatation(crossing, dimension / 2, 1.0 / 1.5e-3); @@ -1165,7 +1284,7 @@ public class RuggedTest { java.lang.reflect.Method getGenerator = GroundOptimizationProblemBuilder.class.getSuperclass().getDeclaredMethod("getGenerator"); getGenerator.setAccessible(true); - + DSGenerator generator = (DSGenerator) getGenerator.invoke(optimizationPbBuilder); double referenceLine = 0.87654 * dimension; @@ -1273,13 +1392,14 @@ public class RuggedTest { // one line sensor // position: 1.5m in front (+X) and 20 cm above (-Z) of the S/C center of mass - // los: swath in the (YZ) plane, looking at 50° roll, 2.6" per pixel + // los: swath in the (YZ) plane, looking at 50° roll, 5.2" per pixel Vector3D position = new Vector3D(1.5, 0, -0.2); TimeDependentLOS los = TestUtils.createLOSPerfectLine(new Rotation(Vector3D.PLUS_I, FastMath.toRadians(50.0), RotationConvention.VECTOR_OPERATOR).applyTo(Vector3D.PLUS_K), Vector3D.PLUS_I, - FastMath.toRadians(dimension * 2.6 / 3600.0), dimension).build(); + FastMath.toRadians((dimension/2.) * 5.2 / 3600.0), dimension).build(); + // In fact the pixel size = 5.2" as we construct the LOS with the full line (dimension) instead of dimension/2 // linear datation model: at reference time we get line 100, and the rate is one line every 1.5ms LineDatation lineDatation = new LinearLineDatation(crossing, dimension / 2, 1.0 / 1.5e-3); @@ -1340,22 +1460,23 @@ public class RuggedTest { } private void checkLineDatation(int dimension, double maxLineError) - throws RuggedException, OrekitException, URISyntaxException { + throws RuggedException, OrekitException, URISyntaxException { - String path = getClass().getClassLoader().getResource("orekit-data").toURI().getPath(); - DataProvidersManager.getInstance().addProvider(new DirectoryCrawler(new File(path))); + String path = getClass().getClassLoader().getResource("orekit-data").toURI().getPath(); + DataProvidersManager.getInstance().addProvider(new DirectoryCrawler(new File(path))); - AbsoluteDate crossing = new AbsoluteDate("2012-01-01T12:30:00.000", TimeScalesFactory.getUTC()); + AbsoluteDate crossing = new AbsoluteDate("2012-01-01T12:30:00.000", TimeScalesFactory.getUTC()); - // one line sensor + // one line sensor // position: 1.5m in front (+X) and 20 cm above (-Z) of the S/C center of mass // los: swath in the (YZ) plane, looking at 50° roll, 2.6" per pixel Vector3D position = new Vector3D(1.5, 0, -0.2); TimeDependentLOS los = TestUtils.createLOSPerfectLine(new Rotation(Vector3D.PLUS_I, FastMath.toRadians(50.0), RotationConvention.VECTOR_OPERATOR).applyTo(Vector3D.PLUS_K), - Vector3D.PLUS_I, - FastMath.toRadians(dimension * 2.6 / 3600.0), dimension).build(); + Vector3D.PLUS_I, + FastMath.toRadians(dimension * 2.6 / 3600.0), dimension).build(); + // In fact the pixel size = 5.2" as we construct the LOS with the full line (dimension) instead of dimension/2 // linear datation model: at reference time we get the middle line, and the rate is one line every 1.5ms LineDatation lineDatation = new LinearLineDatation(crossing, dimension / 2, 1.0 / 1.5e-3); @@ -1372,13 +1493,13 @@ public class RuggedTest { Assert.assertEquals(firstLine, recomputedFirstLine, maxLineError); Assert.assertEquals(lastLine, recomputedLastLine, maxLineError); } - + @Test public void testDistanceBetweenLOS() throws RuggedException { InitInterRefiningTest refiningTest = new InitInterRefiningTest(); refiningTest.initRefiningTest(); - + final SensorPixel realPixelA = new SensorPixel(2005.015883575199, 18004.968656395424); final SensorPixel realPixelB = new SensorPixel(4798.487736488162, 13952.2195710654); @@ -1390,7 +1511,7 @@ public class RuggedTest { Assert.assertEquals(expectedDistanceBetweenLOS, distancesBetweenLOS[0], 1.e-8); Assert.assertEquals(expectedDistanceToTheGround, distancesBetweenLOS[1], 1.e-5); } - + @Test public void testDistanceBetweenLOSDerivatives() throws RuggedException, NoSuchMethodException, SecurityException, IllegalAccessException, IllegalArgumentException, InvocationTargetException { @@ -1430,7 +1551,7 @@ public class RuggedTest { } } - + @Before public void setUp() throws OrekitException, URISyntaxException { TestUtils.clearFactories(); diff --git a/src/test/java/org/orekit/rugged/errors/RuggedMessagesTest.java b/src/test/java/org/orekit/rugged/errors/RuggedMessagesTest.java index 9755c416a493521007cf6f6f4c2cb793cc4294bd..913b3d7de41359d0e1239259e5e207ac104857ac 100644 --- a/src/test/java/org/orekit/rugged/errors/RuggedMessagesTest.java +++ b/src/test/java/org/orekit/rugged/errors/RuggedMessagesTest.java @@ -30,7 +30,7 @@ public class RuggedMessagesTest { private final String[] LANGUAGES_LIST = { "da", "de", "en", "es", "fr", "gl", "it", "no", "ro" } ; @Test public void testMessageNumber() { - Assert.assertEquals(30, RuggedMessages.values().length); + Assert.assertEquals(32, RuggedMessages.values().length); } @Test diff --git a/src/test/java/org/orekit/rugged/raster/RandomLandscapeUpdater.java b/src/test/java/org/orekit/rugged/raster/RandomLandscapeUpdater.java index 0cf94ac1935cb1a4c33bcdd79fc7ce15f29017e4..45e927b3448a41d8ccffdd1a6fc11e5f4abfb1c3 100644 --- a/src/test/java/org/orekit/rugged/raster/RandomLandscapeUpdater.java +++ b/src/test/java/org/orekit/rugged/raster/RandomLandscapeUpdater.java @@ -30,9 +30,18 @@ public class RandomLandscapeUpdater implements TileUpdater { private int n; private double[][] heightMap; + /** + * @param baseH elevation of the base of DEM; unit = m + * @param initialScale + * @param reductionFactor for smoothing for low value (0.1) or not (0.5 for instance) the landscape + * @param seed + * @param size size in latitude / size in longitude (rad) + * @param n number of latitude / number of longitude + * @throws MathIllegalArgumentException + */ public RandomLandscapeUpdater(double baseH, double initialScale, double reductionFactor, long seed, double size, int n) - throws MathIllegalArgumentException { + throws MathIllegalArgumentException { if (!ArithmeticUtils.isPowerOfTwo(n - 1)) { throw new MathIllegalArgumentException(LocalizedCoreFormats.SIMPLE_MESSAGE, @@ -100,7 +109,7 @@ public class RandomLandscapeUpdater implements TileUpdater { @Override public void updateTile(double latitude, double longitude, UpdatableTile tile) throws RuggedException { - + double step = size / (n - 1); double minLatitude = size * FastMath.floor(latitude / size); double minLongitude = size * FastMath.floor(longitude / size); @@ -110,7 +119,6 @@ public class RandomLandscapeUpdater implements TileUpdater { tile.setElevation(i, j, heightMap[i][j]); } } - } private double mean(int i1, int j1, int i2, int j2, int i3, int j3, int i4, int j4) { diff --git a/src/test/java/org/orekit/rugged/refraction/AtmosphericRefractionTest.java b/src/test/java/org/orekit/rugged/refraction/AtmosphericRefractionTest.java new file mode 100644 index 0000000000000000000000000000000000000000..63586b030db54fb5d21e49dfe8edcafc6b374a21 --- /dev/null +++ b/src/test/java/org/orekit/rugged/refraction/AtmosphericRefractionTest.java @@ -0,0 +1,66 @@ +package org.orekit.rugged.refraction; + +import static org.junit.Assert.assertFalse; + +import java.util.ArrayList; + +import org.hipparchus.geometry.euclidean.threed.Vector3D; +import org.junit.Assert; +import org.junit.Test; +import org.orekit.rugged.errors.RuggedException; +import org.orekit.rugged.errors.RuggedMessages; +import org.orekit.rugged.linesensor.LineSensor; +import org.orekit.rugged.los.LOSBuilder; +import org.orekit.rugged.los.TimeDependentLOS; + +public class AtmosphericRefractionTest { + + @Test + public void testBadConfig() throws RuggedException { + + int dimension = 400; + + TimeDependentLOS los = new LOSBuilder(new ArrayList<Vector3D>(dimension)).build(); + LineSensor lineSensor = new LineSensor("line", null, Vector3D.ZERO, los); + + // Defines atmospheric refraction model (with the default multi layers model) + AtmosphericRefraction atmosphericRefraction = new MultiLayerModel(null); + + // Check the context + atmosphericRefraction.setGridSteps(100, 100); + atmosphericRefraction.configureCorrectionGrid(lineSensor, 0, 300); + assertFalse(atmosphericRefraction.isSameContext("otherSensor", 0, 300)); + assertFalse(atmosphericRefraction.isSameContext("line", 42, 300)); + assertFalse(atmosphericRefraction.isSameContext("line", 0, 42)); + + // Check the test of validity of min / max line vs line step + try { + atmosphericRefraction.setGridSteps(100, 100); + atmosphericRefraction.configureCorrectionGrid(lineSensor, 0, 100); + Assert.fail("An exception should have been thrown"); + + } catch (RuggedException re) { + Assert.assertEquals(RuggedMessages.INVALID_RANGE_FOR_LINES,re.getSpecifier()); + } + + // Bad pixel step + try { + atmosphericRefraction.setGridSteps(-5, 100); + atmosphericRefraction.configureCorrectionGrid(lineSensor, 0, 100); + Assert.fail("An exception should have been thrown"); + + } catch (RuggedException re) { + Assert.assertEquals(RuggedMessages.INVALID_STEP,re.getSpecifier()); + } + + // Bad line step + try { + atmosphericRefraction.setGridSteps(10, -42); + atmosphericRefraction.configureCorrectionGrid(lineSensor, 0, 100); + Assert.fail("An exception should have been thrown"); + + } catch (RuggedException re) { + Assert.assertEquals(RuggedMessages.INVALID_STEP,re.getSpecifier()); + } + } +} diff --git a/src/test/java/org/orekit/rugged/utils/GeodeticUtilities.java b/src/test/java/org/orekit/rugged/utils/GeodeticUtilities.java new file mode 100644 index 0000000000000000000000000000000000000000..cabacab8ffd0dc73be00624d947e79c371d4bb7f --- /dev/null +++ b/src/test/java/org/orekit/rugged/utils/GeodeticUtilities.java @@ -0,0 +1,270 @@ +/* Copyright 2013-2018 CS Systèmes d'Information + * Licensed to CS Systèmes d'Information (CS) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * CS licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.orekit.rugged.utils; + +import java.text.NumberFormat; + +import org.hipparchus.exception.MathRuntimeException; +import org.hipparchus.geometry.euclidean.threed.Vector3D; +import org.hipparchus.util.CompositeFormat; +import org.hipparchus.util.FastMath; +import org.orekit.bodies.GeodeticPoint; +import org.orekit.errors.OrekitException; + +public class GeodeticUtilities { + + /** Compute an (approximate) geodetic distance in meters between geodetic points (long1, lat1) and (long2, lat2). + * TBN: Orekit does not have such a method + * @param earthRadius Earth radius (m) + * @param long1 longitude of first geodetic point (rad) + * @param lat1 latitude of first geodetic point (rad) + * @param long2 longitude of second geodetic point (rad) + * @param lat2 latitude of second geodetic point (rad) + * @return distance in meters + */ + public static double computeDistanceInMeter(double earthRadius, final double long1, final double lat1, + final double long2, final double lat2) { + + // get vectors on unit sphere from angular coordinates + final Vector3D p1 = new Vector3D(long1, lat1); + final Vector3D p2 = new Vector3D(long2, lat2); + return earthRadius * Vector3D.angle(p1, p2); + } + + /** Compute an (approximate) geodetic distance in meters between two geodetic points + * TBN: Orekit does not have such a method + * @param earthRadius Earth radius (m) + * @param gp1 first geodetic point + * @param gp2 second geodetic point + * @return distance in meters + */ + public static double computeDistanceInMeter(double earthRadius, final GeodeticPoint gp1, final GeodeticPoint gp2) { + + return computeDistanceInMeter(earthRadius, gp1.getLongitude(), gp1.getLatitude(), gp2.getLongitude(), gp2.getLatitude()); + } + + + + public static String toStringDMS(GeodeticPoint gp) { + + final NumberFormat format = CompositeFormat.getDefaultNumberFormat(); + format.setMaximumFractionDigits(1); + DMSangle latDMS = convertLatitudeToDMS(FastMath.toDegrees(gp.getLatitude())); + DMSangle lonDMS = convertLongitudeToDMS(FastMath.toDegrees(gp.getLongitude())); + + String latSign = ""; + if (latDMS.getCardinalPoint() == CardinalDirection.South) latSign = "-"; + String lonSign = ""; + if (lonDMS.getCardinalPoint() == CardinalDirection.West) lonSign = "-"; + + return "{lat: " + latSign + + format.format(latDMS.getDegrees()) + "° " + + format.format(latDMS.getMinutes()) + "' " + + format.format(latDMS.getSeconds()) + "'' " + + " lon: " + lonSign + + format.format(lonDMS.getDegrees()) + "° " + + format.format(lonDMS.getMinutes()) + "' " + + format.format(lonDMS.getSeconds()) + "'' " + + "}"; + } + + /** + * Convert longitude (in decimal degrees) in degrees/minutes/seconds + * @param longitudeInDecimalDegrees + * @return the DMS angle + * @throws S2GeolibException + */ + static DMSangle convertLongitudeToDMS(double longitudeInDecimalDegrees) { + + String cardinalDirection; + // Get the cardinal direction + if (longitudeInDecimalDegrees >= 0.0){ + cardinalDirection= "E"; + } else { + cardinalDirection= "W"; + } + + return convertToDMS(longitudeInDecimalDegrees, cardinalDirection); + + } + + /** + * Convert latitude (in decimal degrees) in degrees/minutes/seconds + * @param latitudeInDecimalDegrees + * @return the DMSangle + * @throws S2GeolibException + */ + public static DMSangle convertLatitudeToDMS(double latitudeInDecimalDegrees){ + + String cardinalDirection; + // Get the cardinal direction + if (latitudeInDecimalDegrees >= 0.0){ + cardinalDirection= "N"; + } else { + cardinalDirection= "S"; + } + + return convertToDMS(latitudeInDecimalDegrees, cardinalDirection); + + } + + /** + * Convert angle (in decimal degrees) in degrees/minutes/seconds and add the associated cardinal direction + * @param angleInDecimalDegrees angle in decimal degrees + * @param cardinalDirection the associated cardinal direction + * @return the DMSangle + * @throws S2GeolibException + */ + private static DMSangle convertToDMS(double angleInDecimalDegrees, String cardinalDirection) throws OrekitException{ + + // We know the cardinal direction so we work on |angleInDecimalDegrees| the positive value + double angleInDD = FastMath.abs(angleInDecimalDegrees); + // Get the degrees part + int degreesPart = (int) FastMath.floor(angleInDD); + + // Get the minutes part (always positive value) + double minutesInDecimal = 60.*(angleInDD - degreesPart); + int minutesPart = (int) FastMath.floor(minutesInDecimal); + + // Get the seconds (in decimal) + double secondsInDecimal = 60.*(minutesInDecimal - minutesPart); + + // Due to rounding problem (at equator around 30 micrometre > diameter of human hair) + if (secondsInDecimal > (60. - 1.e-8)) { + secondsInDecimal = 0.0; + minutesPart++; + if (minutesPart == 60) { + minutesPart = 0; + degreesPart++; + } + } + + return new DMSangle(degreesPart, minutesPart, secondsInDecimal, cardinalDirection); + } +} +class DMSangle { + + /** + * Degrees part of the angle + */ + private int degrees; + + /** + * Minutes part of the angle + */ + private int minutes; + + /** + * Seconds part of the angle + */ + private double seconds; + + /** + * Cardinal direction + */ + private CardinalDirection cardinalPoint; + + /** + * Create a DMS angle for a GeodeticPoint (for instance) + * @param degrees degrees part + * @param minutes minutes part + * @param seconds seconds part + * @param cardinalName cardinal direction + * @throws S2GeolibException + */ + public DMSangle(int degrees, int minutes, double seconds, String cardinalName) throws OrekitException { + this.degrees = degrees; + this.minutes = minutes; + this.seconds = seconds; + this.cardinalPoint = CardinalDirection.getCardinalDirectionFromName(cardinalName); + if (this.cardinalPoint == null){ + throw new OrekitException(new MathRuntimeException(null, this.cardinalPoint)); + } + } + + /** + * @return the degrees part of the angle + */ + public int getDegrees() { + return degrees; + } + /** + * @return the minutes part of the angle + */ + public int getMinutes() { + return minutes; + } + /** + * @return the seconds part of the angle + */ + public double getSeconds() { + return seconds; + } + /** + * @return the cardinal direction of the latitude or the longitude + */ + public CardinalDirection getCardinalPoint() { + return cardinalPoint; + } +} + +enum CardinalDirection { + North("North","N"), + South("South","S"), + West("West","W"), + East("East","E"); + /** + * Cardinal direction full name + */ + private String cardinalFullName = null; + + /** + * Cardinal direction short name + */ + private String cardinalShortName = null; + + /** + * Private constructor + * @param fullName + * @param shortName + */ + private CardinalDirection(String fullName, String shortName){ + this.cardinalFullName = fullName; + this.cardinalShortName = shortName; + } + + /** + * Get the cardinal direction from name (long or short) + * @param cardinalName cardinal name (long or short) + * @return the CardinalDirection + */ + public static CardinalDirection getCardinalDirectionFromName(String cardinalName){ + // Get the cardinal direction from full name ... + for (CardinalDirection currentName : CardinalDirection.values()) { + if (currentName.cardinalFullName.equals(cardinalName)) { + return currentName; + } + } + // otherwise get for short name ... + for (CardinalDirection currentName : CardinalDirection.values()) { + if (currentName.cardinalShortName.equals(cardinalName)) { + return currentName; + } + } + return null; + } +} \ No newline at end of file