Date: Mon, 4 Dec 2023 08:51:27 +0100
Subject: [PATCH 012/359] =?UTF-8?q?Renamed=20SaastamoinenModel=20=E2=86=92?=
=?UTF-8?q?=20ModifiedSaastamoinenModel?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
---
.../EstimatedTroposphericModel.java | 6 +-
.../ModifiedSaastamoinenModel.java | 456 ++++++++++++++++++
.../earth/troposphere/SaastamoinenModel.java | 401 +--------------
.../common/AbstractOrbitDetermination.java | 4 +-
.../measurements/BistaticRangeRateTest.java | 6 +-
.../measurements/BistaticRangeTest.java | 6 +-
.../estimation/measurements/Range2Test.java | 6 +-
.../measurements/RangeAnalyticTest.java | 6 +-
.../measurements/RangeRateTest.java | 6 +-
.../estimation/measurements/RangeTest.java | 6 +-
.../estimation/measurements/TDOATest.java | 6 +-
.../TurnAroundRangeAnalyticTest.java | 6 +-
.../measurements/TurnAroundRangeTest.java | 6 +-
.../modifiers/TropoModifierTest.java | 18 +-
.../ModifiedSaastamoinenModelTest.java | 437 +++++++++++++++++
.../troposphere/SaastamoinenModelTest.java | 2 +-
16 files changed, 960 insertions(+), 418 deletions(-)
create mode 100644 src/main/java/org/orekit/models/earth/troposphere/ModifiedSaastamoinenModel.java
create mode 100644 src/test/java/org/orekit/models/earth/troposphere/ModifiedSaastamoinenModelTest.java
diff --git a/src/main/java/org/orekit/models/earth/troposphere/EstimatedTroposphericModel.java b/src/main/java/org/orekit/models/earth/troposphere/EstimatedTroposphericModel.java
index ece77cc7c0..90749f2e7e 100644
--- a/src/main/java/org/orekit/models/earth/troposphere/EstimatedTroposphericModel.java
+++ b/src/main/java/org/orekit/models/earth/troposphere/EstimatedTroposphericModel.java
@@ -45,7 +45,7 @@
* the {@link GlobalMappingFunctionModel Global Mapping Function}, or
* the {@link NiellMappingFunctionModel Niell Mapping Function}
*
- * The tropospheric zenith delay δh is computed empirically with a {@link SaastamoinenModel}
+ * The tropospheric zenith delay δh is computed empirically with a {@link ModifiedSaastamoinenModel}
* while the tropospheric total zenith delay δt is estimated as a {@link ParameterDriver}
*/
public class EstimatedTroposphericModel implements DiscreteTroposphericModel {
@@ -105,7 +105,7 @@ public List getParametersDrivers() {
public double pathDelay(final double elevation, final GeodeticPoint point,
final double[] parameters, final AbsoluteDate date) {
// Use an empirical model for tropospheric zenith hydro-static delay : Saastamoinen model
- final SaastamoinenModel saastamoinen = new SaastamoinenModel(t0, p0, 0.0);
+ final ModifiedSaastamoinenModel saastamoinen = new ModifiedSaastamoinenModel(t0, p0, 0.0);
// Zenith delays. elevation = pi/2 because we compute the delay in the zenith direction
final double zhd = saastamoinen.pathDelay(0.5 * FastMath.PI, point, parameters, date);
final double ztd = parameters[0];
@@ -120,7 +120,7 @@ public double pathDelay(final double elevation, final GeodeticPoint point,
public > T pathDelay(final T elevation, final FieldGeodeticPoint point,
final T[] parameters, final FieldAbsoluteDate date) {
// Use an empirical model for tropospheric zenith hydro-static delay : Saastamoinen model
- final SaastamoinenModel saastamoinen = new SaastamoinenModel(t0, p0, 0.0);
+ final ModifiedSaastamoinenModel saastamoinen = new ModifiedSaastamoinenModel(t0, p0, 0.0);
// Zenith delays. elevation = pi/2 because we compute the delay in the zenith direction
final T zhd = saastamoinen.pathDelay(elevation.getPi().multiply(0.5), point, parameters, date);
final T ztd = parameters[0];
diff --git a/src/main/java/org/orekit/models/earth/troposphere/ModifiedSaastamoinenModel.java b/src/main/java/org/orekit/models/earth/troposphere/ModifiedSaastamoinenModel.java
new file mode 100644
index 0000000000..ab03d5ac66
--- /dev/null
+++ b/src/main/java/org/orekit/models/earth/troposphere/ModifiedSaastamoinenModel.java
@@ -0,0 +1,456 @@
+/* Copyright 2011-2012 Space Applications Services
+ * Licensed to CS Communication & Systèmes (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.models.earth.troposphere;
+
+import java.util.Collections;
+import java.util.List;
+import java.util.regex.Pattern;
+
+import org.hipparchus.CalculusFieldElement;
+import org.hipparchus.Field;
+import org.hipparchus.analysis.interpolation.BilinearInterpolatingFunction;
+import org.hipparchus.analysis.interpolation.LinearInterpolator;
+import org.hipparchus.analysis.polynomials.PolynomialFunction;
+import org.hipparchus.analysis.polynomials.PolynomialSplineFunction;
+import org.hipparchus.util.FastMath;
+import org.orekit.annotation.DefaultDataContext;
+import org.orekit.bodies.FieldGeodeticPoint;
+import org.orekit.bodies.GeodeticPoint;
+import org.orekit.data.DataContext;
+import org.orekit.data.DataProvidersManager;
+import org.orekit.errors.OrekitException;
+import org.orekit.errors.OrekitMessages;
+import org.orekit.time.AbsoluteDate;
+import org.orekit.time.FieldAbsoluteDate;
+import org.orekit.utils.InterpolationTableLoader;
+import org.orekit.utils.ParameterDriver;
+
+/** The modified Saastamoinen model. Estimates the path delay imposed to
+ * electro-magnetic signals by the troposphere according to the formula:
+ *
+ * δ = 2.277e-3 / cos z * (P + (1255 / T + 0.05) * e - B * tan²
+ * z) + δR
+ *
+ * with the following input data provided to the model:
+ *
+ * - z: zenith angle
+ * - P: atmospheric pressure
+ * - T: temperature
+ * - e: partial pressure of water vapour
+ * - B, δR: correction terms
+ *
+ *
+ * The model supports custom δR correction terms to be read from a
+ * configuration file (saastamoinen-correction.txt) via the
+ * {@link DataProvidersManager}.
+ *
+ * @author Thomas Neidhart
+ * @see "Guochang Xu, GPS - Theory, Algorithms and Applications, Springer, 2007"
+ */
+public class ModifiedSaastamoinenModel implements DiscreteTroposphericModel {
+
+ /** Default file name for δR correction term table. */
+ public static final String DELTA_R_FILE_NAME = "^saastamoinen-correction\\.txt$";
+
+ /** Default lowest acceptable elevation angle [rad]. */
+ public static final double DEFAULT_LOW_ELEVATION_THRESHOLD = 0.05;
+
+ /** First pattern for δR correction term table. */
+ private static final Pattern FIRST_DELTA_R_PATTERN = Pattern.compile("^\\^");
+
+ /** Second pattern for δR correction term table. */
+ private static final Pattern SECOND_DELTA_R_PATTERN = Pattern.compile("\\$$");
+
+ /** X values for the B function. */
+ private static final double[] X_VALUES_FOR_B = {
+ 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0
+ };
+
+ /** E values for the B function. */
+ private static final double[] Y_VALUES_FOR_B = {
+ 1.156, 1.079, 1.006, 0.938, 0.874, 0.813, 0.757, 0.654, 0.563
+ };
+
+ /** Coefficients for the partial pressure of water vapor polynomial. */
+ private static final double[] E_COEFFICIENTS = {
+ -37.2465, 0.213166, -0.000256908
+ };
+
+ /** Interpolation function for the B correction term. */
+ private final PolynomialSplineFunction bFunction;
+
+ /** Polynomial function for the e term. */
+ private final PolynomialFunction eFunction;
+
+ /** Interpolation function for the delta R correction term. */
+ private final BilinearInterpolatingFunction deltaRFunction;
+
+ /** The temperature at the station [K]. */
+ private double t0;
+
+ /** The atmospheric pressure [mbar]. */
+ private double p0;
+
+ /** The humidity [percent]. */
+ private double r0;
+
+ /** Lowest acceptable elevation angle [rad]. */
+ private double lowElevationThreshold;
+
+ /**
+ * Create a new Saastamoinen model for the troposphere using the given environmental
+ * conditions and table from the reference book.
+ *
+ * @param t0 the temperature at the station [K]
+ * @param p0 the atmospheric pressure at the station [mbar]
+ * @param r0 the humidity at the station [fraction] (50% -> 0.5)
+ * @see #ModifiedSaastamoinenModel(double, double, double, String, DataProvidersManager)
+ * @since 10.1
+ */
+ public ModifiedSaastamoinenModel(final double t0, final double p0, final double r0) {
+ this(t0, p0, r0, defaultDeltaR());
+ }
+
+ /** Create a new Saastamoinen model for the troposphere using the given
+ * environmental conditions. This constructor uses the {@link DataContext#getDefault()
+ * default data context} if {@code deltaRFileName != null}.
+ *
+ * @param t0 the temperature at the station [K]
+ * @param p0 the atmospheric pressure at the station [mbar]
+ * @param r0 the humidity at the station [fraction] (50% -> 0.5)
+ * @param deltaRFileName regular expression for filename containing δR
+ * correction term table (typically {@link #DELTA_R_FILE_NAME}), if null
+ * default values from the reference book are used
+ * @since 7.1
+ * @see #ModifiedSaastamoinenModel(double, double, double, String, DataProvidersManager)
+ */
+ @DefaultDataContext
+ public ModifiedSaastamoinenModel(final double t0, final double p0, final double r0,
+ final String deltaRFileName) {
+ this(t0, p0, r0, deltaRFileName,
+ DataContext.getDefault().getDataProvidersManager());
+ }
+
+ /** Create a new Saastamoinen model for the troposphere using the given
+ * environmental conditions. This constructor allows the user to specify the source of
+ * of the δR file.
+ *
+ * @param t0 the temperature at the station [K]
+ * @param p0 the atmospheric pressure at the station [mbar]
+ * @param r0 the humidity at the station [fraction] (50% -> 0.5)
+ * @param deltaRFileName regular expression for filename containing δR
+ * correction term table (typically {@link #DELTA_R_FILE_NAME}), if null
+ * default values from the reference book are used
+ * @param dataProvidersManager provides access to auxiliary data.
+ * @since 10.1
+ */
+ public ModifiedSaastamoinenModel(final double t0,
+ final double p0,
+ final double r0,
+ final String deltaRFileName,
+ final DataProvidersManager dataProvidersManager) {
+ this(t0, p0, r0,
+ deltaRFileName == null ?
+ defaultDeltaR() :
+ loadDeltaR(deltaRFileName, dataProvidersManager));
+ }
+
+ /** Create a new Saastamoinen model.
+ *
+ * @param t0 the temperature at the station [K]
+ * @param p0 the atmospheric pressure at the station [mbar]
+ * @param r0 the humidity at the station [fraction] (50% -> 0.5)
+ * @param deltaR δR correction term function
+ * @since 7.1
+ */
+ private ModifiedSaastamoinenModel(final double t0, final double p0, final double r0,
+ final BilinearInterpolatingFunction deltaR) {
+ checkParameterRangeInclusive("humidity", r0, 0.0, 1.0);
+ this.t0 = t0;
+ this.p0 = p0;
+ this.r0 = r0;
+ this.bFunction = new LinearInterpolator().interpolate(X_VALUES_FOR_B, Y_VALUES_FOR_B);
+ this.eFunction = new PolynomialFunction(E_COEFFICIENTS);
+ this.deltaRFunction = deltaR;
+ this.lowElevationThreshold = DEFAULT_LOW_ELEVATION_THRESHOLD;
+ }
+
+ /** Create a new Saastamoinen model using a standard atmosphere model.
+ *
+ *
+ * - temperature: 18 degree Celsius
+ *
- pressure: 1013.25 mbar
+ *
- humidity: 50%
+ *
+ *
+ * @return a Saastamoinen model with standard environmental values
+ */
+ public static ModifiedSaastamoinenModel getStandardModel() {
+ return new ModifiedSaastamoinenModel(273.16 + 18, 1013.25, 0.5);
+ }
+
+ /** Check if the given parameter is within an acceptable range.
+ * The bounds are inclusive: an exception is raised when either of those conditions are met:
+ *
+ * - The parameter is strictly greater than upperBound
+ * - The parameter is strictly lower than lowerBound
+ *
+ *
+ * In either of these cases, an OrekitException is raised.
+ *
+ * @param parameterName name of the parameter
+ * @param parameter value of the parameter
+ * @param lowerBound lower bound of the acceptable range (inclusive)
+ * @param upperBound upper bound of the acceptable range (inclusive)
+ */
+ private void checkParameterRangeInclusive(final String parameterName, final double parameter,
+ final double lowerBound, final double upperBound) {
+ if (parameter < lowerBound || parameter > upperBound) {
+ throw new OrekitException(OrekitMessages.INVALID_PARAMETER_RANGE, parameterName,
+ parameter, lowerBound, upperBound);
+ }
+ }
+
+ /** {@inheritDoc}
+ *
+ * The Saastamoinen model is not defined for altitudes below 0.0. for continuity
+ * reasons, we use the value for h = 0 when altitude is negative.
+ *
+ *
+ * There are also numerical issues for elevation angles close to zero. For continuity reasons,
+ * elevations lower than a threshold will use the value obtained
+ * for the threshold itself.
+ *
+ * @see #getLowElevationThreshold()
+ * @see #setLowElevationThreshold(double)
+ */
+ @Override
+ public double pathDelay(final double elevation, final GeodeticPoint point,
+ final double[] parameters, final AbsoluteDate date) {
+
+ // there are no data in the model for negative altitudes and altitude bigger than 5000 m
+ // limit the height to a range of [0, 5000] m
+ final double fixedHeight = FastMath.min(FastMath.max(0, point.getAltitude()), 5000);
+
+ // the corrected temperature using a temperature gradient of -6.5 K/km
+ final double T = t0 - 6.5e-3 * fixedHeight;
+ // the corrected pressure
+ final double P = p0 * FastMath.pow(1.0 - 2.26e-5 * fixedHeight, 5.225);
+ // the corrected humidity
+ final double R = r0 * FastMath.exp(-6.396e-4 * fixedHeight);
+
+ // interpolate the b correction term
+ final double B = bFunction.value(fixedHeight / 1e3);
+ // calculate e
+ final double e = R * FastMath.exp(eFunction.value(T));
+
+ // calculate the zenith angle from the elevation
+ final double z = FastMath.abs(0.5 * FastMath.PI - FastMath.max(elevation, lowElevationThreshold));
+
+ // get correction factor
+ final double deltaR = getDeltaR(fixedHeight, z);
+
+ // calculate the path delay in m
+ final double tan = FastMath.tan(z);
+ final double delta = 2.277e-3 / FastMath.cos(z) *
+ (P + (1255d / T + 5e-2) * e - B * tan * tan) + deltaR;
+
+ return delta;
+ }
+
+ /** {@inheritDoc}
+ *
+ * The Saastamoinen model is not defined for altitudes below 0.0. for continuity
+ * reasons, we use the value for h = 0 when altitude is negative.
+ *
+ *
+ * There are also numerical issues for elevation angles close to zero. For continuity reasons,
+ * elevations lower than a threshold will use the value obtained
+ * for the threshold itself.
+ *
+ * @see #getLowElevationThreshold()
+ * @see #setLowElevationThreshold(double)
+ */
+ @Override
+ public > T pathDelay(final T elevation, final FieldGeodeticPoint point,
+ final T[] parameters, final FieldAbsoluteDate date) {
+
+ final Field field = date.getField();
+ final T zero = field.getZero();
+ // there are no data in the model for negative altitudes and altitude bigger than 5000 m
+ // limit the height to a range of [0, 5000] m
+ final T fixedHeight = FastMath.min(FastMath.max(zero, point.getAltitude()), zero.add(5000));
+
+ // the corrected temperature using a temperature gradient of -6.5 K/km
+ final T T = fixedHeight.multiply(6.5e-3).negate().add(t0);
+ // the corrected pressure
+ final T P = fixedHeight.multiply(2.26e-5).negate().add(1.0).pow(5.225).multiply(p0);
+ // the corrected humidity
+ final T R = FastMath.exp(fixedHeight.multiply(-6.396e-4)).multiply(r0);
+
+ // interpolate the b correction term
+ final T B = bFunction.value(fixedHeight.divide(1e3));
+ // calculate e
+ final T e = R.multiply(FastMath.exp(eFunction.value(T)));
+
+ // calculate the zenith angle from the elevation
+ final T z = FastMath.abs(FastMath.max(elevation, zero.add(lowElevationThreshold)).negate().add(zero.getPi().multiply(0.5)));
+
+ // get correction factor
+ final T deltaR = getDeltaR(fixedHeight, z, field);
+
+ // calculate the path delay in m
+ final T tan = FastMath.tan(z);
+ final T delta = FastMath.cos(z).divide(2.277e-3).reciprocal().
+ multiply(P.add(T.divide(1255d).reciprocal().add(5e-2).multiply(e)).subtract(B.multiply(tan).multiply(tan))).add(deltaR);
+
+ return delta;
+ }
+
+ /** Calculates the delta R correction term using linear interpolation.
+ * @param height the height of the station in m
+ * @param zenith the zenith angle of the satellite
+ * @return the delta R correction term in m
+ */
+ private double getDeltaR(final double height, final double zenith) {
+ // limit the height to a range of [0, 5000] m
+ final double h = FastMath.min(FastMath.max(0, height), 5000);
+ // limit the zenith angle to 90 degree
+ // Note: the function is symmetric for negative zenith angles
+ final double z = FastMath.min(Math.abs(zenith), 0.5 * FastMath.PI);
+ return deltaRFunction.value(h, z);
+ }
+
+ /** Calculates the delta R correction term using linear interpolation.
+ * @param type of the elements
+ * @param height the height of the station in m
+ * @param zenith the zenith angle of the satellite
+ * @param field field used by default
+ * @return the delta R correction term in m
+ */
+ private > T getDeltaR(final T height, final T zenith,
+ final Field field) {
+ final T zero = field.getZero();
+ // limit the height to a range of [0, 5000] m
+ final T h = FastMath.min(FastMath.max(zero, height), zero.add(5000));
+ // limit the zenith angle to 90 degree
+ // Note: the function is symmetric for negative zenith angles
+ final T z = FastMath.min(zenith.abs(), zero.getPi().multiply(0.5));
+ return deltaRFunction.value(h, z);
+ }
+
+ /** Load δR function.
+ * @param deltaRFileName regular expression for filename containing δR
+ * correction term table
+ * @param dataProvidersManager provides access to auxiliary data.
+ * @return δR function
+ */
+ private static BilinearInterpolatingFunction loadDeltaR(
+ final String deltaRFileName,
+ final DataProvidersManager dataProvidersManager) {
+
+ // read the δR interpolation function from the config file
+ final InterpolationTableLoader loader = new InterpolationTableLoader();
+ dataProvidersManager.feed(deltaRFileName, loader);
+ if (!loader.stillAcceptsData()) {
+ final double[] elevations = loader.getOrdinateGrid();
+ for (int i = 0; i < elevations.length; ++i) {
+ elevations[i] = FastMath.toRadians(elevations[i]);
+ }
+ return new BilinearInterpolatingFunction(loader.getAbscissaGrid(), elevations,
+ loader.getValuesSamples());
+ }
+ throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_FILE,
+ SECOND_DELTA_R_PATTERN.
+ matcher(FIRST_DELTA_R_PATTERN.matcher(deltaRFileName).replaceAll("")).
+ replaceAll(""));
+ }
+
+ /** Create the default δR function.
+ * @return δR function
+ */
+ private static BilinearInterpolatingFunction defaultDeltaR() {
+
+ // the correction table in the referenced book only contains values for an angle of 60 - 80
+ // degree, thus for 0 degree, the correction term is assumed to be 0, for degrees > 80 it
+ // is assumed to be the same value as for 80.
+
+ // the height in m
+ final double[] xValForR = {
+ 0, 500, 1000, 1500, 2000, 3000, 4000, 5000
+ };
+
+ // the zenith angle
+ final double[] yValForR = {
+ FastMath.toRadians( 0.00), FastMath.toRadians(60.00), FastMath.toRadians(66.00), FastMath.toRadians(70.00),
+ FastMath.toRadians(73.00), FastMath.toRadians(75.00), FastMath.toRadians(76.00), FastMath.toRadians(77.00),
+ FastMath.toRadians(78.00), FastMath.toRadians(78.50), FastMath.toRadians(79.00), FastMath.toRadians(79.50),
+ FastMath.toRadians(79.75), FastMath.toRadians(80.00), FastMath.toRadians(90.00)
+ };
+
+ final double[][] fval = new double[][] {
+ {
+ 0.000, 0.003, 0.006, 0.012, 0.020, 0.031, 0.039, 0.050, 0.065, 0.075, 0.087, 0.102, 0.111, 0.121, 0.121
+ }, {
+ 0.000, 0.003, 0.006, 0.011, 0.018, 0.028, 0.035, 0.045, 0.059, 0.068, 0.079, 0.093, 0.101, 0.110, 0.110
+ }, {
+ 0.000, 0.002, 0.005, 0.010, 0.017, 0.025, 0.032, 0.041, 0.054, 0.062, 0.072, 0.085, 0.092, 0.100, 0.100
+ }, {
+ 0.000, 0.002, 0.005, 0.009, 0.015, 0.023, 0.029, 0.037, 0.049, 0.056, 0.065, 0.077, 0.083, 0.091, 0.091
+ }, {
+ 0.000, 0.002, 0.004, 0.008, 0.013, 0.021, 0.026, 0.033, 0.044, 0.051, 0.059, 0.070, 0.076, 0.083, 0.083
+ }, {
+ 0.000, 0.002, 0.003, 0.006, 0.011, 0.017, 0.021, 0.027, 0.036, 0.042, 0.049, 0.058, 0.063, 0.068, 0.068
+ }, {
+ 0.000, 0.001, 0.003, 0.005, 0.009, 0.014, 0.017, 0.022, 0.030, 0.034, 0.040, 0.047, 0.052, 0.056, 0.056
+ }, {
+ 0.000, 0.001, 0.002, 0.004, 0.007, 0.011, 0.014, 0.018, 0.024, 0.028, 0.033, 0.039, 0.043, 0.047, 0.047
+ }
+ };
+
+ // the actual delta R is interpolated using a a bilinear interpolator
+ return new BilinearInterpolatingFunction(xValForR, yValForR, fval);
+
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public List getParametersDrivers() {
+ return Collections.emptyList();
+ }
+
+ /** Get the low elevation threshold value for path delay computation.
+ * @return low elevation threshold, in rad.
+ * @see #pathDelay(double, GeodeticPoint, double[], AbsoluteDate)
+ * @see #pathDelay(CalculusFieldElement, FieldGeodeticPoint, CalculusFieldElement[], FieldAbsoluteDate)
+ * @since 10.2
+ */
+ public double getLowElevationThreshold() {
+ return lowElevationThreshold;
+ }
+
+ /** Set the low elevation threshold value for path delay computation.
+ * @param lowElevationThreshold The new value for the threshold [rad]
+ * @see #pathDelay(double, GeodeticPoint, double[], AbsoluteDate)
+ * @see #pathDelay(CalculusFieldElement, FieldGeodeticPoint, CalculusFieldElement[], FieldAbsoluteDate)
+ * @since 10.2
+ */
+ public void setLowElevationThreshold(final double lowElevationThreshold) {
+ this.lowElevationThreshold = lowElevationThreshold;
+ }
+}
+
diff --git a/src/main/java/org/orekit/models/earth/troposphere/SaastamoinenModel.java b/src/main/java/org/orekit/models/earth/troposphere/SaastamoinenModel.java
index 78c9c99883..b4fb6da29e 100644
--- a/src/main/java/org/orekit/models/earth/troposphere/SaastamoinenModel.java
+++ b/src/main/java/org/orekit/models/earth/troposphere/SaastamoinenModel.java
@@ -1,4 +1,4 @@
-/* Copyright 2011-2012 Space Applications Services
+/* Copyright 2023 Thales Alenia Space
* Licensed to CS Communication & Systèmes (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
@@ -16,100 +16,22 @@
*/
package org.orekit.models.earth.troposphere;
-import java.util.Collections;
-import java.util.List;
-import java.util.regex.Pattern;
-
-import org.hipparchus.CalculusFieldElement;
-import org.hipparchus.Field;
-import org.hipparchus.analysis.interpolation.BilinearInterpolatingFunction;
-import org.hipparchus.analysis.interpolation.LinearInterpolator;
-import org.hipparchus.analysis.polynomials.PolynomialFunction;
-import org.hipparchus.analysis.polynomials.PolynomialSplineFunction;
-import org.hipparchus.util.FastMath;
import org.orekit.annotation.DefaultDataContext;
-import org.orekit.bodies.FieldGeodeticPoint;
-import org.orekit.bodies.GeodeticPoint;
import org.orekit.data.DataContext;
import org.orekit.data.DataProvidersManager;
-import org.orekit.errors.OrekitException;
-import org.orekit.errors.OrekitMessages;
-import org.orekit.time.AbsoluteDate;
-import org.orekit.time.FieldAbsoluteDate;
-import org.orekit.utils.InterpolationTableLoader;
-import org.orekit.utils.ParameterDriver;
-/** The modified Saastamoinen model. Estimates the path delay imposed to
- * electro-magnetic signals by the troposphere according to the formula:
- *
- * δ = 2.277e-3 / cos z * (P + (1255 / T + 0.05) * e - B * tan²
- * z) + δR
- *
- * with the following input data provided to the model:
- *
- * - z: zenith angle
- * - P: atmospheric pressure
- * - T: temperature
- * - e: partial pressure of water vapour
- * - B, δR: correction terms
- *
- *
- * The model supports custom δR correction terms to be read from a
- * configuration file (saastamoinen-correction.txt) via the
- * {@link DataProvidersManager}.
- *
- * @author Thomas Neidhart
- * @see "Guochang Xu, GPS - Theory, Algorithms and Applications, Springer, 2007"
+/** The modified Saastamoinen model.
+ * @author Luc Maisonobe
+ * @deprecated as of 12.1, replaced by {@link ModifiedSaastamoinenModel}
*/
-public class SaastamoinenModel implements DiscreteTroposphericModel {
+@Deprecated
+public class SaastamoinenModel extends ModifiedSaastamoinenModel {
/** Default file name for δR correction term table. */
- public static final String DELTA_R_FILE_NAME = "^saastamoinen-correction\\.txt$";
+ public static final String DELTA_R_FILE_NAME = ModifiedSaastamoinenModel.DELTA_R_FILE_NAME;
/** Default lowest acceptable elevation angle [rad]. */
- public static final double DEFAULT_LOW_ELEVATION_THRESHOLD = 0.05;
-
- /** First pattern for δR correction term table. */
- private static final Pattern FIRST_DELTA_R_PATTERN = Pattern.compile("^\\^");
-
- /** Second pattern for δR correction term table. */
- private static final Pattern SECOND_DELTA_R_PATTERN = Pattern.compile("\\$$");
-
- /** X values for the B function. */
- private static final double[] X_VALUES_FOR_B = {
- 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0
- };
-
- /** E values for the B function. */
- private static final double[] Y_VALUES_FOR_B = {
- 1.156, 1.079, 1.006, 0.938, 0.874, 0.813, 0.757, 0.654, 0.563
- };
-
- /** Coefficients for the partial pressure of water vapor polynomial. */
- private static final double[] E_COEFFICIENTS = {
- -37.2465, 0.213166, -0.000256908
- };
-
- /** Interpolation function for the B correction term. */
- private final PolynomialSplineFunction bFunction;
-
- /** Polynomial function for the e term. */
- private final PolynomialFunction eFunction;
-
- /** Interpolation function for the delta R correction term. */
- private final BilinearInterpolatingFunction deltaRFunction;
-
- /** The temperature at the station [K]. */
- private double t0;
-
- /** The atmospheric pressure [mbar]. */
- private double p0;
-
- /** The humidity [percent]. */
- private double r0;
-
- /** Lowest acceptable elevation angle [rad]. */
- private double lowElevationThreshold;
+ public static final double DEFAULT_LOW_ELEVATION_THRESHOLD = ModifiedSaastamoinenModel.DEFAULT_LOW_ELEVATION_THRESHOLD;
/**
* Create a new Saastamoinen model for the troposphere using the given environmental
@@ -118,11 +40,11 @@ public class SaastamoinenModel implements DiscreteTroposphericModel {
* @param t0 the temperature at the station [K]
* @param p0 the atmospheric pressure at the station [mbar]
* @param r0 the humidity at the station [fraction] (50% -> 0.5)
- * @see #SaastamoinenModel(double, double, double, String, DataProvidersManager)
+ * @see #ModifiedSaastamoinenModel(double, double, double, String, DataProvidersManager)
* @since 10.1
*/
public SaastamoinenModel(final double t0, final double p0, final double r0) {
- this(t0, p0, r0, defaultDeltaR());
+ super(t0, p0, r0);
}
/** Create a new Saastamoinen model for the troposphere using the given
@@ -136,13 +58,12 @@ public SaastamoinenModel(final double t0, final double p0, final double r0) {
* correction term table (typically {@link #DELTA_R_FILE_NAME}), if null
* default values from the reference book are used
* @since 7.1
- * @see #SaastamoinenModel(double, double, double, String, DataProvidersManager)
+ * @see #ModifiedSaastamoinenModel(double, double, double, String, DataProvidersManager)
*/
@DefaultDataContext
public SaastamoinenModel(final double t0, final double p0, final double r0,
final String deltaRFileName) {
- this(t0, p0, r0, deltaRFileName,
- DataContext.getDefault().getDataProvidersManager());
+ super(t0, p0, r0, deltaRFileName);
}
/** Create a new Saastamoinen model for the troposphere using the given
@@ -163,294 +84,22 @@ public SaastamoinenModel(final double t0,
final double r0,
final String deltaRFileName,
final DataProvidersManager dataProvidersManager) {
- this(t0, p0, r0,
- deltaRFileName == null ?
- defaultDeltaR() :
- loadDeltaR(deltaRFileName, dataProvidersManager));
- }
-
- /** Create a new Saastamoinen model.
- *
- * @param t0 the temperature at the station [K]
- * @param p0 the atmospheric pressure at the station [mbar]
- * @param r0 the humidity at the station [fraction] (50% -> 0.5)
- * @param deltaR δR correction term function
- * @since 7.1
- */
- private SaastamoinenModel(final double t0, final double p0, final double r0,
- final BilinearInterpolatingFunction deltaR) {
- checkParameterRangeInclusive("humidity", r0, 0.0, 1.0);
- this.t0 = t0;
- this.p0 = p0;
- this.r0 = r0;
- this.bFunction = new LinearInterpolator().interpolate(X_VALUES_FOR_B, Y_VALUES_FOR_B);
- this.eFunction = new PolynomialFunction(E_COEFFICIENTS);
- this.deltaRFunction = deltaR;
- this.lowElevationThreshold = DEFAULT_LOW_ELEVATION_THRESHOLD;
+ super(t0, p0, r0, deltaRFileName, dataProvidersManager);
}
/** Create a new Saastamoinen model using a standard atmosphere model.
- *
- *
- * - temperature: 18 degree Celsius
- *
- pressure: 1013.25 mbar
- *
- humidity: 50%
- *
- *
- * @return a Saastamoinen model with standard environmental values
- */
- public static SaastamoinenModel getStandardModel() {
- return new SaastamoinenModel(273.16 + 18, 1013.25, 0.5);
- }
-
- /** Check if the given parameter is within an acceptable range.
- * The bounds are inclusive: an exception is raised when either of those conditions are met:
- *
- * - The parameter is strictly greater than upperBound
- * - The parameter is strictly lower than lowerBound
- *
- *
- * In either of these cases, an OrekitException is raised.
- *
- * @param parameterName name of the parameter
- * @param parameter value of the parameter
- * @param lowerBound lower bound of the acceptable range (inclusive)
- * @param upperBound upper bound of the acceptable range (inclusive)
- */
- private void checkParameterRangeInclusive(final String parameterName, final double parameter,
- final double lowerBound, final double upperBound) {
- if (parameter < lowerBound || parameter > upperBound) {
- throw new OrekitException(OrekitMessages.INVALID_PARAMETER_RANGE, parameterName,
- parameter, lowerBound, upperBound);
- }
- }
-
- /** {@inheritDoc}
- *
- * The Saastamoinen model is not defined for altitudes below 0.0. for continuity
- * reasons, we use the value for h = 0 when altitude is negative.
- *
- *
- * There are also numerical issues for elevation angles close to zero. For continuity reasons,
- * elevations lower than a threshold will use the value obtained
- * for the threshold itself.
- *
- * @see #getLowElevationThreshold()
- * @see #setLowElevationThreshold(double)
- */
- @Override
- public double pathDelay(final double elevation, final GeodeticPoint point,
- final double[] parameters, final AbsoluteDate date) {
+ *
+ *
+ * - temperature: 18 degree Celsius
+ *
- pressure: 1013.25 mbar
+ *
- humidity: 50%
+ *
+ *
+ * @return a Saastamoinen model with standard environmental values
+ */
+ public static SaastamoinenModel getStandardModel() {
+ return new SaastamoinenModel(273.16 + 18, 1013.25, 0.5);
+ }
- // there are no data in the model for negative altitudes and altitude bigger than 5000 m
- // limit the height to a range of [0, 5000] m
- final double fixedHeight = FastMath.min(FastMath.max(0, point.getAltitude()), 5000);
-
- // the corrected temperature using a temperature gradient of -6.5 K/km
- final double T = t0 - 6.5e-3 * fixedHeight;
- // the corrected pressure
- final double P = p0 * FastMath.pow(1.0 - 2.26e-5 * fixedHeight, 5.225);
- // the corrected humidity
- final double R = r0 * FastMath.exp(-6.396e-4 * fixedHeight);
-
- // interpolate the b correction term
- final double B = bFunction.value(fixedHeight / 1e3);
- // calculate e
- final double e = R * FastMath.exp(eFunction.value(T));
-
- // calculate the zenith angle from the elevation
- final double z = FastMath.abs(0.5 * FastMath.PI - FastMath.max(elevation, lowElevationThreshold));
-
- // get correction factor
- final double deltaR = getDeltaR(fixedHeight, z);
-
- // calculate the path delay in m
- final double tan = FastMath.tan(z);
- final double delta = 2.277e-3 / FastMath.cos(z) *
- (P + (1255d / T + 5e-2) * e - B * tan * tan) + deltaR;
-
- return delta;
- }
-
- /** {@inheritDoc}
- *
- * The Saastamoinen model is not defined for altitudes below 0.0. for continuity
- * reasons, we use the value for h = 0 when altitude is negative.
- *
- *
- * There are also numerical issues for elevation angles close to zero. For continuity reasons,
- * elevations lower than a threshold will use the value obtained
- * for the threshold itself.
- *
- * @see #getLowElevationThreshold()
- * @see #setLowElevationThreshold(double)
- */
- @Override
- public > T pathDelay(final T elevation, final FieldGeodeticPoint point,
- final T[] parameters, final FieldAbsoluteDate date) {
-
- final Field field = date.getField();
- final T zero = field.getZero();
- // there are no data in the model for negative altitudes and altitude bigger than 5000 m
- // limit the height to a range of [0, 5000] m
- final T fixedHeight = FastMath.min(FastMath.max(zero, point.getAltitude()), zero.add(5000));
-
- // the corrected temperature using a temperature gradient of -6.5 K/km
- final T T = fixedHeight.multiply(6.5e-3).negate().add(t0);
- // the corrected pressure
- final T P = fixedHeight.multiply(2.26e-5).negate().add(1.0).pow(5.225).multiply(p0);
- // the corrected humidity
- final T R = FastMath.exp(fixedHeight.multiply(-6.396e-4)).multiply(r0);
-
- // interpolate the b correction term
- final T B = bFunction.value(fixedHeight.divide(1e3));
- // calculate e
- final T e = R.multiply(FastMath.exp(eFunction.value(T)));
-
- // calculate the zenith angle from the elevation
- final T z = FastMath.abs(FastMath.max(elevation, zero.add(lowElevationThreshold)).negate().add(zero.getPi().multiply(0.5)));
-
- // get correction factor
- final T deltaR = getDeltaR(fixedHeight, z, field);
-
- // calculate the path delay in m
- final T tan = FastMath.tan(z);
- final T delta = FastMath.cos(z).divide(2.277e-3).reciprocal().
- multiply(P.add(T.divide(1255d).reciprocal().add(5e-2).multiply(e)).subtract(B.multiply(tan).multiply(tan))).add(deltaR);
-
- return delta;
- }
-
- /** Calculates the delta R correction term using linear interpolation.
- * @param height the height of the station in m
- * @param zenith the zenith angle of the satellite
- * @return the delta R correction term in m
- */
- private double getDeltaR(final double height, final double zenith) {
- // limit the height to a range of [0, 5000] m
- final double h = FastMath.min(FastMath.max(0, height), 5000);
- // limit the zenith angle to 90 degree
- // Note: the function is symmetric for negative zenith angles
- final double z = FastMath.min(Math.abs(zenith), 0.5 * FastMath.PI);
- return deltaRFunction.value(h, z);
- }
-
- /** Calculates the delta R correction term using linear interpolation.
- * @param type of the elements
- * @param height the height of the station in m
- * @param zenith the zenith angle of the satellite
- * @param field field used by default
- * @return the delta R correction term in m
- */
- private > T getDeltaR(final T height, final T zenith,
- final Field field) {
- final T zero = field.getZero();
- // limit the height to a range of [0, 5000] m
- final T h = FastMath.min(FastMath.max(zero, height), zero.add(5000));
- // limit the zenith angle to 90 degree
- // Note: the function is symmetric for negative zenith angles
- final T z = FastMath.min(zenith.abs(), zero.getPi().multiply(0.5));
- return deltaRFunction.value(h, z);
- }
-
- /** Load δR function.
- * @param deltaRFileName regular expression for filename containing δR
- * correction term table
- * @param dataProvidersManager provides access to auxiliary data.
- * @return δR function
- */
- private static BilinearInterpolatingFunction loadDeltaR(
- final String deltaRFileName,
- final DataProvidersManager dataProvidersManager) {
-
- // read the δR interpolation function from the config file
- final InterpolationTableLoader loader = new InterpolationTableLoader();
- dataProvidersManager.feed(deltaRFileName, loader);
- if (!loader.stillAcceptsData()) {
- final double[] elevations = loader.getOrdinateGrid();
- for (int i = 0; i < elevations.length; ++i) {
- elevations[i] = FastMath.toRadians(elevations[i]);
- }
- return new BilinearInterpolatingFunction(loader.getAbscissaGrid(), elevations,
- loader.getValuesSamples());
- }
- throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_FILE,
- SECOND_DELTA_R_PATTERN.
- matcher(FIRST_DELTA_R_PATTERN.matcher(deltaRFileName).replaceAll("")).
- replaceAll(""));
- }
-
- /** Create the default δR function.
- * @return δR function
- */
- private static BilinearInterpolatingFunction defaultDeltaR() {
-
- // the correction table in the referenced book only contains values for an angle of 60 - 80
- // degree, thus for 0 degree, the correction term is assumed to be 0, for degrees > 80 it
- // is assumed to be the same value as for 80.
-
- // the height in m
- final double[] xValForR = {
- 0, 500, 1000, 1500, 2000, 3000, 4000, 5000
- };
-
- // the zenith angle
- final double[] yValForR = {
- FastMath.toRadians( 0.00), FastMath.toRadians(60.00), FastMath.toRadians(66.00), FastMath.toRadians(70.00),
- FastMath.toRadians(73.00), FastMath.toRadians(75.00), FastMath.toRadians(76.00), FastMath.toRadians(77.00),
- FastMath.toRadians(78.00), FastMath.toRadians(78.50), FastMath.toRadians(79.00), FastMath.toRadians(79.50),
- FastMath.toRadians(79.75), FastMath.toRadians(80.00), FastMath.toRadians(90.00)
- };
-
- final double[][] fval = new double[][] {
- {
- 0.000, 0.003, 0.006, 0.012, 0.020, 0.031, 0.039, 0.050, 0.065, 0.075, 0.087, 0.102, 0.111, 0.121, 0.121
- }, {
- 0.000, 0.003, 0.006, 0.011, 0.018, 0.028, 0.035, 0.045, 0.059, 0.068, 0.079, 0.093, 0.101, 0.110, 0.110
- }, {
- 0.000, 0.002, 0.005, 0.010, 0.017, 0.025, 0.032, 0.041, 0.054, 0.062, 0.072, 0.085, 0.092, 0.100, 0.100
- }, {
- 0.000, 0.002, 0.005, 0.009, 0.015, 0.023, 0.029, 0.037, 0.049, 0.056, 0.065, 0.077, 0.083, 0.091, 0.091
- }, {
- 0.000, 0.002, 0.004, 0.008, 0.013, 0.021, 0.026, 0.033, 0.044, 0.051, 0.059, 0.070, 0.076, 0.083, 0.083
- }, {
- 0.000, 0.002, 0.003, 0.006, 0.011, 0.017, 0.021, 0.027, 0.036, 0.042, 0.049, 0.058, 0.063, 0.068, 0.068
- }, {
- 0.000, 0.001, 0.003, 0.005, 0.009, 0.014, 0.017, 0.022, 0.030, 0.034, 0.040, 0.047, 0.052, 0.056, 0.056
- }, {
- 0.000, 0.001, 0.002, 0.004, 0.007, 0.011, 0.014, 0.018, 0.024, 0.028, 0.033, 0.039, 0.043, 0.047, 0.047
- }
- };
-
- // the actual delta R is interpolated using a a bilinear interpolator
- return new BilinearInterpolatingFunction(xValForR, yValForR, fval);
-
- }
-
- /** {@inheritDoc} */
- @Override
- public List getParametersDrivers() {
- return Collections.emptyList();
- }
-
- /** Get the low elevation threshold value for path delay computation.
- * @return low elevation threshold, in rad.
- * @see #pathDelay(double, GeodeticPoint, double[], AbsoluteDate)
- * @see #pathDelay(CalculusFieldElement, FieldGeodeticPoint, CalculusFieldElement[], FieldAbsoluteDate)
- * @since 10.2
- */
- public double getLowElevationThreshold() {
- return lowElevationThreshold;
- }
-
- /** Set the low elevation threshold value for path delay computation.
- * @param lowElevationThreshold The new value for the threshold [rad]
- * @see #pathDelay(double, GeodeticPoint, double[], AbsoluteDate)
- * @see #pathDelay(CalculusFieldElement, FieldGeodeticPoint, CalculusFieldElement[], FieldAbsoluteDate)
- * @since 10.2
- */
- public void setLowElevationThreshold(final double lowElevationThreshold) {
- this.lowElevationThreshold = lowElevationThreshold;
- }
}
diff --git a/src/test/java/org/orekit/estimation/common/AbstractOrbitDetermination.java b/src/test/java/org/orekit/estimation/common/AbstractOrbitDetermination.java
index 828255b177..b84ca47992 100644
--- a/src/test/java/org/orekit/estimation/common/AbstractOrbitDetermination.java
+++ b/src/test/java/org/orekit/estimation/common/AbstractOrbitDetermination.java
@@ -137,7 +137,7 @@
import org.orekit.models.earth.troposphere.MappingFunction;
import org.orekit.models.earth.troposphere.MendesPavlisModel;
import org.orekit.models.earth.troposphere.NiellMappingFunctionModel;
-import org.orekit.models.earth.troposphere.SaastamoinenModel;
+import org.orekit.models.earth.troposphere.ModifiedSaastamoinenModel;
import org.orekit.models.earth.troposphere.TimeSpanEstimatedTroposphericModel;
import org.orekit.models.earth.weather.GlobalPressureTemperatureModel;
import org.orekit.orbits.CartesianOrbit;
@@ -1678,7 +1678,7 @@ private Map createStationsData(final KeyValueFileParser measurement : measurements) {
@@ -324,7 +324,7 @@ public void testParameterDerivativesWithModifier() {
1.0, 3.0, 300.0);
propagator.clearStepHandlers();
- final BistaticRangeRateTroposphericDelayModifier modifier = new BistaticRangeRateTroposphericDelayModifier(SaastamoinenModel.getStandardModel());
+ final BistaticRangeRateTroposphericDelayModifier modifier = new BistaticRangeRateTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
double maxRelativeError = 0;
for (final ObservedMeasurement> measurement : measurements) {
diff --git a/src/test/java/org/orekit/estimation/measurements/BistaticRangeTest.java b/src/test/java/org/orekit/estimation/measurements/BistaticRangeTest.java
index 93602f795d..cc70a99746 100644
--- a/src/test/java/org/orekit/estimation/measurements/BistaticRangeTest.java
+++ b/src/test/java/org/orekit/estimation/measurements/BistaticRangeTest.java
@@ -26,7 +26,7 @@
import org.orekit.estimation.Context;
import org.orekit.estimation.EstimationTestUtils;
import org.orekit.estimation.measurements.modifiers.BistaticRangeTroposphericDelayModifier;
-import org.orekit.models.earth.troposphere.SaastamoinenModel;
+import org.orekit.models.earth.troposphere.ModifiedSaastamoinenModel;
import org.orekit.orbits.OrbitType;
import org.orekit.orbits.PositionAngleType;
import org.orekit.propagation.Propagator;
@@ -169,7 +169,7 @@ public void testStateDerivativesWithModifier() {
1.0, 3.0, 300.0);
propagator.clearStepHandlers();
- final BistaticRangeTroposphericDelayModifier modifier = new BistaticRangeTroposphericDelayModifier(SaastamoinenModel.getStandardModel());
+ final BistaticRangeTroposphericDelayModifier modifier = new BistaticRangeTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
double maxRelativeError = 0;
for (final ObservedMeasurement> measurement : measurements) {
@@ -325,7 +325,7 @@ public void testParameterDerivativesWithModifier() {
1.0, 3.0, 300.0);
propagator.clearStepHandlers();
- final BistaticRangeTroposphericDelayModifier modifier = new BistaticRangeTroposphericDelayModifier(SaastamoinenModel.getStandardModel());
+ final BistaticRangeTroposphericDelayModifier modifier = new BistaticRangeTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
double maxRelativeError = 0;
for (final ObservedMeasurement> measurement : measurements) {
diff --git a/src/test/java/org/orekit/estimation/measurements/Range2Test.java b/src/test/java/org/orekit/estimation/measurements/Range2Test.java
index 64c9502ca8..6b97fb5a91 100644
--- a/src/test/java/org/orekit/estimation/measurements/Range2Test.java
+++ b/src/test/java/org/orekit/estimation/measurements/Range2Test.java
@@ -28,7 +28,7 @@
import org.orekit.estimation.Context;
import org.orekit.estimation.EstimationTestUtils;
import org.orekit.estimation.measurements.modifiers.RangeTroposphericDelayModifier;
-import org.orekit.models.earth.troposphere.SaastamoinenModel;
+import org.orekit.models.earth.troposphere.ModifiedSaastamoinenModel;
import org.orekit.orbits.OrbitType;
import org.orekit.orbits.PositionAngleType;
import org.orekit.propagation.Propagator;
@@ -330,7 +330,7 @@ void genericTestStateDerivatives(final boolean isModifier, final boolean printRe
) {
// Add modifiers if test implies it
- final RangeTroposphericDelayModifier modifier = new RangeTroposphericDelayModifier(SaastamoinenModel.getStandardModel());
+ final RangeTroposphericDelayModifier modifier = new RangeTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
if (isModifier) {
((Range) measurement).addModifier(modifier);
}
@@ -481,7 +481,7 @@ void genericTestParameterDerivatives(final boolean isModifier, final boolean pri
) {
// Add modifiers if test implies it
- final RangeTroposphericDelayModifier modifier = new RangeTroposphericDelayModifier(SaastamoinenModel.getStandardModel());
+ final RangeTroposphericDelayModifier modifier = new RangeTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
if (isModifier) {
((Range) measurement).addModifier(modifier);
}
diff --git a/src/test/java/org/orekit/estimation/measurements/RangeAnalyticTest.java b/src/test/java/org/orekit/estimation/measurements/RangeAnalyticTest.java
index a03097bd43..7bb17a32be 100644
--- a/src/test/java/org/orekit/estimation/measurements/RangeAnalyticTest.java
+++ b/src/test/java/org/orekit/estimation/measurements/RangeAnalyticTest.java
@@ -31,7 +31,7 @@
import org.orekit.estimation.Context;
import org.orekit.estimation.EstimationTestUtils;
import org.orekit.estimation.measurements.modifiers.RangeTroposphericDelayModifier;
-import org.orekit.models.earth.troposphere.SaastamoinenModel;
+import org.orekit.models.earth.troposphere.ModifiedSaastamoinenModel;
import org.orekit.orbits.OrbitType;
import org.orekit.orbits.PositionAngleType;
import org.orekit.propagation.Propagator;
@@ -361,7 +361,7 @@ void genericTestStateDerivatives(final boolean isModifier, final boolean isFinit
) {
// Add modifiers if test implies it
- final RangeTroposphericDelayModifier modifier = new RangeTroposphericDelayModifier(SaastamoinenModel.getStandardModel());
+ final RangeTroposphericDelayModifier modifier = new RangeTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
if (isModifier) {
((Range) measurement).addModifier(modifier);
}
@@ -536,7 +536,7 @@ void genericTestParameterDerivatives(final boolean isModifier, final boolean isF
) {
// Add modifiers if test implies it
- final RangeTroposphericDelayModifier modifier = new RangeTroposphericDelayModifier(SaastamoinenModel.getStandardModel());
+ final RangeTroposphericDelayModifier modifier = new RangeTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
if (isModifier) {
((Range) measurement).addModifier(modifier);
}
diff --git a/src/test/java/org/orekit/estimation/measurements/RangeRateTest.java b/src/test/java/org/orekit/estimation/measurements/RangeRateTest.java
index 92b5a2b083..dc5c4d480e 100644
--- a/src/test/java/org/orekit/estimation/measurements/RangeRateTest.java
+++ b/src/test/java/org/orekit/estimation/measurements/RangeRateTest.java
@@ -27,7 +27,7 @@
import org.orekit.estimation.measurements.modifiers.RangeRateTroposphericDelayModifier;
import org.orekit.models.earth.troposphere.EstimatedTroposphericModel;
import org.orekit.models.earth.troposphere.GlobalMappingFunctionModel;
-import org.orekit.models.earth.troposphere.SaastamoinenModel;
+import org.orekit.models.earth.troposphere.ModifiedSaastamoinenModel;
import org.orekit.orbits.OrbitType;
import org.orekit.orbits.PositionAngleType;
import org.orekit.propagation.Propagator;
@@ -463,7 +463,7 @@ public void testStateDerivativesWithModifier() {
double maxRelativeError = 0;
for (final ObservedMeasurement> measurement : measurements) {
- final RangeRateTroposphericDelayModifier modifier = new RangeRateTroposphericDelayModifier(SaastamoinenModel.getStandardModel(), true);
+ final RangeRateTroposphericDelayModifier modifier = new RangeRateTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel(), true);
((RangeRate) measurement).addModifier(modifier);
//
@@ -610,7 +610,7 @@ public void testParameterDerivativesWithModifier() {
double maxRelativeError = 0;
for (final ObservedMeasurement> measurement : measurements) {
- final RangeRateTroposphericDelayModifier modifier = new RangeRateTroposphericDelayModifier(SaastamoinenModel.getStandardModel(), true);
+ final RangeRateTroposphericDelayModifier modifier = new RangeRateTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel(), true);
((RangeRate) measurement).addModifier(modifier);
// parameter corresponding to station position offset
diff --git a/src/test/java/org/orekit/estimation/measurements/RangeTest.java b/src/test/java/org/orekit/estimation/measurements/RangeTest.java
index 1c2eef63ed..c4c8fd34fd 100644
--- a/src/test/java/org/orekit/estimation/measurements/RangeTest.java
+++ b/src/test/java/org/orekit/estimation/measurements/RangeTest.java
@@ -33,7 +33,7 @@
import org.orekit.estimation.measurements.modifiers.RangeTroposphericDelayModifier;
import org.orekit.models.earth.troposphere.EstimatedTroposphericModel;
import org.orekit.models.earth.troposphere.NiellMappingFunctionModel;
-import org.orekit.models.earth.troposphere.SaastamoinenModel;
+import org.orekit.models.earth.troposphere.ModifiedSaastamoinenModel;
import org.orekit.orbits.OrbitType;
import org.orekit.orbits.PositionAngleType;
import org.orekit.propagation.Propagator;
@@ -350,7 +350,7 @@ void genericTestStateDerivatives(final boolean isModifier, final boolean printRe
) {
// Add modifiers if test implies it
- final RangeTroposphericDelayModifier modifier = new RangeTroposphericDelayModifier(SaastamoinenModel.getStandardModel());
+ final RangeTroposphericDelayModifier modifier = new RangeTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
if (isModifier) {
((Range) measurement).addModifier(modifier);
}
@@ -501,7 +501,7 @@ void genericTestParameterDerivatives(final boolean isModifier, final boolean pri
) {
// Add modifiers if test implies it
- final RangeTroposphericDelayModifier modifier = new RangeTroposphericDelayModifier(SaastamoinenModel.getStandardModel());
+ final RangeTroposphericDelayModifier modifier = new RangeTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
if (isModifier) {
((Range) measurement).addModifier(modifier);
}
diff --git a/src/test/java/org/orekit/estimation/measurements/TDOATest.java b/src/test/java/org/orekit/estimation/measurements/TDOATest.java
index 12e06855e0..7d22a02194 100644
--- a/src/test/java/org/orekit/estimation/measurements/TDOATest.java
+++ b/src/test/java/org/orekit/estimation/measurements/TDOATest.java
@@ -26,7 +26,7 @@
import org.orekit.estimation.Context;
import org.orekit.estimation.EstimationTestUtils;
import org.orekit.estimation.measurements.modifiers.TDOATroposphericDelayModifier;
-import org.orekit.models.earth.troposphere.SaastamoinenModel;
+import org.orekit.models.earth.troposphere.ModifiedSaastamoinenModel;
import org.orekit.orbits.OrbitType;
import org.orekit.orbits.PositionAngleType;
import org.orekit.propagation.Propagator;
@@ -170,7 +170,7 @@ public void testStateDerivativesWithModifier() {
// create a modifier
final TDOATroposphericDelayModifier modifier =
- new TDOATroposphericDelayModifier(SaastamoinenModel.getStandardModel());
+ new TDOATroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
double maxRelativeError = 0;
for (final ObservedMeasurement> measurement : measurements) {
@@ -328,7 +328,7 @@ public void testParameterDerivativesWithModifier() {
// create a modifier
final TDOATroposphericDelayModifier modifier =
- new TDOATroposphericDelayModifier(SaastamoinenModel.getStandardModel());
+ new TDOATroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
double maxRelativeError = 0;
for (final ObservedMeasurement> measurement : measurements) {
diff --git a/src/test/java/org/orekit/estimation/measurements/TurnAroundRangeAnalyticTest.java b/src/test/java/org/orekit/estimation/measurements/TurnAroundRangeAnalyticTest.java
index 4a2d52d6a0..80a9c81040 100644
--- a/src/test/java/org/orekit/estimation/measurements/TurnAroundRangeAnalyticTest.java
+++ b/src/test/java/org/orekit/estimation/measurements/TurnAroundRangeAnalyticTest.java
@@ -31,7 +31,7 @@
import org.orekit.estimation.Context;
import org.orekit.estimation.EstimationTestUtils;
import org.orekit.estimation.measurements.modifiers.TurnAroundRangeTroposphericDelayModifier;
-import org.orekit.models.earth.troposphere.SaastamoinenModel;
+import org.orekit.models.earth.troposphere.ModifiedSaastamoinenModel;
import org.orekit.orbits.OrbitType;
import org.orekit.orbits.PositionAngleType;
import org.orekit.propagation.Propagator;
@@ -366,7 +366,7 @@ void genericTestStateDerivatives(final boolean isModifier, final boolean isFinit
// Add modifiers if test implies it
final TurnAroundRangeTroposphericDelayModifier modifier =
- new TurnAroundRangeTroposphericDelayModifier(SaastamoinenModel.getStandardModel());
+ new TurnAroundRangeTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
if (isModifier) {
((TurnAroundRange) measurement).addModifier(modifier);
}
@@ -538,7 +538,7 @@ void genericTestParameterDerivatives(final boolean isModifier, final boolean isF
for (final ObservedMeasurement> measurement : measurements) {
// Add modifiers if test implies it
- final TurnAroundRangeTroposphericDelayModifier modifier = new TurnAroundRangeTroposphericDelayModifier(SaastamoinenModel.getStandardModel());
+ final TurnAroundRangeTroposphericDelayModifier modifier = new TurnAroundRangeTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
if (isModifier) {
((TurnAroundRange) measurement).addModifier(modifier);
}
diff --git a/src/test/java/org/orekit/estimation/measurements/TurnAroundRangeTest.java b/src/test/java/org/orekit/estimation/measurements/TurnAroundRangeTest.java
index 2b6187d55d..a9c0eed8fc 100644
--- a/src/test/java/org/orekit/estimation/measurements/TurnAroundRangeTest.java
+++ b/src/test/java/org/orekit/estimation/measurements/TurnAroundRangeTest.java
@@ -31,7 +31,7 @@
import org.orekit.estimation.Context;
import org.orekit.estimation.EstimationTestUtils;
import org.orekit.estimation.measurements.modifiers.TurnAroundRangeTroposphericDelayModifier;
-import org.orekit.models.earth.troposphere.SaastamoinenModel;
+import org.orekit.models.earth.troposphere.ModifiedSaastamoinenModel;
import org.orekit.orbits.OrbitType;
import org.orekit.orbits.PositionAngleType;
import org.orekit.propagation.Propagator;
@@ -309,7 +309,7 @@ void genericTestStateDerivatives(final boolean isModifier, final boolean printRe
// Add modifiers if test implies it
final TurnAroundRangeTroposphericDelayModifier modifier =
- new TurnAroundRangeTroposphericDelayModifier(SaastamoinenModel.getStandardModel());
+ new TurnAroundRangeTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
if (isModifier) {
((TurnAroundRange) measurement).addModifier(modifier);
}
@@ -456,7 +456,7 @@ void genericTestParameterDerivatives(final boolean isModifier, final boolean pri
for (final ObservedMeasurement> measurement : measurements) {
// Add modifiers if test implies it
- final TurnAroundRangeTroposphericDelayModifier modifier = new TurnAroundRangeTroposphericDelayModifier(SaastamoinenModel.getStandardModel());
+ final TurnAroundRangeTroposphericDelayModifier modifier = new TurnAroundRangeTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
if (isModifier) {
((TurnAroundRange) measurement).addModifier(modifier);
}
diff --git a/src/test/java/org/orekit/estimation/measurements/modifiers/TropoModifierTest.java b/src/test/java/org/orekit/estimation/measurements/modifiers/TropoModifierTest.java
index 67a5b4b684..a9610c86cf 100644
--- a/src/test/java/org/orekit/estimation/measurements/modifiers/TropoModifierTest.java
+++ b/src/test/java/org/orekit/estimation/measurements/modifiers/TropoModifierTest.java
@@ -52,7 +52,7 @@
import org.orekit.models.earth.EarthITU453AtmosphereRefraction;
import org.orekit.models.earth.troposphere.EstimatedTroposphericModel;
import org.orekit.models.earth.troposphere.NiellMappingFunctionModel;
-import org.orekit.models.earth.troposphere.SaastamoinenModel;
+import org.orekit.models.earth.troposphere.ModifiedSaastamoinenModel;
import org.orekit.orbits.OrbitType;
import org.orekit.orbits.PositionAngleType;
import org.orekit.propagation.Propagator;
@@ -97,7 +97,7 @@ public void testRangeTropoModifier() {
1.0, 3.0, 300.0);
propagator.clearStepHandlers();
- final RangeTroposphericDelayModifier modifier = new RangeTroposphericDelayModifier(SaastamoinenModel.getStandardModel());
+ final RangeTroposphericDelayModifier modifier = new RangeTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
for (final ObservedMeasurement> measurement : measurements) {
final AbsoluteDate date = measurement.getDate();
@@ -198,7 +198,7 @@ public void testPhaseTropoModifier() {
1.0, 3.0, 300.0);
propagator.clearStepHandlers();
- final PhaseTroposphericDelayModifier modifier = new PhaseTroposphericDelayModifier(SaastamoinenModel.getStandardModel());
+ final PhaseTroposphericDelayModifier modifier = new PhaseTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
for (final ObservedMeasurement> measurement : measurements) {
final AbsoluteDate date = measurement.getDate();
@@ -306,7 +306,7 @@ public void testTurnAroundRangeTropoModifier() {
1.0, 3.0, 300.0);
propagator.clearStepHandlers();
- final TurnAroundRangeTroposphericDelayModifier modifier = new TurnAroundRangeTroposphericDelayModifier(SaastamoinenModel.getStandardModel());
+ final TurnAroundRangeTroposphericDelayModifier modifier = new TurnAroundRangeTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
for (final ObservedMeasurement> measurement : measurements) {
final AbsoluteDate date = measurement.getDate();
@@ -357,7 +357,7 @@ public void testBistaticRangeTropoModifier() {
propagator.clearStepHandlers();
final BistaticRangeTroposphericDelayModifier modifier =
- new BistaticRangeTroposphericDelayModifier(SaastamoinenModel.getStandardModel());
+ new BistaticRangeTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
for (final ObservedMeasurement> measurement : measurements) {
BistaticRange biRange = (BistaticRange) measurement;
@@ -408,7 +408,7 @@ public void testBistaticRangeRateTropoModifier() {
propagator.clearStepHandlers();
final BistaticRangeRateTroposphericDelayModifier modifier =
- new BistaticRangeRateTroposphericDelayModifier(SaastamoinenModel.getStandardModel());
+ new BistaticRangeRateTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
for (final ObservedMeasurement> measurement : measurements) {
BistaticRangeRate biRangeRate = (BistaticRangeRate) measurement;
@@ -520,7 +520,7 @@ public void testTDOATropoModifier() {
propagator.clearStepHandlers();
final TDOATroposphericDelayModifier modifier =
- new TDOATroposphericDelayModifier(SaastamoinenModel.getStandardModel());
+ new TDOATroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
for (final ObservedMeasurement> measurement : measurements) {
TDOA tdoa = (TDOA) measurement;
@@ -625,7 +625,7 @@ public void testRangeRateTropoModifier() {
1.0, 3.0, 300.0);
propagator.clearStepHandlers();
- final RangeRateTroposphericDelayModifier modifier = new RangeRateTroposphericDelayModifier(SaastamoinenModel.getStandardModel(), false);
+ final RangeRateTroposphericDelayModifier modifier = new RangeRateTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel(), false);
for (final ObservedMeasurement> measurement : measurements) {
final AbsoluteDate date = measurement.getDate();
@@ -722,7 +722,7 @@ public void testAngularTropoModifier() {
1.0, 3.0, 300.0);
propagator.clearStepHandlers();
- final AngularTroposphericDelayModifier modifier = new AngularTroposphericDelayModifier(SaastamoinenModel.getStandardModel());
+ final AngularTroposphericDelayModifier modifier = new AngularTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
for (final ObservedMeasurement> measurement : measurements) {
final AbsoluteDate date = measurement.getDate();
diff --git a/src/test/java/org/orekit/models/earth/troposphere/ModifiedSaastamoinenModelTest.java b/src/test/java/org/orekit/models/earth/troposphere/ModifiedSaastamoinenModelTest.java
new file mode 100644
index 0000000000..02c2008ba0
--- /dev/null
+++ b/src/test/java/org/orekit/models/earth/troposphere/ModifiedSaastamoinenModelTest.java
@@ -0,0 +1,437 @@
+/* Copyright 2011-2012 Space Applications Services
+ * Licensed to CS Communication & Systèmes (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.models.earth.troposphere;
+
+import org.hipparchus.CalculusFieldElement;
+import org.hipparchus.Field;
+import org.hipparchus.util.Binary64Field;
+import org.hipparchus.util.FastMath;
+import org.hipparchus.util.Precision;
+import org.junit.jupiter.api.Assertions;
+import org.junit.jupiter.api.BeforeEach;
+import org.junit.jupiter.api.Test;
+import org.orekit.Utils;
+import org.orekit.bodies.FieldGeodeticPoint;
+import org.orekit.bodies.GeodeticPoint;
+import org.orekit.errors.OrekitException;
+import org.orekit.errors.OrekitMessages;
+import org.orekit.time.AbsoluteDate;
+import org.orekit.time.FieldAbsoluteDate;
+
+
+public class ModifiedSaastamoinenModelTest {
+
+ private static double epsilon = 1e-6;
+
+ private double[][] expectedValues;
+
+ private double[] elevations;
+
+ private double[] heights;
+
+ @Test
+ public void testIssue1078() {
+ try {
+ new ModifiedSaastamoinenModel(273.16 + 18, 1013.25, 50.0);
+ } catch (OrekitException oe) {
+ Assertions.assertEquals(OrekitMessages.INVALID_PARAMETER_RANGE, oe.getSpecifier());
+ }
+ try {
+ new ModifiedSaastamoinenModel(273.16 + 18, 1013.25, -50.0);
+ } catch (OrekitException oe) {
+ Assertions.assertEquals(OrekitMessages.INVALID_PARAMETER_RANGE, oe.getSpecifier());
+ }
+ }
+
+ @Test
+ public void testFixedElevation() {
+ Utils.setDataRoot("atmosphere");
+ ModifiedSaastamoinenModel model = ModifiedSaastamoinenModel.getStandardModel();
+ double lastDelay = Double.MAX_VALUE;
+ // delay shall decline with increasing height of the station
+ for (double height = 0; height < 5000; height += 100) {
+ final double delay = model.pathDelay(FastMath.toRadians(5), new GeodeticPoint(0.0, 0.0, height), null, AbsoluteDate.J2000_EPOCH);
+ Assertions.assertTrue(Precision.compareTo(delay, lastDelay, epsilon) < 0);
+ lastDelay = delay;
+ }
+ }
+
+ @Test
+ public void testFieldFixedElevation() {
+ doTestFieldFixedElevation(Binary64Field.getInstance());
+ }
+
+ private > void doTestFieldFixedElevation(final Field field) {
+ final T zero = field.getZero();
+ Utils.setDataRoot("atmosphere");
+ ModifiedSaastamoinenModel model = ModifiedSaastamoinenModel.getStandardModel();
+ T lastDelay = zero.add(Double.MAX_VALUE);
+ // delay shall decline with increasing height of the station
+ for (double height = 0; height < 5000; height += 100) {
+ final T delay = model.pathDelay(zero.add(FastMath.toRadians(5)), new FieldGeodeticPoint<>(zero, zero, zero.add(height)), null, FieldAbsoluteDate.getJ2000Epoch(field));
+ Assertions.assertTrue(Precision.compareTo(delay.getReal(), lastDelay.getReal(), epsilon) < 0);
+ lastDelay = delay;
+ }
+ }
+
+ @Test
+ public void testFixedHeight() {
+ Utils.setDataRoot("atmosphere");
+ ModifiedSaastamoinenModel model = ModifiedSaastamoinenModel.getStandardModel();
+ double lastDelay = Double.MAX_VALUE;
+ // delay shall decline with increasing elevation angle
+ for (double elev = 10d; elev < 90d; elev += 8d) {
+ final double delay = model.pathDelay(FastMath.toRadians(elev), new GeodeticPoint(0.0, 0.0, 350.0), null, AbsoluteDate.J2000_EPOCH);
+ Assertions.assertTrue(Precision.compareTo(delay, lastDelay, epsilon) < 0);
+ lastDelay = delay;
+ }
+ }
+
+ @Test
+ public void testFieldFixedHeight() {
+ doTestFieldFixedHeight(Binary64Field.getInstance());
+ }
+
+ private > void doTestFieldFixedHeight(final Field field) {
+ final T zero = field.getZero();
+ Utils.setDataRoot("atmosphere");
+ ModifiedSaastamoinenModel model = ModifiedSaastamoinenModel.getStandardModel();
+ T lastDelay = zero.add(Double.MAX_VALUE);
+ // delay shall decline with increasing elevation angle
+ for (double elev = 10d; elev < 90d; elev += 8d) {
+ final T delay = model.pathDelay(zero.add(FastMath.toRadians(elev)), new FieldGeodeticPoint<>(zero, zero, zero.add(350.0)), null, FieldAbsoluteDate.getJ2000Epoch(field));
+ Assertions.assertTrue(Precision.compareTo(delay.getReal(), lastDelay.getReal(), epsilon) < 0);
+ lastDelay = delay;
+ }
+ }
+
+ @Test
+ public void NoFile() {
+ Utils.setDataRoot("atmosphere");
+ try {
+ new ModifiedSaastamoinenModel(273.16 + 18, 1013.25, 0.5, "^non-existent-file$");
+ Assertions.fail("an exception should have been thrown");
+ } catch (OrekitException oe) {
+ Assertions.assertEquals(OrekitMessages.UNABLE_TO_FIND_FILE, oe.getSpecifier());
+ Assertions.assertEquals("non-existent-file", oe.getParts()[0]);
+ }
+ }
+
+ @Test
+ public void compareDefaultAndLoaded() {
+ Utils.setDataRoot("atmosphere");
+ ModifiedSaastamoinenModel defaultModel = new ModifiedSaastamoinenModel(273.16 + 18, 1013.25, 0.5, null);
+ ModifiedSaastamoinenModel loadedModel = new ModifiedSaastamoinenModel(273.16 + 18, 1013.25, 0.5, ModifiedSaastamoinenModel.DELTA_R_FILE_NAME);
+ double[] heights = new double[] {
+ 0.0, 250.0, 500.0, 750.0, 1000.0, 1250.0, 1500.0, 1750.0, 2000.0, 2250.0, 2500.0, 2750.0, 3000.0, 3250.0,
+ 3500.0, 3750.0, 4000.0, 4250.0, 4500.0, 4750.0, 5000.0
+ };
+ double[] elevations = new double[] {
+ FastMath.toRadians(10.0), FastMath.toRadians(15.0), FastMath.toRadians(20.0),
+ FastMath.toRadians(25.0), FastMath.toRadians(30.0), FastMath.toRadians(35.0),
+ FastMath.toRadians(40.0), FastMath.toRadians(45.0), FastMath.toRadians(50.0),
+ FastMath.toRadians(55.0), FastMath.toRadians(60.0), FastMath.toRadians(65.0),
+ FastMath.toRadians(70.0), FastMath.toRadians(75.0), FastMath.toRadians(80.0),
+ FastMath.toRadians(85.0), FastMath.toRadians(90.0)
+ };
+ for (int h = 0; h < heights.length; h++) {
+ for (int e = 0; e < elevations.length; e++) {
+ double height = heights[h];
+ double elevation = elevations[e];
+ double expectedValue = defaultModel.pathDelay(elevation, new GeodeticPoint(0.0, 0.0, height), null, AbsoluteDate.J2000_EPOCH);
+ double actualValue = loadedModel.pathDelay(elevation, new GeodeticPoint(0.0, 0.0, height), null, AbsoluteDate.J2000_EPOCH);
+ Assertions.assertEquals(expectedValue, actualValue, epsilon, "For height=" + height + " elevation = " +
+ FastMath.toDegrees(elevation) + " precision not met");
+ }
+ }
+ }
+
+ @Test
+ public void testNegativeHeight() {
+ Utils.setDataRoot("atmosphere");
+ ModifiedSaastamoinenModel model = ModifiedSaastamoinenModel.getStandardModel();
+ final double height = -500.0;
+ for (double elevation = 0; elevation < FastMath.PI; elevation += 0.1) {
+ Assertions.assertEquals(model.pathDelay(elevation, new GeodeticPoint(0.0, 0.0, 0.0), null, AbsoluteDate.J2000_EPOCH),
+ model.pathDelay(elevation, new GeodeticPoint(0.0, 0.0, height), null, AbsoluteDate.J2000_EPOCH),
+ 1.e-10);
+ }
+ }
+
+ @Test
+ public void testFieldNegativeHeight() {
+ doTestFieldNegativeHeight(Binary64Field.getInstance());
+ }
+
+ private > void doTestFieldNegativeHeight(final Field field) {
+ final T zero = field.getZero();
+ Utils.setDataRoot("atmosphere");
+ ModifiedSaastamoinenModel model = ModifiedSaastamoinenModel.getStandardModel();
+ final T height = zero.subtract(500.0);
+ for (double elevation = 0; elevation < FastMath.PI; elevation += 0.1) {
+ Assertions.assertEquals(model.pathDelay(zero.add(elevation), new FieldGeodeticPoint<>(zero, zero, zero), null, FieldAbsoluteDate.getJ2000Epoch(field)).getReal(),
+ model.pathDelay(zero.add(elevation), new FieldGeodeticPoint<>(zero, zero, zero.add(height)), null, FieldAbsoluteDate.getJ2000Epoch(field)).getReal(),
+ 1.e-10);
+ }
+ }
+
+ @Test
+ public void testIssue654LowElevation() {
+ Utils.setDataRoot("atmosphere");
+ ModifiedSaastamoinenModel model = ModifiedSaastamoinenModel.getStandardModel();
+
+ // Test model new setter/getter
+ model.setLowElevationThreshold(1e-3);
+ Assertions.assertEquals(1.e-3, model.getLowElevationThreshold(), 0.);
+
+ // Reset to default value
+ model.setLowElevationThreshold(ModifiedSaastamoinenModel.DEFAULT_LOW_ELEVATION_THRESHOLD);
+ double lowElevationPathDelay = model.pathDelay(0.001, new GeodeticPoint(0.0, 0.0, 0.0), null, AbsoluteDate.J2000_EPOCH);
+ Assertions.assertTrue(lowElevationPathDelay > 0.);
+ Assertions.assertEquals(model.pathDelay(model.getLowElevationThreshold(), new GeodeticPoint(0.0, 0.0, 0.0), null, AbsoluteDate.J2000_EPOCH),
+ lowElevationPathDelay, 1.e-10);
+ }
+
+ @Test
+ public void testIssue654FieldLowElevation() { doTestFieldLowElevation(Binary64Field.getInstance()); }
+
+ private > void doTestFieldLowElevation(final Field field) {
+ final T zero = field.getZero();
+ Utils.setDataRoot("atmosphere");
+ ModifiedSaastamoinenModel model = ModifiedSaastamoinenModel.getStandardModel();
+ final T elevation = zero.add(0.001);
+ double lowElevationPathDelay = model.pathDelay(zero.add(elevation), new FieldGeodeticPoint<>(zero, zero, zero), null,
+ FieldAbsoluteDate.getJ2000Epoch(field)).getReal();
+ double thresholdElevationPathDelay = model.pathDelay(zero.add(model.getLowElevationThreshold()), new FieldGeodeticPoint<>(zero, zero, zero),
+ null, FieldAbsoluteDate.getJ2000Epoch(field)).getReal();
+ Assertions.assertTrue(lowElevationPathDelay > 0.);
+ Assertions.assertEquals(thresholdElevationPathDelay, lowElevationPathDelay, 1.e-10);
+ }
+
+ @Test
+ public void compareExpectedValues() {
+ Utils.setDataRoot("atmosphere");
+ ModifiedSaastamoinenModel model = ModifiedSaastamoinenModel.getStandardModel();
+
+ for (int h = 0; h < heights.length; h++) {
+ for (int e = 0; e < elevations.length; e++) {
+ double height = heights[h];
+ double elevation = elevations[e];
+ double expectedValue = expectedValues[h][e];
+ double actualValue = model.pathDelay(elevation, new GeodeticPoint(0.0, 0.0, height), null, AbsoluteDate.J2000_EPOCH);
+ Assertions.assertEquals(expectedValue, actualValue, epsilon, "For height=" + height + " elevation = " +
+ FastMath.toDegrees(elevation) + " precision not met");
+ }
+ }
+ }
+
+ @Test
+ public void compareFieldExpectedValues() {
+ doCompareFieldExpectedValues(Binary64Field.getInstance());
+ }
+
+ private > void doCompareFieldExpectedValues(final Field field) {
+ final T zero = field.getZero();
+ Utils.setDataRoot("atmosphere");
+ ModifiedSaastamoinenModel model = ModifiedSaastamoinenModel.getStandardModel();
+
+ for (int h = 0; h < heights.length; h++) {
+ for (int e = 0; e < elevations.length; e++) {
+ T height = zero.add(heights[h]);
+ T elevation = zero.add(elevations[e]);
+ double expectedValue = expectedValues[h][e];
+ T actualValue = model.pathDelay(elevation, new FieldGeodeticPoint<>(zero, zero, zero.add(height)), null, FieldAbsoluteDate.getJ2000Epoch(field));
+ Assertions.assertEquals(expectedValue, actualValue.getReal(), epsilon, "For height=" + height + " elevation = " +
+ FastMath.toDegrees(elevation.getReal()) + " precision not met");
+ }
+ }
+ }
+
+ @Test
+ public void testIssue572() {
+ Utils.setDataRoot("atmosphere");
+ ModifiedSaastamoinenModel model = ModifiedSaastamoinenModel.getStandardModel();
+ final double height = 6000.0;
+ for (double elevation = 0; elevation < FastMath.PI; elevation += 0.1) {
+ Assertions.assertEquals(model.pathDelay(elevation, new GeodeticPoint(0.0, 0.0, 5000.0), null, AbsoluteDate.J2000_EPOCH), model.pathDelay(elevation, new GeodeticPoint(0.0, 0.0, height), null, AbsoluteDate.J2000_EPOCH), 1.e-10);
+ }
+ }
+
+ @Test
+ public void testFieldIssue572() {
+ doTestFieldIssue572(Binary64Field.getInstance());
+ }
+
+ private > void doTestFieldIssue572(final Field field) {
+ final T zero = field.getZero();
+ Utils.setDataRoot("atmosphere");
+ ModifiedSaastamoinenModel model = ModifiedSaastamoinenModel.getStandardModel();
+ final T height = zero.add(6000.0);
+ for (double elevation = 0; elevation < FastMath.PI; elevation += 0.1) {
+ Assertions.assertEquals(model.pathDelay(zero.add(elevation),new FieldGeodeticPoint<>(zero, zero, zero.add(5000.0)), null, FieldAbsoluteDate.getJ2000Epoch(field)).getReal(),
+ model.pathDelay(zero.add(elevation), new FieldGeodeticPoint<>(zero, zero, zero.add(height)), null, FieldAbsoluteDate.getJ2000Epoch(field)).getReal(),
+ 1.e-10);
+ }
+ }
+
+ @BeforeEach
+ public void setUp() throws Exception {
+ heights = new double[] {
+ 0.0, 250.0, 500.0, 750.0, 1000.0, 1250.0, 1500.0, 1750.0, 2000.0, 2250.0, 2500.0, 2750.0, 3000.0, 3250.0,
+ 3500.0, 3750.0, 4000.0, 4250.0, 4500.0, 4750.0, 5000.0
+ };
+
+ elevations = new double[] {
+ FastMath.toRadians(10.0), FastMath.toRadians(15.0), FastMath.toRadians(20.0),
+ FastMath.toRadians(25.0), FastMath.toRadians(30.0), FastMath.toRadians(35.0),
+ FastMath.toRadians(40.0), FastMath.toRadians(45.0), FastMath.toRadians(50.0),
+ FastMath.toRadians(55.0), FastMath.toRadians(60.0), FastMath.toRadians(65.0),
+ FastMath.toRadians(70.0), FastMath.toRadians(75.0), FastMath.toRadians(80.0),
+ FastMath.toRadians(85.0), FastMath.toRadians(90.0)
+ };
+
+ expectedValues = new double[][] {
+ {
+ 13.517414068807756, 9.204443522241771, 7.0029750138616835, 5.681588299211439, 4.8090544808193805,
+ 4.196707503563898, 3.7474156937027994, 3.408088733958258, 3.1468182787091985, 2.943369134588668,
+ 2.784381959210485, 2.6607786449639343, 2.5662805210588195, 2.496526411065139, 2.4485331622035984,
+ 2.420362989334734, 2.4109238764096896
+ },
+ {
+ 13.004543691989646, 8.85635958366884, 6.7385672069739835, 5.467398852004842, 4.627733401816586,
+ 4.038498891518074, 3.606157486803878, 3.279627416773643, 3.0282066103153564, 2.8324244661904134,
+ 2.6794262487842104, 2.56047665894495, 2.4695340768552505, 2.402401959167246, 2.356209754788866,
+ 2.329092816778967, 2.3200003434082928
+ },
+ {
+ 12.531363728988735, 8.534904937493899, 6.494310751299656, 5.269517663702665, 4.460196658874199,
+ 3.892306407739391, 3.475621589928269, 3.1609130970914956, 2.9185920283074447, 2.729893581384905,
+ 2.582428928493012, 2.4677793389442715, 2.380122124673593, 2.3154128046624067, 2.2708848382055904,
+ 2.2447411392156, 2.2359689784370986
+ },
+ {
+ 12.09063673875363, 8.235209195291242, 6.266604991324561, 5.084562550097057, 4.303522687504001,
+ 3.755568409663704, 3.353518885593948, 3.0498713508982442, 2.8160742841783275, 2.6340206738065692,
+ 2.4917559301110956, 2.3811563829818176, 2.2966034070866277, 2.234194258980332, 2.191259347593356,
+ 2.1660645619383274, 2.1576326612520003
+ },
+ {
+ 11.67727244511714, 7.953871771343415, 6.052791636499061, 4.910850401669602, 4.156351680934609,
+ 3.6271143683534492, 3.2388081756428404, 2.9455492155570195, 2.7197591598098305, 2.5439482552015917,
+ 2.4065694710150236, 2.2997761077823076, 2.218141111456481, 2.157894809849032, 2.1164586386563973,
+ 2.0921576169523175, 2.0840478264673044
+ },
+ {
+ 11.286432636721148, 7.688212194124222, 5.850570088713712, 4.747059102172618, 4.017661723553029,
+ 3.506085459183713, 3.1307362765144857, 2.8472616350388877, 2.6290036896596245, 2.459056474904878,
+ 2.3262584011701235, 2.223024699820715, 2.1441094810550236, 2.085868912585518, 2.0458105091517886,
+ 2.022315205377469, 2.0144705937765144
+ },
+ {
+ 10.915292255872545, 7.435769456080562, 5.6583502024136285, 4.591362033342799, 3.8858133055608204,
+ 3.391020479976734, 3.027986150306355, 2.7538117534167346, 2.5427137172039873, 2.3783406833254412,
+ 2.249897295933348, 2.1500476936060946, 2.0737181577274097, 2.0173844567055803, 1.978635920228296,
+ 1.9559066302818124, 1.9483141307804097
+ },
+ {
+ 10.560838722447652, 7.194507733765155, 5.474676340193435, 4.442196283017724, 3.759852141271757,
+ 3.281095182764508, 2.9298266864995415, 2.6645376694677676, 2.4602800254909183, 2.3012323491024245,
+ 2.1769492208902093, 2.080332602007238, 2.0064732734566384, 1.9519612730357911, 1.9144640973301363,
+ 1.8924666054100323, 1.8851149566358782
+ },
+ {
+ 10.221303680396087, 6.9632552008410915, 5.298576794548074, 4.299160340784349, 3.6390721146637293,
+ 3.175686404496119, 2.8356974323630437, 2.5789272031073223, 2.381228081225778, 2.2272865154903485,
+ 2.106992477081925, 2.0134758868645615, 1.9419852149640153, 1.8892200435803028, 1.8529228069725603,
+ 1.8316270449308856, 1.8245063513318645
+ },
+ {
+ 9.894629211990411, 6.740704058398638, 5.129125614634508, 4.161737017461952, 3.5228709097629975,
+ 3.0742748111390386, 2.7451382632192436, 2.4965641131187946, 2.305174987172354, 2.156146000844558,
+ 2.039689834231423, 1.949155743414788, 1.8799439192071106, 1.8288593407337934, 1.7937165420599064,
+ 1.7730958998798743, 1.7661974033814984
+ },
+ {
+ 9.579864280378851, 6.526143322382175, 4.965721065018929, 4.029207163135991, 3.4108058435845767,
+ 2.9764687866827866, 2.6577964388561925, 2.417125714868741, 2.2318215659378273, 2.087530132722662,
+ 1.9747751921856533, 1.8871174620329965, 1.820103416896286, 1.770639660836324, 1.7366102498117952,
+ 1.7166407238790062, 1.709956524792288
+ },
+ {
+ 9.274887896199909, 6.318551036055798, 4.807695974007752, 3.901070574943598, 3.302472329989199,
+ 2.881925175414427, 2.5733712370891264, 2.3403420235759187, 2.1609208036663294, 2.021209397507298,
+ 1.9120324913823707, 1.8271553141955852, 1.7622658032319802, 1.7143688314998151, 1.681415677692901,
+ 1.662075550126913, 1.6555984999945992
+ },
+ {
+ 8.979896361522131, 6.117657835260275, 4.654740323946152, 3.777036627543426, 3.197606518234132,
+ 2.7904044409768964, 2.4916434285107703, 2.266010367769514, 2.0922834230246177, 1.9570053034359607,
+ 1.851291869170061, 1.7691062591782465, 1.706273315177691, 1.6598930167213255, 1.627981703939375,
+ 1.6092508503238145, 1.6029743261170657
+ },
+ {
+ 8.69399537567174, 5.922971478822413, 4.5066379960473855, 3.6569305214071046, 3.0956861718504136,
+ 2.701448680145781, 2.412205508374445, 2.193765237935503, 2.025580422503155, 1.8946215440529706,
+ 1.7922869252688411, 1.7127316699243513, 1.6519133992903239, 1.6070243276625142, 1.5761438889147779,
+ 1.5580245420477885, 1.5519632546752067
+ },
+ {
+ 8.416815377025097, 5.734136251055505, 4.362963429392922, 3.5404077519897172, 2.996794592537471,
+ 2.615133166436779, 2.3351235507871273, 2.1236617698358717, 1.9608543092876316, 1.8340865056076205,
+ 1.735030640851246, 1.658028018024244, 1.5991650567003197, 1.555723443805895, 1.5258438192326151,
+ 1.5083184020062343, 1.5024665667687351
+ },
+ {
+ 8.1478867312736, 5.550837062543208, 4.223478184395355, 3.4272753528502906, 2.9007686780115858,
+ 2.531315719772937, 2.2602706846943197, 2.055584632739777, 1.8979986259230126, 1.7753006325382452,
+ 1.679428848769866, 1.60490532174873, 1.5479415024951926, 1.5059059371968162, 1.4769986856932136,
+ 1.4600505675451787, 1.4544027112557927
+ },
+ {
+ 7.886816833128103, 5.372810504566989, 4.087982930125399, 3.3173720077392095, 2.807472077886753,
+ 2.449877480332473, 2.187540848323867, 1.9894374123646559, 1.8369243760154002, 1.718180698301642,
+ 1.6254028270926364, 1.5532883580952295, 1.4981701861499142, 1.457501227681867, 1.4295392613913276,
+ 1.4131526030534576, 1.4077035129433761
+ },
+ {
+ 7.632757849842542, 5.199465832889435, 3.956158045029617, 3.2103200618121055, 2.7169983015566497,
+ 2.370922636048298, 2.1170374348830436, 1.9253162741272756, 1.7777165574325338, 1.6627979450991903,
+ 1.5730082599508786, 1.5032159322097494, 1.4498720192717967, 1.4105115179325056, 1.3834483541077152,
+ 1.3675873164708097, 1.3623112195283247
+ },
+ {
+ 7.385939247167236, 5.03097891384785, 3.8280091975097514, 3.1062430912053003, 2.629039083023718,
+ 2.294159790631063, 2.048490000181759, 1.8629731967599608, 1.720149999888605, 1.608950046027118,
+ 1.5220654734302106, 1.454530760098946, 1.402911820651357, 1.364823439078982, 1.3386341212762891,
+ 1.3232841113878862, 1.3181762050118586
+ },
+ {
+ 7.146111766814672, 4.867182513816216, 3.7034098358166165, 3.005038679030282, 2.5435078557931674,
+ 2.219513482041795, 1.981831207440737, 1.8023469685072222, 1.664168201116958, 1.5565841619968868,
+ 1.472524488341563, 1.407185083983383, 1.3572435292187972, 1.320392181006258, 1.295052611935796,
+ 1.2801995392223198, 1.2752551861465835
+ },
+ {
+ 6.913054711276373, 4.707928561310275, 3.58224790888857, 2.906616143639732, 2.460327972424674,
+ 2.14691689491338, 1.9170014355353862, 1.7433833914442447, 1.6097211330539392, 1.5056535083847744,
+ 1.4243410522646305, 1.3611366183164608, 1.3128263617229614, 1.277178068079704, 1.2526649111609156,
+ 1.238295129863562, 1.2335098392123367
+ }
+ };
+ }
+
+}
diff --git a/src/test/java/org/orekit/models/earth/troposphere/SaastamoinenModelTest.java b/src/test/java/org/orekit/models/earth/troposphere/SaastamoinenModelTest.java
index f3676b0b5f..e02b8b1c77 100644
--- a/src/test/java/org/orekit/models/earth/troposphere/SaastamoinenModelTest.java
+++ b/src/test/java/org/orekit/models/earth/troposphere/SaastamoinenModelTest.java
@@ -32,7 +32,7 @@
import org.orekit.time.AbsoluteDate;
import org.orekit.time.FieldAbsoluteDate;
-
+@Deprecated
public class SaastamoinenModelTest {
private static double epsilon = 1e-6;
From 18728a1003b6fd3470a2a354cee1ee06172a56c7 Mon Sep 17 00:00:00 2001
From: Luc Maisonobe
Date: Tue, 5 Dec 2023 10:27:13 +0100
Subject: [PATCH 013/359] Work In Progress on tropospheric models.
---
.../CanonicalSaastamoinenModel.java | 260 ++++++++++++++++++
.../earth/troposphere/GiacomoDavis.java | 92 +++++++
.../earth/troposphere/MendesPavlisModel.java | 25 +-
.../ModifiedSaastamoinenModel.java | 2 +-
.../WaterVaporPressureProvider.java | 45 +++
.../FieldPressureTemperatureHumidity.java | 71 +++++
.../weather/PressureTemperatureHumidity.java | 68 +++++
7 files changed, 538 insertions(+), 25 deletions(-)
create mode 100644 src/main/java/org/orekit/models/earth/troposphere/CanonicalSaastamoinenModel.java
create mode 100644 src/main/java/org/orekit/models/earth/troposphere/GiacomoDavis.java
create mode 100644 src/main/java/org/orekit/models/earth/troposphere/WaterVaporPressureProvider.java
create mode 100644 src/main/java/org/orekit/models/earth/weather/FieldPressureTemperatureHumidity.java
create mode 100644 src/main/java/org/orekit/models/earth/weather/PressureTemperatureHumidity.java
diff --git a/src/main/java/org/orekit/models/earth/troposphere/CanonicalSaastamoinenModel.java b/src/main/java/org/orekit/models/earth/troposphere/CanonicalSaastamoinenModel.java
new file mode 100644
index 0000000000..4d3e204344
--- /dev/null
+++ b/src/main/java/org/orekit/models/earth/troposphere/CanonicalSaastamoinenModel.java
@@ -0,0 +1,260 @@
+/* Copyright 2023 Thales Alenia Space
+ * Licensed to CS Communication & Systèmes (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.models.earth.troposphere;
+
+import java.util.Collections;
+import java.util.List;
+import java.util.regex.Pattern;
+
+import org.hipparchus.CalculusFieldElement;
+import org.hipparchus.Field;
+import org.hipparchus.analysis.interpolation.BilinearInterpolatingFunction;
+import org.hipparchus.analysis.interpolation.LinearInterpolator;
+import org.hipparchus.analysis.polynomials.PolynomialFunction;
+import org.hipparchus.analysis.polynomials.PolynomialSplineFunction;
+import org.hipparchus.util.FastMath;
+import org.orekit.annotation.DefaultDataContext;
+import org.orekit.bodies.FieldGeodeticPoint;
+import org.orekit.bodies.GeodeticPoint;
+import org.orekit.data.DataContext;
+import org.orekit.data.DataProvidersManager;
+import org.orekit.errors.OrekitException;
+import org.orekit.errors.OrekitMessages;
+import org.orekit.time.AbsoluteDate;
+import org.orekit.time.FieldAbsoluteDate;
+import org.orekit.utils.InterpolationTableLoader;
+import org.orekit.utils.ParameterDriver;
+
+/** The canonical Saastamoinen model.
+ *
+ * Estimates the path delay imposed to
+ * electro-magnetic signals by the troposphere according to the formula:
+ * \[
+ * \delta = \frac{0.002277}{\cos z (1 - 0.00266\cos 2\varphi - 0.00028 h})}
+ * \left[P+(\frac{1255}{T}+0.005)e - B(h) \tan^2 z\right]
+ * \]
+ * with the following input data provided to the model:
+ *
+ * - z: zenith angle
+ * - P: atmospheric pressure
+ * - T: temperature
+ * - e: partial pressure of water vapor
+ *
+ * @author Luc Maisonobe
+ * @see "J Saastamoinen, Atmospheric Correction for the Troposphere and Stratosphere in Radio
+ * Ranging of Satellites"
+ * @since 12.1
+ */
+public class CanonicalSaastamoinenModel implements DiscreteTroposphericModel {
+
+ /** Default lowest acceptable elevation angle [rad]. */
+ public static final double DEFAULT_LOW_ELEVATION_THRESHOLD = 0.05;
+
+ /** X values for the B function (table 1 in reference paper). */
+ private static final double[] X_VALUES_FOR_B = {
+ 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0
+ };
+
+ /** E values for the B function (table 1 in reference paper). */
+ private static final double[] Y_VALUES_FOR_B = {
+ 1.16, 1.13, 1.10, 1.07, 1.04, 1.01, 0.94, 0.88, 0.82, 0.76, 0.66, 0.57, 0.49
+ };
+
+ /** Interpolation function for the B correction term. */
+ private final PolynomialSplineFunction bFunction;
+
+ /** The temperature at the station [K]. */
+ private double t0;
+
+ /** The atmospheric pressure [mbar]. */
+ private double p0;
+
+ /** The humidity [percent]. */
+ private double r0;
+
+ /** Lowest acceptable elevation angle [rad]. */
+ private double lowElevationThreshold;
+
+ /**
+ * Create a new Saastamoinen model for the troposphere using the given environmental
+ * conditions and table from the reference book.
+ *
+ * @param t0 the temperature at the station [K]
+ * @param p0 the atmospheric pressure at the station [mbar]
+ * @param r0 the humidity at the station [fraction] (50% -> 0.5)
+ */
+ public CanonicalSaastamoinenModel(final double t0, final double p0, final double r0) {
+ checkParameterRangeInclusive("humidity", r0, 0.0, 1.0);
+ this.t0 = t0;
+ this.p0 = p0;
+ this.r0 = r0;
+ this.bFunction = new LinearInterpolator().interpolate(X_VALUES_FOR_B, Y_VALUES_FOR_B);
+ this.lowElevationThreshold = DEFAULT_LOW_ELEVATION_THRESHOLD;
+ }
+
+ /** Create a new Saastamoinen model using a standard atmosphere model.
+ *
+ *
+ * - temperature: 18 degree Celsius
+ *
- pressure: 1013.25 mbar
+ *
- humidity: 50%
+ *
+ *
+ * @return a Saastamoinen model with standard environmental values
+ */
+ public static CanonicalSaastamoinenModel getStandardModel() {
+ return new CanonicalSaastamoinenModel(273.16 + 18, 1013.25, 0.5);
+ }
+
+ /** Check if the given parameter is within an acceptable range.
+ * The bounds are inclusive: an exception is raised when either of those conditions are met:
+ *
+ * - The parameter is strictly greater than upperBound
+ * - The parameter is strictly lower than lowerBound
+ *
+ *
+ * In either of these cases, an OrekitException is raised.
+ *
+ * @param parameterName name of the parameter
+ * @param parameter value of the parameter
+ * @param lowerBound lower bound of the acceptable range (inclusive)
+ * @param upperBound upper bound of the acceptable range (inclusive)
+ */
+ private void checkParameterRangeInclusive(final String parameterName, final double parameter,
+ final double lowerBound, final double upperBound) {
+ if (parameter < lowerBound || parameter > upperBound) {
+ throw new OrekitException(OrekitMessages.INVALID_PARAMETER_RANGE, parameterName,
+ parameter, lowerBound, upperBound);
+ }
+ }
+
+ /** {@inheritDoc}
+ *
+ * The Saastamoinen model is not defined for altitudes below 0.0. for continuity
+ * reasons, we use the value for h = 0 when altitude is negative.
+ *
+ *
+ * There are also numerical issues for elevation angles close to zero. For continuity reasons,
+ * elevations lower than a threshold will use the value obtained
+ * for the threshold itself.
+ *
+ * @see #getLowElevationThreshold()
+ * @see #setLowElevationThreshold(double)
+ */
+ @Override
+ public double pathDelay(final double elevation, final GeodeticPoint point,
+ final double[] parameters, final AbsoluteDate date) {
+
+ // there are no data in the model for negative altitudes and altitude bigger than 5000 m
+ // limit the height to a range of [0, 5000] m
+ final double fixedHeight = FastMath.min(FastMath.max(0, point.getAltitude()), 5000);
+
+ // the corrected temperature using a temperature gradient of -6.5 K/km
+ final double T = t0 - 6.5e-3 * fixedHeight;
+ // the corrected pressure
+ final double P = p0 * FastMath.pow(1.0 - 2.26e-5 * fixedHeight, 5.225);
+ // the corrected humidity
+ final double R = r0 * FastMath.exp(-6.396e-4 * fixedHeight);
+
+ // interpolate the b correction term
+ final double B = bFunction.value(fixedHeight / 1e3);
+ // calculate e
+ final double e = R * FastMath.exp(eFunction.value(T));
+
+ // calculate the zenith angle from the elevation
+ final double z = FastMath.abs(0.5 * FastMath.PI - FastMath.max(elevation, lowElevationThreshold));
+
+ // calculate the path delay in m
+ final double tan = FastMath.tan(z);
+ final double delta = 2.277e-3 / FastMath.cos(z) *
+ (P + (1255d / T + 5e-2) * e - B * tan * tan);
+
+ return delta;
+ }
+
+ /** {@inheritDoc}
+ *
+ * The Saastamoinen model is not defined for altitudes below 0.0. for continuity
+ * reasons, we use the value for h = 0 when altitude is negative.
+ *
+ *
+ * There are also numerical issues for elevation angles close to zero. For continuity reasons,
+ * elevations lower than a threshold will use the value obtained
+ * for the threshold itself.
+ *
+ * @see #getLowElevationThreshold()
+ * @see #setLowElevationThreshold(double)
+ */
+ @Override
+ public > T pathDelay(final T elevation, final FieldGeodeticPoint point,
+ final T[] parameters, final FieldAbsoluteDate date) {
+
+ final Field field = date.getField();
+ final T zero = field.getZero();
+ // there are no data in the model for negative altitudes and altitude bigger than 5000 m
+ // limit the height to a range of [0, 5000] m
+ final T fixedHeight = FastMath.min(FastMath.max(zero, point.getAltitude()), zero.add(5000));
+
+ // the corrected temperature using a temperature gradient of -6.5 K/km
+ final T T = fixedHeight.multiply(6.5e-3).negate().add(t0);
+ // the corrected pressure
+ final T P = fixedHeight.multiply(2.26e-5).negate().add(1.0).pow(5.225).multiply(p0);
+ // the corrected humidity
+ final T R = FastMath.exp(fixedHeight.multiply(-6.396e-4)).multiply(r0);
+
+ // interpolate the b correction term
+ final T B = bFunction.value(fixedHeight.divide(1e3));
+ // calculate e
+ final T e = R.multiply(FastMath.exp(eFunction.value(T)));
+
+ // calculate the zenith angle from the elevation
+ final T z = FastMath.abs(FastMath.max(elevation, zero.add(lowElevationThreshold)).negate().add(zero.getPi().multiply(0.5)));
+
+ // calculate the path delay in m
+ final T tan = FastMath.tan(z);
+ final T delta = FastMath.cos(z).divide(2.277e-3).reciprocal().
+ multiply(P.add(T.divide(1255d).reciprocal().add(5e-2).multiply(e)).subtract(B.multiply(tan).multiply(tan)));
+
+ return delta;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public List getParametersDrivers() {
+ return Collections.emptyList();
+ }
+
+ /** Get the low elevation threshold value for path delay computation.
+ * @return low elevation threshold, in rad.
+ * @see #pathDelay(double, GeodeticPoint, double[], AbsoluteDate)
+ * @see #pathDelay(CalculusFieldElement, FieldGeodeticPoint, CalculusFieldElement[], FieldAbsoluteDate)
+ */
+ public double getLowElevationThreshold() {
+ return lowElevationThreshold;
+ }
+
+ /** Set the low elevation threshold value for path delay computation.
+ * @param lowElevationThreshold The new value for the threshold [rad]
+ * @see #pathDelay(double, GeodeticPoint, double[], AbsoluteDate)
+ * @see #pathDelay(CalculusFieldElement, FieldGeodeticPoint, CalculusFieldElement[], FieldAbsoluteDate)
+ */
+ public void setLowElevationThreshold(final double lowElevationThreshold) {
+ this.lowElevationThreshold = lowElevationThreshold;
+ }
+
+}
+
diff --git a/src/main/java/org/orekit/models/earth/troposphere/GiacomoDavis.java b/src/main/java/org/orekit/models/earth/troposphere/GiacomoDavis.java
new file mode 100644
index 0000000000..3eedf53240
--- /dev/null
+++ b/src/main/java/org/orekit/models/earth/troposphere/GiacomoDavis.java
@@ -0,0 +1,92 @@
+/* Copyright 2023 Thales Alenia Space
+ * Licensed to CS Communication & Systèmes (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.models.earth.troposphere;
+
+import org.hipparchus.CalculusFieldElement;
+import org.hipparchus.util.FastMath;
+
+/** Giacomo and Davis water vapor mode.
+ *
+ * See: Giacomo, P., Equation for the determination of the density of moist air, Metrologia, V. 18, 1982
+ *
+ * @author Bryan Cazabonne
+ * @author Luc Maisonobe
+ * @since 12.1
+ */
+public class GiacomoDavis implements WaterVaporPressureProvider {
+
+ /** Saturation water vapor coefficient. */
+ private static final double E = 0.01;
+
+ /** Laurent series coefficient for degree +2. */
+ private static final double L_P2 = 1.2378847e-5;
+
+ /** Laurent series coefficient for degree +1. */
+ private static final double L_P1 = -1.9121316e-2;
+
+ /** Laurent series coefficient for degree 0. */
+ private static final double L_0 = 33.93711047;
+
+ /** Laurent series coefficient for degree -1. */
+ private static final double L_M1 = -6343.1645;
+
+ /** Celsius temperature offset. */
+ private static final double CELSIUS = 273.15;
+
+ /** Constant enhancement factor. */
+ private static final double F_0 = 1.00062;
+
+ /** Pressure enhancement factor. */
+ private static final double F_P = 3.14e-6;
+
+ /** Temperature enhancement factor. */
+ private static final double F_T2 = 5.6e-7;
+
+ /** {@inheritDoc} */
+ @Override
+ public double waterVaporPressure(final double t, final double p, final double rh) {
+
+ // saturation water vapor, equation (3) of reference paper, in mbar
+ // with amended 1991 values (see reference paper)
+ final double es = FastMath.exp(t * (t * L_P2 + L_P1) + L_0 + L_M1 / t) * E;
+
+ // enhancement factor, equation (4) of reference paper
+ final double tC = t - CELSIUS;
+ final double fw = p * F_P + tC * tC * F_T2 + F_0;
+
+ return rh * fw * es;
+
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public > T waterVaporPressure(final T t, final T p, final T rh) {
+
+ // saturation water vapor, equation (3) of reference paper, in mbar
+ // with amended 1991 values (see reference paper)
+ final T es = FastMath.exp(t.multiply(t.multiply(L_P2).add(L_P1)).add(L_0).add(t.reciprocal().multiply(L_M1))).
+ multiply(E);
+
+ // enhancement factor, equation (4) of reference paper
+ final T tC = t.subtract(CELSIUS);
+ final T fw = p.multiply(F_P).add(tC.multiply(tC).multiply(F_T2)).add(F_0);
+
+ return rh.multiply(fw).multiply(es);
+
+ }
+
+}
diff --git a/src/main/java/org/orekit/models/earth/troposphere/MendesPavlisModel.java b/src/main/java/org/orekit/models/earth/troposphere/MendesPavlisModel.java
index 4644e256b2..6e21891b80 100644
--- a/src/main/java/org/orekit/models/earth/troposphere/MendesPavlisModel.java
+++ b/src/main/java/org/orekit/models/earth/troposphere/MendesPavlisModel.java
@@ -86,7 +86,7 @@ public MendesPavlisModel(final double t0, final double p0,
final double rh, final double lambda) {
this.P0 = p0;
this.T0 = t0;
- this.e0 = getWaterVapor(rh);
+ this.e0 = new GiacomoDavis().waterVaporPressure(t0, p0, rh);
this.lambda = lambda;
}
@@ -378,27 +378,4 @@ private > T computeMFCoeffient(final double a0
return point.getAltitude().multiply(a3).add(FastMath.cos(point.getLatitude()).multiply(a2)).add(a0 + a1 * T);
}
- /** Get the water vapor.
- * The water vapor model is the one of Giacomo and Davis as indicated in IERS TN 32, chap. 9.
- *
- * See: Giacomo, P., Equation for the dertermination of the density of moist air, Metrologia, V. 18, 1982
- *
- * @param rh relative humidity, in percent (50% → 0.5).
- * @return the water vapor, in mbar (1 mbar = 1 hPa).
- */
- private double getWaterVapor(final double rh) {
-
- // saturation water vapor, equation (3) of reference paper, in mbar
- // with amended 1991 values (see reference paper)
- final double es = 0.01 * FastMath.exp((1.2378847 * 1e-5) * T0 * T0 -
- (1.9121316 * 1e-2) * T0 +
- 33.93711047 -
- (6.3431645 * 1e3) * 1. / T0);
-
- // enhancement factor, equation (4) of reference paper
- final double fw = 1.00062 + (3.14 * 1e-6) * P0 + (5.6 * 1e-7) * FastMath.pow(T0 - 273.15, 2);
-
- final double e = rh * fw * es;
- return e;
- }
}
diff --git a/src/main/java/org/orekit/models/earth/troposphere/ModifiedSaastamoinenModel.java b/src/main/java/org/orekit/models/earth/troposphere/ModifiedSaastamoinenModel.java
index ab03d5ac66..fa626a6fe6 100644
--- a/src/main/java/org/orekit/models/earth/troposphere/ModifiedSaastamoinenModel.java
+++ b/src/main/java/org/orekit/models/earth/troposphere/ModifiedSaastamoinenModel.java
@@ -80,7 +80,7 @@ public class ModifiedSaastamoinenModel implements DiscreteTroposphericModel {
0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0
};
- /** E values for the B function. */
+ /** Y values for the B function. */
private static final double[] Y_VALUES_FOR_B = {
1.156, 1.079, 1.006, 0.938, 0.874, 0.813, 0.757, 0.654, 0.563
};
diff --git a/src/main/java/org/orekit/models/earth/troposphere/WaterVaporPressureProvider.java b/src/main/java/org/orekit/models/earth/troposphere/WaterVaporPressureProvider.java
new file mode 100644
index 0000000000..6d3e80471b
--- /dev/null
+++ b/src/main/java/org/orekit/models/earth/troposphere/WaterVaporPressureProvider.java
@@ -0,0 +1,45 @@
+/* Copyright 2023 Thales Alenia Space
+ * Licensed to CS Communication & Systèmes (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.models.earth.troposphere;
+
+import org.hipparchus.CalculusFieldElement;
+
+/** Interface for converting relative humidity into water vapor pressure.
+ * @author Luc Maisonobe
+ * @since 12.1
+ */
+public interface WaterVaporPressureProvider {
+
+ /** Compute water vapor pressure.
+ * @param t temperature (Kelvin)
+ * @param p pressure (hPa)
+ * @param rh relative humidity, in percent (50% → 0.5)
+ * @return water vapor pressure (hPa)
+ */
+ double waterVaporPressure(double t, double p, double rh);
+
+ /** Compute water vapor pressure.
+ * @param type of the field elements
+ * @param t temperature (Kelvin)
+ * @param p pressure (hPa)
+ * @param rh relative humidity, in percent (50% → 0.5)
+ * @return water vapor pressure (hPa)
+ */
+ > T waterVaporPressure(T t, T p, T rh);
+
+}
+
diff --git a/src/main/java/org/orekit/models/earth/weather/FieldPressureTemperatureHumidity.java b/src/main/java/org/orekit/models/earth/weather/FieldPressureTemperatureHumidity.java
new file mode 100644
index 0000000000..40003dcfcc
--- /dev/null
+++ b/src/main/java/org/orekit/models/earth/weather/FieldPressureTemperatureHumidity.java
@@ -0,0 +1,71 @@
+/* Copyright 2023 Thales Alenia Space
+ * Licensed to CS GROUP (CS) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * CS licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.orekit.models.earth.weather;
+
+import org.hipparchus.CalculusFieldElement;
+
+/** Container for pressure, temperature, and humidity.
+ * @param type of the field elements
+ * @author Luc Maisonobe
+ * @since 12.1
+ */
+public class FieldPressureTemperatureHumidity> {
+
+ /** Pressure (hPa). */
+ private final T pressure;
+
+ /** Temperature (Kelvin). */
+ private final T temperature;
+
+ /** Humidity as water vapor pressure (hPa). */
+ private final T waterVaporPressure;
+
+ /** Simple constructor.
+ * @param pressure pressure (hPa)
+ * @param temperature temperature (Kelvin)
+ * @param humidity humidity as water vapor pressure (hPa)
+ */
+ public FieldPressureTemperatureHumidity(final T pressure,
+ final T temperature,
+ final T waterVaporPressure) {
+ this.pressure = pressure;
+ this.temperature = temperature;
+ this.waterVaporPressure = waterVaporPressure;
+ }
+
+ /** Get pressure.
+ * @return pressure hPa
+ */
+ public T getPressure() {
+ return pressure;
+ }
+
+ /** Get temperature.
+ * @return temperature (Kelvin)
+ */
+ public T getTemperature() {
+ return temperature;
+ }
+
+ /** Get humidity as water vapor pressure.
+ * @return humidity as water vapor pressure (hPa)
+ */
+ public T getWaterVaporPressure() {
+ return waterVaporPressure;
+ }
+
+}
diff --git a/src/main/java/org/orekit/models/earth/weather/PressureTemperatureHumidity.java b/src/main/java/org/orekit/models/earth/weather/PressureTemperatureHumidity.java
new file mode 100644
index 0000000000..a5bf8541c5
--- /dev/null
+++ b/src/main/java/org/orekit/models/earth/weather/PressureTemperatureHumidity.java
@@ -0,0 +1,68 @@
+/* Copyright 2023 Thales Alenia Space
+ * Licensed to CS GROUP (CS) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * CS licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.orekit.models.earth.weather;
+
+/** Container for pressure, temperature, and humidity.
+ * @author Luc Maisonobe
+ * @since 12.1
+ */
+public class PressureTemperatureHumidity {
+
+ /** Pressure (hPa). */
+ private final double pressure;
+
+ /** Temperature (Kelvin). */
+ private final double temperature;
+
+ /** Humidity as water vapor pressure (hPa). */
+ private final double waterVaporPressure;
+
+ /** Simple constructor.
+ * @param pressure pressure (hPa)
+ * @param temperature temperature (Kelvin)
+ * @param humidity humidity as water vapor pressure (hPa)
+ */
+ public PressureTemperatureHumidity(final double pressure,
+ final double temperature,
+ final double waterVaporPressure) {
+ this.pressure = pressure;
+ this.temperature = temperature;
+ this.waterVaporPressure = waterVaporPressure;
+ }
+
+ /** Get pressure.
+ * @return pressure hPa
+ */
+ public double getPressure() {
+ return pressure;
+ }
+
+ /** Get temperature.
+ * @return temperature (Kelvin)
+ */
+ public double getTemperature() {
+ return temperature;
+ }
+
+ /** Get humidity as water vapor pressure.
+ * @return humidity as water vapor pressure (hPa)
+ */
+ public double getWaterVaporPressure() {
+ return waterVaporPressure;
+ }
+
+}
From 47a51cef832e42b9a79302d19b95d2516673fdc0 Mon Sep 17 00:00:00 2001
From: Luc Maisonobe
Date: Tue, 5 Dec 2023 10:55:34 +0100
Subject: [PATCH 014/359] Typo.
---
src/changes/changes.xml | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index 1cb39fe26a..ebe14aace0 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -2189,7 +2189,7 @@
view detector, measurement generation feature, possibility to use Marshall Solar Activity Future Estimation
to feed NRL MSISE 2000 atmosphere model, new tropospheric models: Mendes-Pavlis, Vienna 1, Vienna 3, estimated model,
new mapping functions for tropospheric effect: Global Mapping Function, Niell Mapping Function, Global
- Pression Temperature Models GPT and GPT2, possibility to estimate tropospheric zenith delay,
+ Pressure Temperature Models GPT and GPT2, possibility to estimate tropospheric zenith delay,
clock offset that can be estimated (both for ground station and satellite clocks).">
Added a way to manage clock corrections from GPSPropagator.
From 43ba0252a94bb435add87ac1b12ba41875e6f4b9 Mon Sep 17 00:00:00 2001
From: Serrof
Date: Tue, 5 Dec 2023 21:03:14 +0100
Subject: [PATCH 015/359] Fixed #1276: replace use of scalar multiply on Field
one
---
src/changes/changes.xml | 5 +-
.../orekit/attitudes/AttitudesSequence.java | 4 +-
.../org/orekit/attitudes/NadirPointing.java | 8 +-
.../org/orekit/bodies/PosVelChebyshev.java | 12 +--
.../EstimatedEarthFrameProvider.java | 2 +-
.../HolmesFeatherstoneAttractionModel.java | 2 +-
.../forces/gravity/NewtonianAttraction.java | 2 +-
.../ScaledConstantThrustPropulsionModel.java | 2 +-
.../radiation/KnockeRediffusedForceModel.java | 4 +-
.../java/org/orekit/frames/EOPHistory.java | 10 +--
.../orekit/frames/L1TransformProvider.java | 8 +-
.../orekit/frames/L2TransformProvider.java | 8 +-
.../org/orekit/frames/TopocentricFrame.java | 6 +-
.../java/org/orekit/frames/VEISProvider.java | 2 +-
.../org/orekit/gnss/attitude/GPSBlockIIA.java | 8 +-
.../org/orekit/gnss/attitude/GPSBlockIIF.java | 6 +-
.../org/orekit/gnss/attitude/GPSBlockIIR.java | 4 +-
.../org/orekit/gnss/attitude/Galileo.java | 2 +-
.../org/orekit/gnss/attitude/Glonass.java | 14 +--
.../java/org/orekit/models/earth/Geoid.java | 6 +-
.../models/earth/atmosphere/DTM2000.java | 12 +--
.../earth/atmosphere/HarrisPriester.java | 4 +-
.../models/earth/atmosphere/JB2008.java | 10 +--
.../models/earth/atmosphere/NRLMSISE00.java | 88 +++++++++----------
.../ionosphere/FieldNeQuickParameters.java | 31 +++----
.../earth/ionosphere/KlobucharIonoModel.java | 4 +-
.../models/earth/ionosphere/NeQuickModel.java | 6 +-
.../troposphere/FixedTroposphericDelay.java | 2 +-
.../GlobalMappingFunctionModel.java | 16 ++--
.../earth/troposphere/MendesPavlisModel.java | 2 +-
.../NiellMappingFunctionModel.java | 4 +-
.../earth/troposphere/SaastamoinenModel.java | 6 +-
.../troposphere/TroposphericModelUtils.java | 2 +-
.../earth/troposphere/ViennaOneModel.java | 24 ++---
.../earth/troposphere/ViennaThreeModel.java | 8 +-
.../orekit/orbits/FieldCartesianOrbit.java | 4 +-
.../org/orekit/orbits/FieldCircularOrbit.java | 30 +++----
.../orekit/orbits/FieldEquinoctialOrbit.java | 30 +++----
.../orbits/FieldKeplerianAnomalyUtility.java | 12 +--
.../orekit/orbits/FieldKeplerianOrbit.java | 22 ++---
.../propagation/FieldSpacecraftState.java | 24 ++---
.../FieldSpacecraftStateInterpolator.java | 2 +-
.../BrouwerLyddaneGradientConverter.java | 2 +-
.../EcksteinHechlerGradientConverter.java | 2 +-
.../FieldBrouwerLyddanePropagator.java | 2 +-
.../FieldEcksteinHechlerPropagator.java | 6 +-
.../analytical/FieldKeplerianPropagator.java | 9 +-
.../KeplerianGradientConverter.java | 2 +-
.../gnss/GLONASSAnalyticalPropagator.java | 10 +--
.../analytical/tle/FieldDeepSDP4.java | 20 ++---
.../propagation/analytical/tle/FieldSGP4.java | 2 +-
.../propagation/analytical/tle/FieldTLE.java | 10 +--
.../analytical/tle/FieldTLEPropagator.java | 10 +--
.../FixedPointTleGenerationAlgorithm.java | 2 +-
.../LeastSquaresTleGenerationAlgorithm.java | 2 +-
...stractFixedStepFieldIntegratorBuilder.java | 2 +-
.../events/FieldAltitudeDetector.java | 6 +-
.../propagation/events/FieldDateDetector.java | 2 +-
.../events/FieldElevationDetector.java | 2 +-
.../propagation/events/FieldEventState.java | 4 +-
.../events/FieldFunctionalDetector.java | 2 +-
.../events/FieldLatitudeCrossingDetector.java | 2 +-
.../FieldLongitudeCrossingDetector.java | 2 +-
.../propagation/events/FieldNodeDetector.java | 2 +-
.../numerical/FieldNumericalPropagator.java | 4 +-
.../numerical/cr3bp/CR3BPForceModel.java | 4 +-
.../dsst/FieldDSSTPropagator.java | 6 +-
.../forces/AbstractGaussianContribution.java | 6 +-
.../forces/DSSTSolarRadiationPressure.java | 6 +-
.../dsst/forces/DSSTTesseral.java | 10 +--
.../dsst/forces/DSSTThirdBody.java | 16 ++--
.../semianalytical/dsst/forces/DSSTZonal.java | 16 ++--
.../dsst/forces/FieldDSSTTesseralContext.java | 2 +-
.../dsst/utilities/CoefficientsFactory.java | 2 +-
.../dsst/utilities/FieldCjSjCoefficient.java | 2 +-
.../dsst/utilities/FieldLnsCoefficients.java | 2 +-
.../dsst/utilities/UpperBounds.java | 8 +-
.../hansen/FieldHansenThirdBodyLinear.java | 2 +-
.../FieldShortTermEncounter2DDefinition.java | 4 +-
.../probability/twod/Laas2015.java | 10 +--
.../probability/twod/Patera2005.java | 2 +-
src/main/java/org/orekit/time/BDTScale.java | 2 +-
src/main/java/org/orekit/time/GPSScale.java | 2 +-
.../java/org/orekit/time/GalileoScale.java | 2 +-
src/main/java/org/orekit/time/IRNSSScale.java | 2 +-
src/main/java/org/orekit/time/QZSSScale.java | 2 +-
src/main/java/org/orekit/time/TTScale.java | 2 +-
src/main/java/org/orekit/time/UTCScale.java | 2 +-
.../java/org/orekit/time/UTCTAIOffset.java | 2 +-
.../orekit/utils/FieldAngularCoordinates.java | 4 +-
.../java/org/orekit/utils/Fieldifier.java | 76 ++++++++--------
.../org/orekit/utils/IERSConventions.java | 4 +-
.../utils/ParameterDriversProvider.java | 6 +-
.../TimeStampedFieldAngularCoordinates.java | 2 +-
...AngularCoordinatesHermiteInterpolator.java | 2 +-
.../attitudes/FrameAlignedProviderTest.java | 2 +-
96 files changed, 389 insertions(+), 386 deletions(-)
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index 0481c0d1f2..de1a4eb91a 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -20,7 +20,10 @@
Orekit Changes
-
+
+
+ Removed uses of scalar multiply on Field one.
+
Added utility classes for position angle conversions for (Field) CircularOrbit and EquinoctialOrbit.
diff --git a/src/main/java/org/orekit/attitudes/AttitudesSequence.java b/src/main/java/org/orekit/attitudes/AttitudesSequence.java
index 3daab4b3dc..6001b285b1 100644
--- a/src/main/java/org/orekit/attitudes/AttitudesSequence.java
+++ b/src/main/java/org/orekit/attitudes/AttitudesSequence.java
@@ -157,13 +157,13 @@ public void init(final FieldSpacecraftState s0,
/** {@inheritDoc} */
@Override
public T g(final FieldSpacecraftState s) {
- return field.getZero().add(sw.g(s.toSpacecraftState()));
+ return field.getZero().newInstance(sw.g(s.toSpacecraftState()));
}
/** {@inheritDoc} */
@Override
public T getThreshold() {
- return field.getZero().add(sw.getThreshold());
+ return field.getZero().newInstance(sw.getThreshold());
}
/** {@inheritDoc} */
diff --git a/src/main/java/org/orekit/attitudes/NadirPointing.java b/src/main/java/org/orekit/attitudes/NadirPointing.java
index 9f28620005..56923da920 100644
--- a/src/main/java/org/orekit/attitudes/NadirPointing.java
+++ b/src/main/java/org/orekit/attitudes/NadirPointing.java
@@ -109,11 +109,11 @@ public > TimeStampedFieldPVCoordinates getT
// sample intersection points in current date neighborhood
final double h = 0.01;
final List> sample = new ArrayList<>();
- sample.add(nadirRef(pvProv.getPVCoordinates(date.shiftedBy(-2 * h), frame), refToBody.staticShiftedBy(zero.add(-2 * h))));
- sample.add(nadirRef(pvProv.getPVCoordinates(date.shiftedBy(-h), frame), refToBody.staticShiftedBy(zero.add(-h))));
+ sample.add(nadirRef(pvProv.getPVCoordinates(date.shiftedBy(-2 * h), frame), refToBody.staticShiftedBy(zero.newInstance(-2 * h))));
+ sample.add(nadirRef(pvProv.getPVCoordinates(date.shiftedBy(-h), frame), refToBody.staticShiftedBy(zero.newInstance(-h))));
sample.add(nadirRef(pvProv.getPVCoordinates(date, frame), refToBody));
- sample.add(nadirRef(pvProv.getPVCoordinates(date.shiftedBy(+h), frame), refToBody.staticShiftedBy(zero.add(+h))));
- sample.add(nadirRef(pvProv.getPVCoordinates(date.shiftedBy(+2 * h), frame), refToBody.staticShiftedBy(zero.add(+2 * h))));
+ sample.add(nadirRef(pvProv.getPVCoordinates(date.shiftedBy(+h), frame), refToBody.staticShiftedBy(zero.newInstance(+h))));
+ sample.add(nadirRef(pvProv.getPVCoordinates(date.shiftedBy(+2 * h), frame), refToBody.staticShiftedBy(zero.newInstance(+2 * h))));
// create interpolator
final FieldTimeInterpolator, T> interpolator =
diff --git a/src/main/java/org/orekit/bodies/PosVelChebyshev.java b/src/main/java/org/orekit/bodies/PosVelChebyshev.java
index 1451eeba7f..2c3b7fc120 100644
--- a/src/main/java/org/orekit/bodies/PosVelChebyshev.java
+++ b/src/main/java/org/orekit/bodies/PosVelChebyshev.java
@@ -169,9 +169,9 @@ > FieldVector3D getPosition(final FieldAbso
// initialize Chebyshev polynomials recursion
T pKm1 = one;
T pK = t;
- T xP = zero.add(xCoeffs[0]);
- T yP = zero.add(yCoeffs[0]);
- T zP = zero.add(zCoeffs[0]);
+ T xP = zero.newInstance(xCoeffs[0]);
+ T yP = zero.newInstance(yCoeffs[0]);
+ T zP = zero.newInstance(zCoeffs[0]);
// combine polynomials by applying coefficients
for (int k = 1; k < xCoeffs.length; ++k) {
@@ -281,9 +281,9 @@ > FieldPVCoordinates getPositionVelocityAcc
// initialize Chebyshev polynomials recursion
T pKm1 = one;
T pK = t;
- T xP = zero.add(xCoeffs[0]);
- T yP = zero.add(yCoeffs[0]);
- T zP = zero.add(zCoeffs[0]);
+ T xP = zero.newInstance(xCoeffs[0]);
+ T yP = zero.newInstance(yCoeffs[0]);
+ T zP = zero.newInstance(zCoeffs[0]);
// initialize Chebyshev polynomials derivatives recursion
T qKm1 = zero;
diff --git a/src/main/java/org/orekit/estimation/measurements/EstimatedEarthFrameProvider.java b/src/main/java/org/orekit/estimation/measurements/EstimatedEarthFrameProvider.java
index 9222f17365..d399c22270 100644
--- a/src/main/java/org/orekit/estimation/measurements/EstimatedEarthFrameProvider.java
+++ b/src/main/java/org/orekit/estimation/measurements/EstimatedEarthFrameProvider.java
@@ -279,7 +279,7 @@ public > FieldTransform getTransform(final
// prime meridian shift parameters
final T theta = linearModel(date, primeMeridianOffsetDriver, primeMeridianDriftDriver);
- final T thetaDot = zero.add(primeMeridianDriftDriver.getValue());
+ final T thetaDot = zero.newInstance(primeMeridianDriftDriver.getValue());
// pole shift parameters
final T xpNeg = linearModel(date, polarOffsetXDriver, polarDriftXDriver).negate();
diff --git a/src/main/java/org/orekit/forces/gravity/HolmesFeatherstoneAttractionModel.java b/src/main/java/org/orekit/forces/gravity/HolmesFeatherstoneAttractionModel.java
index 21bf1c3aff..03809eb4ed 100644
--- a/src/main/java/org/orekit/forces/gravity/HolmesFeatherstoneAttractionModel.java
+++ b/src/main/java/org/orekit/forces/gravity/HolmesFeatherstoneAttractionModel.java
@@ -977,7 +977,7 @@ private > int computeTesseral(final int m, fin
// initialize recursion from sectorial terms
int n = FastMath.max(2, m);
if (n == m) {
- pnm0[n] = zero.add(sectorial[n]);
+ pnm0[n] = zero.newInstance(sectorial[n]);
++n;
}
diff --git a/src/main/java/org/orekit/forces/gravity/NewtonianAttraction.java b/src/main/java/org/orekit/forces/gravity/NewtonianAttraction.java
index 6a891e860b..bba1b4c458 100644
--- a/src/main/java/org/orekit/forces/gravity/NewtonianAttraction.java
+++ b/src/main/java/org/orekit/forces/gravity/NewtonianAttraction.java
@@ -83,7 +83,7 @@ public double getMu(final AbsoluteDate date) {
*/
public > T getMu(final Field field, final FieldAbsoluteDate date) {
final T zero = field.getZero();
- return zero.add(gmParameterDriver.getValue(date.toAbsoluteDate()));
+ return zero.newInstance(gmParameterDriver.getValue(date.toAbsoluteDate()));
}
/** {@inheritDoc} */
diff --git a/src/main/java/org/orekit/forces/maneuvers/propulsion/ScaledConstantThrustPropulsionModel.java b/src/main/java/org/orekit/forces/maneuvers/propulsion/ScaledConstantThrustPropulsionModel.java
index 2f0242340e..999ae687a0 100644
--- a/src/main/java/org/orekit/forces/maneuvers/propulsion/ScaledConstantThrustPropulsionModel.java
+++ b/src/main/java/org/orekit/forces/maneuvers/propulsion/ScaledConstantThrustPropulsionModel.java
@@ -154,6 +154,6 @@ public > FieldVector3D getThrustVector(fina
/** {@inheritDoc} */
@Override
public > T getFlowRate(final T[] parameters) {
- return parameters[0].getField().getZero().add(getInitialFlowRate());
+ return parameters[0].getField().getZero().newInstance(getInitialFlowRate());
}
}
diff --git a/src/main/java/org/orekit/forces/radiation/KnockeRediffusedForceModel.java b/src/main/java/org/orekit/forces/radiation/KnockeRediffusedForceModel.java
index 756cde4a7f..870495a304 100644
--- a/src/main/java/org/orekit/forces/radiation/KnockeRediffusedForceModel.java
+++ b/src/main/java/org/orekit/forces/radiation/KnockeRediffusedForceModel.java
@@ -293,8 +293,8 @@ public > FieldVector3D acceleration(final F
// Compute current elementary crown sector area, it results of the integration of an elementary crown sector
// over the angular resolution
- final T sectorArea = zero.add(equatorialRadius * equatorialRadius *
- 2.0 * angularResolution * FastMath.sin(0.5 * angularResolution) * FastMath.sin(eastAxisOffset));
+ final T sectorArea = zero.newInstance(equatorialRadius * equatorialRadius *
+ 2.0 * angularResolution * FastMath.sin(0.5 * angularResolution) * FastMath.sin(eastAxisOffset));
// Add current sector contribution to total rediffused flux
rediffusedFlux = rediffusedFlux.add(computeElementaryFlux(s, currentCenter, sunPosition, earth, sectorArea));
diff --git a/src/main/java/org/orekit/frames/EOPHistory.java b/src/main/java/org/orekit/frames/EOPHistory.java
index cf787b1630..e93264a491 100644
--- a/src/main/java/org/orekit/frames/EOPHistory.java
+++ b/src/main/java/org/orekit/frames/EOPHistory.java
@@ -405,7 +405,7 @@ public void accept(final EOPEntry neighbor) {
dut = neighbor.getUT1MinusUTC();
}
final T[] array = MathArrays.buildArray(date.getField(), 1);
- array[0] = date.getField().getZero().add(dut);
+ array[0] = date.getField().getZero().newInstance(dut);
interpolator.addSamplePoint(date.durationFrom(neighbor.getDate()).negate(),
array);
}
@@ -745,7 +745,7 @@ private > T interpolate(final FieldAbsoluteDat
// if T was DerivativeStructure
getNeighbors(aDate, interpolationDegree + 1).
forEach(entry -> {
- y[0] = zero.add(selector.apply(entry));
+ y[0] = zero.newInstance(selector.apply(entry));
interpolator.addSamplePoint(central.durationFrom(entry.getDate()).negate(), y);
});
return interpolator.value(date.durationFrom(central))[0]; // here, we introduce derivatives again (in DerivativeStructure case)
@@ -842,8 +842,8 @@ private > T[] interpolate(final FieldAbsoluteD
// if T was DerivativeStructure
getNeighbors(aDate, interpolationDegree + 1).
forEach(entry -> {
- y[0] = zero.add(selector1.apply(entry));
- y[1] = zero.add(selector2.apply(entry));
+ y[0] = zero.newInstance(selector1.apply(entry));
+ y[1] = zero.newInstance(selector2.apply(entry));
interpolator.addSamplePoint(central.durationFrom(entry.getDate()).negate(), y);
});
return interpolator.value(date.durationFrom(central)); // here, we introduce derivatives again (in DerivativeStructure case)
@@ -1052,7 +1052,7 @@ public > T[] value(final FieldAbsoluteDate
// if T was DerivativeStructure
cache.getNeighbors(aDate).forEach(entry -> {
for (int i = 0; i < y.length; ++i) {
- y[i] = zero.add(entry.correction[i]);
+ y[i] = zero.newInstance(entry.correction[i]);
}
interpolator.addSamplePoint(central.durationFrom(entry.getDate()).negate(), y);
});
diff --git a/src/main/java/org/orekit/frames/L1TransformProvider.java b/src/main/java/org/orekit/frames/L1TransformProvider.java
index 1e6710772c..78d5f0d816 100644
--- a/src/main/java/org/orekit/frames/L1TransformProvider.java
+++ b/src/main/java/org/orekit/frames/L1TransformProvider.java
@@ -195,13 +195,13 @@ private Vector3D getL1(final Vector3D primaryToSecondary) {
final T zero = primaryToSecondary.getX().getField().getZero();
final T[] searchInterval = UnivariateSolverUtils.bracket(l1Equation,
baseR, zero, bigR.multiply(2),
- bigR.multiply(0.01), zero.add(1),
+ bigR.multiply(0.01), zero.newInstance(1),
MAX_EVALUATIONS);
- final T relativeAccuracy = zero.add(RELATIVE_ACCURACY);
- final T absoluteAccuracy = zero.add(ABSOLUTE_ACCURACY);
- final T functionAccuracy = zero.add(FUNCTION_ACCURACY);
+ final T relativeAccuracy = zero.newInstance(RELATIVE_ACCURACY);
+ final T absoluteAccuracy = zero.newInstance(ABSOLUTE_ACCURACY);
+ final T functionAccuracy = zero.newInstance(FUNCTION_ACCURACY);
final FieldBracketingNthOrderBrentSolver solver =
new FieldBracketingNthOrderBrentSolver<>(relativeAccuracy,
diff --git a/src/main/java/org/orekit/frames/L2TransformProvider.java b/src/main/java/org/orekit/frames/L2TransformProvider.java
index f36b28141f..1364e0311b 100644
--- a/src/main/java/org/orekit/frames/L2TransformProvider.java
+++ b/src/main/java/org/orekit/frames/L2TransformProvider.java
@@ -195,13 +195,13 @@ private Vector3D getL2(final Vector3D primaryToSecondary) {
final T zero = primaryToSecondary.getX().getField().getZero();
final T[] searchInterval = UnivariateSolverUtils.bracket(l2Equation,
baseR, zero, bigR.multiply(2),
- bigR.multiply(0.01), zero.add(1),
+ bigR.multiply(0.01), zero,
MAX_EVALUATIONS);
- final T relativeAccuracy = zero.add(RELATIVE_ACCURACY);
- final T absoluteAccuracy = zero.add(ABSOLUTE_ACCURACY);
- final T functionAccuracy = zero.add(FUNCTION_ACCURACY);
+ final T relativeAccuracy = zero.newInstance(RELATIVE_ACCURACY);
+ final T absoluteAccuracy = zero.newInstance(ABSOLUTE_ACCURACY);
+ final T functionAccuracy = zero.newInstance(FUNCTION_ACCURACY);
final FieldBracketingNthOrderBrentSolver solver =
new FieldBracketingNthOrderBrentSolver<>(relativeAccuracy,
diff --git a/src/main/java/org/orekit/frames/TopocentricFrame.java b/src/main/java/org/orekit/frames/TopocentricFrame.java
index 4722bd37e7..ccb84ae3f5 100644
--- a/src/main/java/org/orekit/frames/TopocentricFrame.java
+++ b/src/main/java/org/orekit/frames/TopocentricFrame.java
@@ -129,9 +129,9 @@ public Vector3D getCartesianPoint() {
*/
public > FieldGeodeticPoint getPoint(final Field field) {
final T zero = field.getZero();
- return new FieldGeodeticPoint<>(zero.add(point.getLatitude()),
- zero.add(point.getLongitude()),
- zero.add(point.getAltitude()));
+ return new FieldGeodeticPoint<>(zero.newInstance(point.getLatitude()),
+ zero.newInstance(point.getLongitude()),
+ zero.newInstance(point.getAltitude()));
}
/** Get the zenith direction of topocentric frame, expressed in parent shape frame.
diff --git a/src/main/java/org/orekit/frames/VEISProvider.java b/src/main/java/org/orekit/frames/VEISProvider.java
index 7d192d00b5..4a98db53b6 100644
--- a/src/main/java/org/orekit/frames/VEISProvider.java
+++ b/src/main/java/org/orekit/frames/VEISProvider.java
@@ -108,7 +108,7 @@ public > FieldTransform getTransform(final
final T vst = ttd.multiply(VST1).add(rdtt.multiply(MathUtils.TWO_PI)).add(VST0).remainder(MathUtils.TWO_PI);
// compute angular rotation of Earth, in rad/s
- final FieldVector3D rotationRate = new FieldVector3D<>(date.getField().getZero().add(-VSTD),
+ final FieldVector3D rotationRate = new FieldVector3D<>(date.getField().getZero().newInstance(-VSTD),
Vector3D.PLUS_K);
// set up the transform from parent GTOD
diff --git a/src/main/java/org/orekit/gnss/attitude/GPSBlockIIA.java b/src/main/java/org/orekit/gnss/attitude/GPSBlockIIA.java
index eb71b14e22..29676edb9a 100644
--- a/src/main/java/org/orekit/gnss/attitude/GPSBlockIIA.java
+++ b/src/main/java/org/orekit/gnss/attitude/GPSBlockIIA.java
@@ -159,7 +159,7 @@ protected > TimeStampedFieldAngularCoordinates
// noon beta angle limit from yaw rate
final T aNoon = FastMath.atan(context.getMuRate().divide(yawRate));
- final T aNight = field.getZero().add(NIGHT_TURN_LIMIT);
+ final T aNight = field.getZero().newInstance(NIGHT_TURN_LIMIT);
final double cNoon = FastMath.cos(aNoon.getReal());
final double cNight = FastMath.cos(aNight.getReal());
@@ -183,16 +183,16 @@ protected > TimeStampedFieldAngularCoordinates
if (beta.getReal() > 0 && beta.getReal() < yawBias) {
// noon turn problem for small positive beta in block IIA
// rotation is in the wrong direction for these spacecrafts
- phiDot = field.getZero().add(FastMath.copySign(yawRate, beta.getReal()));
+ phiDot = field.getZero().newInstance(FastMath.copySign(yawRate, beta.getReal()));
linearPhi = phiStart.add(phiDot.multiply(dtStart));
} else {
// regular noon turn
- phiDot = field.getZero().add(-FastMath.copySign(yawRate, beta.getReal()));
+ phiDot = field.getZero().newInstance(-FastMath.copySign(yawRate, beta.getReal()));
linearPhi = phiStart.add(phiDot.multiply(dtStart));
}
} else {
// midnight turn
- phiDot = field.getZero().add(yawRate);
+ phiDot = field.getZero().newInstance(yawRate);
linearPhi = phiStart.add(phiDot.multiply(dtStart));
}
diff --git a/src/main/java/org/orekit/gnss/attitude/GPSBlockIIF.java b/src/main/java/org/orekit/gnss/attitude/GPSBlockIIF.java
index 15bda16b27..67a3056660 100644
--- a/src/main/java/org/orekit/gnss/attitude/GPSBlockIIF.java
+++ b/src/main/java/org/orekit/gnss/attitude/GPSBlockIIF.java
@@ -140,7 +140,7 @@ protected > TimeStampedFieldAngularCoordinates
// noon beta angle limit from yaw rate
final T aNoon = FastMath.atan(context.getMuRate().divide(yawRate));
- final T aNight = field.getZero().add(NIGHT_TURN_LIMIT);
+ final T aNight = field.getZero().newInstance(NIGHT_TURN_LIMIT);
final double cNoon = FastMath.cos(aNoon.getReal());
final double cNight = FastMath.cos(aNight.getReal());
@@ -164,11 +164,11 @@ protected > TimeStampedFieldAngularCoordinates
if (beta.getReal() > yawBias && beta.getReal() < 0) {
// noon turn problem for small negative beta in block IIF
// rotation is in the wrong direction for these spacecrafts
- phiDot = field.getZero().add(FastMath.copySign(yawRate, beta.getReal()));
+ phiDot = field.getZero().newInstance(FastMath.copySign(yawRate, beta.getReal()));
linearPhi = phiStart.add(phiDot.multiply(dtStart));
} else {
// regular noon turn
- phiDot = field.getZero().add(-FastMath.copySign(yawRate, beta.getReal()));
+ phiDot = field.getZero().newInstance(-FastMath.copySign(yawRate, beta.getReal()));
linearPhi = phiStart.add(phiDot.multiply(dtStart));
}
} else {
diff --git a/src/main/java/org/orekit/gnss/attitude/GPSBlockIIR.java b/src/main/java/org/orekit/gnss/attitude/GPSBlockIIR.java
index 736e2908db..9d5e49ff02 100644
--- a/src/main/java/org/orekit/gnss/attitude/GPSBlockIIR.java
+++ b/src/main/java/org/orekit/gnss/attitude/GPSBlockIIR.java
@@ -135,11 +135,11 @@ protected > TimeStampedFieldAngularCoordinates
if (context.inSunSide()) {
// noon turn
- phiDot = field.getZero().add(-FastMath.copySign(yawRate, beta.getReal()));
+ phiDot = field.getZero().newInstance(-FastMath.copySign(yawRate, beta.getReal()));
linearPhi = phiStart.add(phiDot.multiply(dtStart));
} else {
// midnight turn
- phiDot = field.getZero().add(FastMath.copySign(yawRate, beta.getReal()));
+ phiDot = field.getZero().newInstance(FastMath.copySign(yawRate, beta.getReal()));
linearPhi = phiStart.add(phiDot.multiply(dtStart));
}
diff --git a/src/main/java/org/orekit/gnss/attitude/Galileo.java b/src/main/java/org/orekit/gnss/attitude/Galileo.java
index e5f3fd0a1b..8c7da5b94b 100644
--- a/src/main/java/org/orekit/gnss/attitude/Galileo.java
+++ b/src/main/java/org/orekit/gnss/attitude/Galileo.java
@@ -127,7 +127,7 @@ protected > TimeStampedFieldAngularCoordinates
context.setUpTurnRegion(COS_NIGHT, COS_NOON)) {
final Field field = context.getDate().getField();
- final T betaX = field.getZero().add(BETA_X);
+ final T betaX = field.getZero().newInstance(BETA_X);
context.setHalfSpan(context.inSunSide() ?
betaX :
context.inOrbitPlaneAbsoluteAngle(betaX),
diff --git a/src/main/java/org/orekit/gnss/attitude/Glonass.java b/src/main/java/org/orekit/gnss/attitude/Glonass.java
index bcdf39dcc8..5d4c46d9e3 100644
--- a/src/main/java/org/orekit/gnss/attitude/Glonass.java
+++ b/src/main/java/org/orekit/gnss/attitude/Glonass.java
@@ -152,7 +152,7 @@ protected > TimeStampedFieldAngularCoordinates
// noon beta angle limit from yaw rate
final T realBeta = context.beta(context.getDate());
final T muRate = context.getMuRate();
- final T aNight = field.getZero().add(NIGHT_TURN_LIMIT);
+ final T aNight = field.getZero().newInstance(NIGHT_TURN_LIMIT);
T aNoon = FastMath.atan(muRate.divide(yawRate));
if (FastMath.abs(realBeta).getReal() < aNoon.getReal()) {
final CalculusFieldUnivariateFunction f = yawEnd -> {
@@ -161,11 +161,11 @@ protected > TimeStampedFieldAngularCoordinates
subtract(context.computePhi(realBeta, delta.negate()))).
multiply(0.5));
};
- final T[] bracket = UnivariateSolverUtils.bracket(f, field.getZero().add(YAW_END_ZERO),
+ final T[] bracket = UnivariateSolverUtils.bracket(f, field.getZero().newInstance(YAW_END_ZERO),
field.getZero(), field.getZero().getPi());
- final T yawEnd = new FieldBracketingNthOrderBrentSolver<>(field.getZero().add(1.0e-14),
- field.getZero().add(1.0e-8),
- field.getZero().add(1.0e-15),
+ final T yawEnd = new FieldBracketingNthOrderBrentSolver<>(field.getZero().newInstance(1.0e-14),
+ field.getZero().newInstance(1.0e-8),
+ field.getZero().newInstance(1.0e-15),
5).
solve(50, f, bracket[0], bracket[1], AllowedSolution.ANY_SIDE);
aNoon = muRate.multiply(yawEnd).divide(yawRate);
@@ -192,11 +192,11 @@ protected > TimeStampedFieldAngularCoordinates
final T phiEnd = context.getYawEnd(beta);
if (context.inSunSide()) {
// noon turn
- phiDot = field.getZero().add(-FastMath.copySign(yawRate, beta.getReal()));
+ phiDot = field.getZero().newInstance(-FastMath.copySign(yawRate, beta.getReal()));
linearPhi = phiStart.add(phiDot.multiply(dtStart));
} else {
// midnight turn
- phiDot = field.getZero().add(FastMath.copySign(yawRate, beta.getReal()));
+ phiDot = field.getZero().newInstance(FastMath.copySign(yawRate, beta.getReal()));
linearPhi = phiStart.add(phiDot.multiply(dtStart));
// this turn limitation is only computed for midnight turns in Kouba model
diff --git a/src/main/java/org/orekit/models/earth/Geoid.java b/src/main/java/org/orekit/models/earth/Geoid.java
index bde307050a..0dacf98d11 100644
--- a/src/main/java/org/orekit/models/earth/Geoid.java
+++ b/src/main/java/org/orekit/models/earth/Geoid.java
@@ -520,9 +520,9 @@ public > FieldGeodeticPoint getIntersection
}
// solve line search problem to find the intersection
final FieldBracketingNthOrderBrentSolver solver =
- new FieldBracketingNthOrderBrentSolver<>(field.getZero().add(1.0e-14),
- field.getZero().add(1.0e-6),
- field.getZero().add(1.0e-15),
+ new FieldBracketingNthOrderBrentSolver<>(field.getZero().newInstance(1.0e-14),
+ field.getZero().newInstance(1.0e-6),
+ field.getZero().newInstance(1.0e-15),
5);
try {
final T abscissa = solver.solve(MAX_EVALUATIONS, heightFunction, lowPoint, highPoint,
diff --git a/src/main/java/org/orekit/models/earth/atmosphere/DTM2000.java b/src/main/java/org/orekit/models/earth/atmosphere/DTM2000.java
index 8102af1d65..4f515f57cf 100644
--- a/src/main/java/org/orekit/models/earth/atmosphere/DTM2000.java
+++ b/src/main/java/org/orekit/models/earth/atmosphere/DTM2000.java
@@ -1154,8 +1154,8 @@ private T gFunction(final double[] a, final T[] da,
fmfb[2] = f[2] - fbar[2];
fbm150[1] = fbar[1] - 150.0;
fbm150[2] = fbar[2];
- da[4] = zero.add(fmfb[1]);
- da[6] = zero.add(fbm150[1]);
+ da[4] = zero.newInstance(fmfb[1]);
+ da[6] = zero.newInstance(fbm150[1]);
da[4] = da[4].add(a[70] * fmfb[2]);
da[6] = da[6].add(a[71] * fbm150[2]);
da[70] = da[4].multiply(a[ 5]).multiply(2).
@@ -1200,10 +1200,10 @@ private T gFunction(final double[] a, final T[] da,
da[60] = dkp.multiply(dkp);
da[61] = p20mg.multiply(da[60]);
da[75] = da[60].multiply(da[60]);
- da[64] = zero.add(dkpm);
+ da[64] = zero.newInstance(dkpm);
da[65] = p20mg.multiply(dkpm);
da[72] = p40mg.multiply(dkpm);
- da[66] = zero.add(dkpm * dkpm);
+ da[66] = zero.newInstance(dkpm * dkpm);
da[73] = p20mg.multiply(da[66]);
da[76] = da[66].multiply(da[66]);
@@ -1240,10 +1240,10 @@ private T gFunction(final double[] a, final T[] da,
add(da[78].multiply(a78)).
add(da[79].multiply(a[79]));
// termes annuels symetriques en latitude
- da[9] = zero.add(FastMath.cos(ROT * (day - a[11])));
+ da[9] = zero.newInstance(FastMath.cos(ROT * (day - a[11])));
da[10] = p20.multiply(da[9]);
// termes semi-annuels symetriques en latitude
- da[12] = zero.add(FastMath.cos(ROT2 * (day - a[14])));
+ da[12] = zero.newInstance(FastMath.cos(ROT2 * (day - a[14])));
da[13] = p20.multiply(da[12]);
// termes annuels non symetriques en latitude
final double coste = FastMath.cos(ROT * (day - a[18]));
diff --git a/src/main/java/org/orekit/models/earth/atmosphere/HarrisPriester.java b/src/main/java/org/orekit/models/earth/atmosphere/HarrisPriester.java
index f7ae30c8c6..78eaf2d834 100644
--- a/src/main/java/org/orekit/models/earth/atmosphere/HarrisPriester.java
+++ b/src/main/java/org/orekit/models/earth/atmosphere/HarrisPriester.java
@@ -361,13 +361,13 @@ public > T getDensity(final Vector3D sunInEart
final T dH = posAlt.negate().add(tabAltRho[ia][0]).divide(tabAltRho[ia][0] - tabAltRho[ia + 1][0]);
// Min exponential density interpolation
- final T rhoMin = zero.add(tabAltRho[ia + 1][1] / tabAltRho[ia][1]).pow(dH).multiply(tabAltRho[ia][1]);
+ final T rhoMin = zero.newInstance(tabAltRho[ia + 1][1] / tabAltRho[ia][1]).pow(dH).multiply(tabAltRho[ia][1]);
if (Precision.equals(cosPow.getReal(), 0.)) {
return zero.add(rhoMin);
} else {
// Max exponential density interpolation
- final T rhoMax = zero.add(tabAltRho[ia + 1][2] / tabAltRho[ia][2]).pow(dH).multiply(tabAltRho[ia][2]);
+ final T rhoMax = zero.newInstance(tabAltRho[ia + 1][2] / tabAltRho[ia][2]).pow(dH).multiply(tabAltRho[ia][2]);
return rhoMin.add(rhoMax.subtract(rhoMin).multiply(cosPow));
}
diff --git a/src/main/java/org/orekit/models/earth/atmosphere/JB2008.java b/src/main/java/org/orekit/models/earth/atmosphere/JB2008.java
index 5cb45077b5..5bfd618fda 100644
--- a/src/main/java/org/orekit/models/earth/atmosphere/JB2008.java
+++ b/src/main/java/org/orekit/models/earth/atmosphere/JB2008.java
@@ -576,7 +576,7 @@ public > T getDensity(final T dateMJD, final T
tc[3] = gsubx.divide(tc[2]);
// Equation (5)
- final T z1 = field.getZero().add(90.);
+ final T z1 = field.getZero().newInstance(90.);
final T z2 = min(105.0, altKm);
T al = z2.divide(z1).log();
int n = 1 + (int) FastMath.floor(al.getReal() / R1);
@@ -826,7 +826,7 @@ private static > T dTc(final double f10, final
final T h = satAlt.subtract(200.0).divide(50.0);
dTc = poly1CDTC(fs, st, cs).multiply(h).add(poly2CDTC(fs, st, cs));
} else if (satAlt.getReal() > 240.0 && satAlt.getReal() <= 300.0) {
- final T h = solarTime.getField().getZero().add(0.8);
+ final T h = solarTime.getField().getZero().newInstance(0.8);
final T bb = poly1CDTC(fs, st, cs);
final T aa = bb.multiply(h).add(poly2CDTC(fs, st, cs));
final T p2BDT = poly2BDTC(st);
@@ -998,7 +998,7 @@ private static double mBar(final double z) {
*/
private static > T mBar(final T z) {
final T dz = z.subtract(100.);
- T amb = z.getField().getZero().add(CMB[6]);
+ T amb = z.getField().getZero().newInstance(CMB[6]);
for (int i = 5; i >= 0; --i) {
amb = dz.multiply(amb).add(CMB[i]);
}
@@ -1182,7 +1182,7 @@ private static > T dayOfYear(final T dateMJD)
* @return min value
*/
private > T min(final double d, final T f) {
- return (f.getReal() > d) ? f.getField().getZero().add(d) : f;
+ return (f.getReal() > d) ? f.getField().getZero().newInstance(d) : f;
}
/** Compute max of two values, one double and one field element.
@@ -1192,7 +1192,7 @@ private > T min(final double d, final T f) {
* @return max value
*/
private > T max(final double d, final T f) {
- return (f.getReal() <= d) ? f.getField().getZero().add(d) : f;
+ return (f.getReal() <= d) ? f.getField().getZero().newInstance(d) : f;
}
/** Get the local density.
diff --git a/src/main/java/org/orekit/models/earth/atmosphere/NRLMSISE00.java b/src/main/java/org/orekit/models/earth/atmosphere/NRLMSISE00.java
index bd429b6479..82a37ad305 100644
--- a/src/main/java/org/orekit/models/earth/atmosphere/NRLMSISE00.java
+++ b/src/main/java/org/orekit/models/earth/atmosphere/NRLMSISE00.java
@@ -2825,7 +2825,7 @@ public class FieldOutput> {
temperatures = MathArrays.buildArray(field, 2);
// Calculates latitude variable gravity and effective radius
- final T xlat = (sw[2] == 0) ? zero.add(LAT_REF) : lat;
+ final T xlat = (sw[2] == 0) ? zero.newInstance(LAT_REF) : lat;
final T c2 = xlat.multiply(2 * DEG_TO_RAD).cos();
glat = c2.multiply(-0.0026373).add(1).multiply(G_REF);
rlat = glat.multiply(2).divide(c2.multiply(2.27e-9).add(3.085462e-6)).multiply(1.e-5);
@@ -2921,7 +2921,7 @@ void gts7(final T alt) {
final double xmm = PDM[2][4];
/**** Exospheric temperature ****/
- T tinf = zero.add(PTM[0] * PT[0]);
+ T tinf = zero.newInstance(PTM[0] * PT[0]);
// Tinf variations not important below ZA or ZN[0]
if (alt.getReal() > ZN1[0]) {
tinf = tinf.multiply(globe7(PT).multiply(sw[16]).add(1));
@@ -2929,24 +2929,24 @@ void gts7(final T alt) {
setTemperature(EXOSPHERIC, tinf);
// Gradient variations not important below ZN[4]
- T g0 = zero.add(PTM[3] * PS[0]);
+ T g0 = zero.newInstance(PTM[3] * PS[0]);
if (alt.getReal() > ZN1[4]) {
g0 = g0.multiply(globe7(PS).multiply(sw[19]).add(1));
}
// Temperature at lower boundary
- T tlb = zero.add(PTM[1] * PD[3][0]);
+ T tlb = zero.newInstance(PTM[1] * PD[3][0]);
tlb = tlb.multiply(globe7(PD[3]).multiply(sw[17]).add(1));
// Slope
final T s = g0.divide(tinf.subtract(tlb));
// Lower thermosphere temp variations not significant for density above 300 km
- meso_tn1[1] = zero.add(PTM[6] * PTL[0][0]);
- meso_tn1[2] = zero.add(PTM[2] * PTL[1][0]);
- meso_tn1[3] = zero.add(PTM[7] * PTL[2][0]);
- meso_tn1[4] = zero.add(PTM[4] * PTL[3][0]);
- meso_tgn1[1] = zero.add(PTM[8] * PMA[8][0]);
+ meso_tn1[1] = zero.newInstance(PTM[6] * PTL[0][0]);
+ meso_tn1[2] = zero.newInstance(PTM[2] * PTL[1][0]);
+ meso_tn1[3] = zero.newInstance(PTM[7] * PTL[2][0]);
+ meso_tn1[4] = zero.newInstance(PTM[4] * PTL[3][0]);
+ meso_tgn1[1] = zero.newInstance(PTM[8] * PMA[8][0]);
if (alt.getReal() < 300.0) {
final double r = PTM[4] * PTL[3][0];
meso_tn1[1] = meso_tn1[1].divide(glob7s(PTL[0]).multiply(sw[18] ).negate().add(1));
@@ -2958,7 +2958,7 @@ void gts7(final T alt) {
}
/**** Temperature at altitude ****/
- setTemperature(ALTITUDE, densu(alt, zero.add(1.0), tinf, tlb, 0, 0, PTM[5], s));
+ setTemperature(ALTITUDE, densu(alt, zero.newInstance(1.0), tinf, tlb, 0, 0, PTM[5], s));
/**** N2 density ****/
/* Density variation factor at Zlb */
@@ -2999,7 +2999,7 @@ void gts7(final T alt) {
/* Turbopause */
final double zh04 = PDM[0][2];
/* Mixed density at Zlb */
- final T b04 = densu(zero.add(zh04), db04, tinf, tlb, HE_MASS - xmm, alpha[0] - 1., PTM[5], s);
+ final T b04 = densu(zero.newInstance(zh04), db04, tinf, tlb, HE_MASS - xmm, alpha[0] - 1., PTM[5], s);
/* Mixed density at Alt */
final T dm04 = densu(alt, b04, tinf, tlb, xmm, 0., PTM[5], s);
final double zhm04 = zhm28;
@@ -3025,7 +3025,7 @@ void gts7(final T alt) {
/* Turbopause */
final double zh16 = PDM[1][2];
/* Mixed density at Zlb */
- final T b16 = densu(zero.add(zh16), db16, tinf, tlb, O_MASS - xmm, alpha[1] - 1.0, PTM[5], s);
+ final T b16 = densu(zero.newInstance(zh16), db16, tinf, tlb, O_MASS - xmm, alpha[1] - 1.0, PTM[5], s);
/* Mixed density at Alt */
final T dm16 = densu(alt, b16, tinf, tlb, xmm, 0., PTM[5], s);
final double zhm16 = zhm28;
@@ -3041,7 +3041,7 @@ void gts7(final T alt) {
final double zcc16 = PDM[1][6] * PDL[1][12];
final double rc16 = PDM[1][3] * PDL[1][14];
/* Net density corrected at Alt */
- setDensity(ATOMIC_OXYGEN, diffusiveDensity.multiply(ccor(alt, zero.add(rc16), hcc16, zcc16)));
+ setDensity(ATOMIC_OXYGEN, diffusiveDensity.multiply(ccor(alt, zero.newInstance(rc16), hcc16, zcc16)));
}
/**** O2 density ****/
@@ -3057,7 +3057,7 @@ void gts7(final T alt) {
/* Turbopause */
final double zh32 = PDM[3][2];
/* Mixed density at Zlb */
- final T b32 = densu(zero.add(zh32), db32, tinf, tlb, O2_MASS - xmm, alpha[3] - 1., PTM[5], s);
+ final T b32 = densu(zero.newInstance(zh32), db32, tinf, tlb, O2_MASS - xmm, alpha[3] - 1., PTM[5], s);
/* Mixed density at Alt */
final T dm32 = densu(alt, b32, tinf, tlb, xmm, 0., PTM[5], s);
final double zhm32 = zhm28;
@@ -3090,7 +3090,7 @@ void gts7(final T alt) {
/* Turbopause */
final double zh40 = PDM[4][2];
/* Mixed density at Zlb */
- final T b40 = densu(zero.add(zh40), db40, tinf, tlb, AR_MASS - xmm, alpha[4] - 1., PTM[5], s);
+ final T b40 = densu(zero.newInstance(zh40), db40, tinf, tlb, AR_MASS - xmm, alpha[4] - 1., PTM[5], s);
/* Mixed density at Alt */
final T dm40 = densu(alt, b40, tinf, tlb, xmm, 0., PTM[5], s);
final double zhm40 = zhm28;
@@ -3116,7 +3116,7 @@ void gts7(final T alt) {
/* Turbopause */
final double zh01 = PDM[5][2];
/* Mixed density at Zlb */
- final T b01 = densu(zero.add(zh01), db01, tinf, tlb, H_MASS - xmm, alpha[6] - 1., PTM[5], s);
+ final T b01 = densu(zero.newInstance(zh01), db01, tinf, tlb, H_MASS - xmm, alpha[6] - 1., PTM[5], s);
/* Mixed density at Alt */
final T dm01 = densu(alt, b01, tinf, tlb, xmm, 0., PTM[5], s);
final double zhm01 = zhm28;
@@ -3132,7 +3132,7 @@ void gts7(final T alt) {
final double zcc01 = PDM[5][6] * PDL[1][18];
final double rc01 = PDM[5][3] * PDL[1][20];
/* Net density corrected at Alt */
- setDensity(HYDROGEN, diffusiveDensity.multiply(ccor(alt, zero.add(rc01), hcc01, zcc01)));
+ setDensity(HYDROGEN, diffusiveDensity.multiply(ccor(alt, zero.newInstance(rc01), hcc01, zcc01)));
}
/**** N density ****/
@@ -3147,7 +3147,7 @@ void gts7(final T alt) {
/* Turbopause */
final double zh14 = PDM[6][2];
/* Mixed density at Zlb */
- final T b14 = densu(zero.add(zh14), db14, tinf, tlb, N_MASS - xmm, alpha[7] - 1., PTM[5], s);
+ final T b14 = densu(zero.newInstance(zh14), db14, tinf, tlb, N_MASS - xmm, alpha[7] - 1., PTM[5], s);
/* Mixed density at Alt */
final T dm14 = densu(alt, b14, tinf, tlb, xmm, 0., PTM[5], s);
final double zhm14 = zhm28;
@@ -3163,14 +3163,14 @@ void gts7(final T alt) {
final double zcc14 = PDM[6][6] * PDL[0][3];
final double rc14 = PDM[6][3] * PDL[0][5];
/* Net density corrected at Alt */
- setDensity(ATOMIC_NITROGEN, diffusiveDensity.multiply(ccor(alt, zero.add(rc14), hcc14, zcc14)));
+ setDensity(ATOMIC_NITROGEN, diffusiveDensity.multiply(ccor(alt, zero.newInstance(rc14), hcc14, zcc14)));
}
/**** Anomalous O density ****/
final T g16h = globe7(PD[8]).multiply(sw[21]);
final T db16h = g16h.exp().multiply(PDM[7][0] * PD[8][0]);
final double tho = PDM[7][9] * PDL[0][6];
- diffusiveDensity = densu(alt, db16h, zero.add(tho), zero.add(tho), O_MASS, alpha[8], PTM[5], s);
+ diffusiveDensity = densu(alt, db16h, zero.newInstance(tho), zero.newInstance(tho), O_MASS, alpha[8], PTM[5], s);
final double zsht = PDM[7][5];
final double zmho = PDM[7][4];
final T zsho = scalh(zmho, O_MASS, tho);
@@ -3219,7 +3219,7 @@ void gts7(final T alt) {
void gtd7(final T alt) {
// Calculates for thermosphere/mesosphere (above ZN2[0])
- final T altt = (alt.getReal() > ZN2[0]) ? alt : zero.add(ZN2[0]);
+ final T altt = (alt.getReal() > ZN2[0]) ? alt : zero.newInstance(ZN2[0]);
gts7(altt);
if (alt.getReal() >= ZN2[0]) {
return;
@@ -3404,10 +3404,8 @@ private T globe7(final double[] p) {
// F10.7 effect
final double df = f107 - f107a;
final double dfa = f107a - FLUX_REF;
- t[0] = zero.add(p[19] * df * (1.0 + p[59] * dfa) +
- p[20] * df * df +
- p[21] * dfa +
- p[29] * dfa * dfa);
+ t[0] = zero.newInstance(p[19] * df * (1.0 + p[59] * dfa) + p[20] * df * df +
+ p[21] * dfa + p[29] * dfa * dfa);
final double f1 = 1.0 + (p[47] * dfa + p[19] * df + p[20] * df * df) * swc[1];
final double f2 = 1.0 + (p[49] * dfa + p[19] * df + p[20] * df * df) * swc[1];
@@ -3420,7 +3418,7 @@ private T globe7(final double[] p) {
add(plg[0][1].multiply(p[26]));
// Symmetrical annual
- t[2] = zero.add(p[18] * cd32);
+ t[2] = zero.newInstance(p[18] * cd32);
// Symmetrical semiannual
t[3] = plg[0][2].multiply(p[16]).add(p[15]).multiply(cd18);
@@ -3539,7 +3537,7 @@ private T globe7(final double[] p) {
}
// Sum all effects (params not used: 82, 89, 99, 139-149)
- T tinf = zero.add(p[30]);
+ T tinf = zero.newInstance(p[30]);
for (int i = 0; i < 14; i++) {
tinf = tinf.add(t[i].multiply(FastMath.abs(sw[i + 1])));
}
@@ -3562,7 +3560,7 @@ private T glob7s(final double[] p) {
final double cd39 = FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[38]));
// F10.7 effect
- t[0] = zero.add(p[21] * (f107a - FLUX_REF));
+ t[0] = zero.newInstance(p[21] * (f107a - FLUX_REF));
// Time independent
t[1] = plg[0][2].multiply(p[1]).
@@ -3715,7 +3713,7 @@ private T ccor2(final T alt, final double r, final double h1, final double zh, f
if (e1.getReal() > 70. || e2.getReal() > 70.) {
return field.getOne();
} else if (e1.getReal() < -70. && e2.getReal() < -70.) {
- return zero.add(FastMath.exp(r));
+ return zero.newInstance(FastMath.exp(r));
} else {
final T ex1 = e1.exp();
final T ex2 = e2.exp();
@@ -3849,7 +3847,7 @@ private T[] spline(final T[] x, final T[] y, final T yp1, final T ypn) {
final T[] u = MathArrays.buildArray(field, n);
if (yp1.getReal() < 1e+30) {
- y2[0] = zero.add(-0.5);
+ y2[0] = zero.newInstance(-0.5);
final T dx = x[1].subtract(x[0]);
final T dy = y[1].subtract(y[0]);
u[0] = dx.reciprocal().multiply(3.0).multiply(dy.divide(dx).subtract(yp1));
@@ -3896,20 +3894,20 @@ private T densm(final T alt, final T d0, final double xm) {
// stratosphere/mesosphere temperature
int mn = ZN2.length;
- T z = (alt.getReal() > ZN2[mn - 1]) ? alt : zero.add(ZN2[mn - 1]);
+ T z = (alt.getReal() > ZN2[mn - 1]) ? alt : zero.newInstance(ZN2[mn - 1]);
double z1 = ZN2[0];
double z2 = ZN2[mn - 1];
T t1 = meso_tn2[0];
T t2 = meso_tn2[mn - 1];
T zg = zeta(z, z1);
- T zgdif = zeta(zero.add(z2), z1);
+ T zgdif = zeta(zero.newInstance(z2), z1);
/* set up spline nodes */
T[] xs = MathArrays.buildArray(field, mn);
T[] ys = MathArrays.buildArray(field, mn);
for (int k = 0; k < mn; k++) {
- xs[k] = zeta(zero.add(ZN2[k]), z1).divide(zgdif);
+ xs[k] = zeta(zero.newInstance(ZN2[k]), z1).divide(zgdif);
ys[k] = meso_tn2[k].reciprocal();
}
final T qSM = rlat.add(z2).divide(rlat.add(z1));
@@ -3926,7 +3924,7 @@ private T densm(final T alt, final T d0, final double xm) {
if (xm != 0.0) {
/* calculate stratosphere / mesospehere density */
- final T glb = galt(zero.add(z1));
+ final T glb = galt(zero.newInstance(z1));
final T gamm = glb.multiply(zgdif).multiply(xm / R_GAS);
/* Integrate temperature profile */
@@ -3949,13 +3947,13 @@ private T densm(final T alt, final T d0, final double xm) {
t1 = meso_tn3[0];
t2 = meso_tn3[mn - 1];
zg = zeta(z, z1);
- zgdif = zeta(zero.add(z2), z1);
+ zgdif = zeta(zero.newInstance(z2), z1);
/* set up spline nodes */
xs = MathArrays.buildArray(field, mn);
ys = MathArrays.buildArray(field, mn);
for (int k = 0; k < mn; k++) {
- xs[k] = zeta(zero.add(ZN3[k]), z1).divide(zgdif);
+ xs[k] = zeta(zero.newInstance(ZN3[k]), z1).divide(zgdif);
ys[k] = meso_tn3[k].reciprocal();
}
final T qTS = rlat.add(z2) .divide(rlat.add(z1));
@@ -3972,7 +3970,7 @@ private T densm(final T alt, final T d0, final double xm) {
if (xm != 0.0) {
/* calculate tropospheric / stratosphere density */
- final T glb = galt(zero.add(z1));
+ final T glb = galt(zero.newInstance(z1));
final T gamm = glb.multiply(zgdif).multiply(xm / R_GAS);
/* Integrate temperature profile */
@@ -4001,7 +3999,7 @@ private T densu(final T alt, final T dlb, final T tinf,
final T tlb, final double xm, final double alpha,
final double zlb, final T s2) {
/* joining altitudes of Bates and spline */
- T z = (alt.getReal() > ZN1[0]) ? alt : zero.add(ZN1[0]);
+ T z = (alt.getReal() > ZN1[0]) ? alt : zero.newInstance(ZN1[0]);
/* geopotential altitude difference from ZLB */
final T zg2 = zeta(z, zlb);
@@ -4024,7 +4022,7 @@ private T densu(final T alt, final T dlb, final T tinf,
final T dta = tinf.subtract(ta).multiply(s2).multiply(p.multiply(p));
meso_tgn1[0] = dta;
meso_tn1[0] = ta;
- final T tzn1mn1 = zero.add(ZN1[mn - 1]);
+ final T tzn1mn1 = zero.newInstance(ZN1[mn - 1]);
z = (alt.getReal() > ZN1[mn - 1]) ? alt : tzn1mn1;
final T t1 = meso_tn1[0];
@@ -4034,7 +4032,7 @@ private T densu(final T alt, final T dlb, final T tinf,
zgdif = zeta(tzn1mn1, ZN1[0]);
/* set up spline nodes */
for (int k = 0; k < mn; k++) {
- xs[k] = zeta(zero.add(ZN1[k]), ZN1[0]).divide(zgdif);
+ xs[k] = zeta(zero.newInstance(ZN1[k]), ZN1[0]).divide(zgdif);
ys[k] = meso_tn1[k].reciprocal();
}
/* end node derivatives */
@@ -4054,20 +4052,20 @@ private T densu(final T alt, final T dlb, final T tinf,
}
/* calculate density above za */
- T glb = galt(zero.add(zlb));
+ T glb = galt(zero.newInstance(zlb));
T gamma = glb.divide(s2.multiply(tinf)).multiply(xm / R_GAS);
T expl = tt.getReal() <= 0 ?
- zero.add(50) :
+ zero.newInstance(50) :
min(50.0, s2.negate().multiply(gamma).multiply(zg2).exp());
T densu = dlb.multiply(tlb.divide(tt).pow(gamma.add(alpha + 1))).multiply(expl);
/* calculate density below za */
if (alt.getReal() < ZN1[0]) {
- glb = galt(zero.add(ZN1[0]));
+ glb = galt(zero.newInstance(ZN1[0]));
gamma = glb.multiply(zgdif).multiply(xm / R_GAS);
/* integrate spline temperatures */
expl = tz.getReal() <= 0 ?
- zero.add(50.0) :
+ zero.newInstance(50.0) :
min(50.0, gamma.multiply(splini(xs, ys, y2out, x)));
/* correct density at altitude */
densu = densu.multiply(meso_tn1[0].divide(tz).pow(alpha + 1).multiply(expl.negate().exp()));
@@ -4083,7 +4081,7 @@ private T densu(final T alt, final T dlb, final T tinf,
* @return min value
*/
private T min(final double d, final T f) {
- return (f.getReal() > d) ? zero.add(d) : f;
+ return (f.getReal() > d) ? zero.newInstance(d) : f;
}
/** Calculate gravity at altitude.
diff --git a/src/main/java/org/orekit/models/earth/ionosphere/FieldNeQuickParameters.java b/src/main/java/org/orekit/models/earth/ionosphere/FieldNeQuickParameters.java
index ee94cbb057..df37f73596 100644
--- a/src/main/java/org/orekit/models/earth/ionosphere/FieldNeQuickParameters.java
+++ b/src/main/java/org/orekit/models/earth/ionosphere/FieldNeQuickParameters.java
@@ -124,7 +124,7 @@ class FieldNeQuickParameters > {
}
// E layer maximum density height in km (Eq. 78)
- this.hmE = field.getZero().add(120.0);
+ this.hmE = field.getZero().newInstance(120.0);
// E layer critical frequency in MHz
final T foE = computefoE(date.getMonth(), az, xeff, latitude);
// E layer maximum density in 10^11 m-3 (Eq. 36)
@@ -162,8 +162,8 @@ class FieldNeQuickParameters > {
this.b2Bot = nmF2.divide(a).multiply(0.385);
this.b1Top = hmF2.subtract(hmF1).multiply(0.3);
this.b1Bot = hmF1.subtract(hmE).multiply(0.5);
- this.beTop = FastMath.max(b1Bot, zero.add(7.0));
- this.beBot = zero.add(5.0);
+ this.beTop = FastMath.max(b1Bot, zero.newInstance(7.0));
+ this.beBot = zero.newInstance(5.0);
// Layer amplitude coefficients
this.amplitudes = computeLayerAmplitudes(field, nmE, nmF1, foF1);
@@ -344,14 +344,14 @@ private T computeAz(final T modip, final double[] alpha) {
final T zero = field.getZero();
// Particular condition (Eq. 17)
if (alpha[0] == 0.0 && alpha[1] == 0.0 && alpha[2] == 0.0) {
- return zero.add(63.7);
+ return zero.newInstance(63.7);
}
// Az = a0 + modip * a1 + modip^2 * a2 (Eq. 18)
T az = modip.multiply(alpha[2]).add(alpha[1]).multiply(modip).add(alpha[0]);
// If Az < 0 -> Az = 0
az = FastMath.max(zero, az);
// If Az > 400 -> Az = 400
- az = FastMath.min(zero.add(400.0), az);
+ az = FastMath.min(zero.newInstance(400.0), az);
return az;
}
@@ -392,7 +392,8 @@ private T computeEffectiveSolarAngle(final int month,
final T angle = FastMath.atan2(FastMath.sqrt(cZenith.multiply(cZenith).negate().add(1.0)), cZenith);
final T x = FastMath.toDegrees(angle);
// Effective solar zenith angle (Eq. 28)
- final T xeff = join(clipExp(x.multiply(0.2).negate().add(20.0)).multiply(0.24).negate().add(90.0), x, zero.add(12.0), x.subtract(X0));
+ final T xeff = join(clipExp(x.multiply(0.2).negate().add(20.0)).multiply(0.24).negate().add(90.0), x,
+ zero.newInstance(12.0), x.subtract(X0));
return FastMath.toRadians(xeff);
}
@@ -441,7 +442,7 @@ private T computehmF2(final Field field, final T foE, final T foF2, final T m
final T zero = field.getZero();
// Ratio
final T fo = foF2.divide(foE);
- final T ratio = join(fo, zero.add(1.75), zero.add(20.0), fo.subtract(1.75));
+ final T ratio = join(fo, zero.newInstance(1.75), zero.newInstance(20.0), fo.subtract(1.75));
// deltaM parameter
T deltaM = zero.subtract(0.012);
@@ -624,9 +625,9 @@ private T computeMF2(final Field field, final T modip, final T[] cm3,
*/
private T computefoF1(final Field field, final T foE, final T foF2) {
final T zero = field.getZero();
- final T temp = join(foE.multiply(1.4), zero, zero.add(1000.0), foE.subtract(2.0));
- final T temp2 = join(zero, temp, zero.add(1000.0), foE.subtract(temp));
- final T value = join(temp2, temp2.multiply(0.85), zero.add(60.0), foF2.multiply(0.85).subtract(temp2));
+ final T temp = join(foE.multiply(1.4), zero, zero.newInstance(1000.0), foE.subtract(2.0));
+ final T temp2 = join(zero, temp, zero.newInstance(1000.0), foE.subtract(temp));
+ final T value = join(temp2, temp2.multiply(0.85), zero.newInstance(60.0), foF2.multiply(0.85).subtract(temp2));
if (value.getReal() < 1.0E-6) {
return zero;
} else {
@@ -674,7 +675,7 @@ private T[] computeLayerAmplitudes(final Field field, final T nmE, final T nm
a3a = nmE.subtract(epst(a2a, hmF1, b1Bot, hmE)).subtract(epst(a1, hmF2, b2Bot, hmE)).multiply(4.0);
}
amplitude[1] = a2a;
- amplitude[2] = join(a3a, zero.add(0.05), zero.add(60.0), a3a.subtract(0.005));
+ amplitude[2] = join(a3a, zero.newInstance(0.05), zero.newInstance(60.0), a3a.subtract(0.005));
}
return amplitude;
@@ -705,8 +706,8 @@ private T computeH0(final Field field, final int month, final T azr) {
}
// Auxiliary parameter kb (Eq. 101 and 102)
- T kb = join(ka, one.multiply(2.0), one, ka.subtract(2.0));
- kb = join(one.multiply(8.0), kb, one, kb.subtract(8.0));
+ T kb = join(ka, one.newInstance(2.0), one, ka.subtract(2.0));
+ kb = join(one.newInstance(8.0), kb, one, kb.subtract(8.0));
// Auxiliary parameter Ha (Eq. 103)
final T hA = kb.multiply(b2Bot);
@@ -732,9 +733,9 @@ private T computeH0(final Field field, final int month, final T azr) {
private T clipExp(final T power) {
final T zero = power.getField().getZero();
if (power.getReal() > 80.0) {
- return zero.add(5.5406E34);
+ return zero.newInstance(5.5406E34);
} else if (power.getReal() < -80) {
- return zero.add(1.8049E-35);
+ return zero.newInstance(1.8049E-35);
} else {
return FastMath.exp(power);
}
diff --git a/src/main/java/org/orekit/models/earth/ionosphere/KlobucharIonoModel.java b/src/main/java/org/orekit/models/earth/ionosphere/KlobucharIonoModel.java
index 1958f92e5f..16da1e053b 100644
--- a/src/main/java/org/orekit/models/earth/ionosphere/KlobucharIonoModel.java
+++ b/src/main/java/org/orekit/models/earth/ionosphere/KlobucharIonoModel.java
@@ -250,7 +250,7 @@ public > T pathDelay(final FieldAbsoluteDate> T pathDelay(final FieldAbsoluteDate> T bottomElectronDensity(final Field<
// Compute the exponential argument for each layer (Eq. 111 to 113)
// If h < 100km we use h = 100km as recommended in the reference document
- final T hTemp = FastMath.max(zero.add(100.0), h);
+ final T hTemp = FastMath.max(zero.newInstance(100.0), h);
final T exp = clipExp(field, FastMath.abs(hTemp.subtract(parameters.getHmF2())).add(1.0).divide(10.0).reciprocal());
final T[] arguments = MathArrays.buildArray(field, 3);
arguments[0] = hTemp.subtract(parameters.getHmF2()).divide(bf2);
@@ -655,9 +655,9 @@ private double clipExp(final double power) {
private > T clipExp(final Field field, final T power) {
final T zero = field.getZero();
if (power.getReal() > 80.0) {
- return zero.add(5.5406E34);
+ return zero.newInstance(5.5406E34);
} else if (power.getReal() < -80) {
- return zero.add(1.8049E-35);
+ return zero.newInstance(1.8049E-35);
} else {
return FastMath.exp(power);
}
diff --git a/src/main/java/org/orekit/models/earth/troposphere/FixedTroposphericDelay.java b/src/main/java/org/orekit/models/earth/troposphere/FixedTroposphericDelay.java
index 12e2aee5bb..5c7a5a34e3 100644
--- a/src/main/java/org/orekit/models/earth/troposphere/FixedTroposphericDelay.java
+++ b/src/main/java/org/orekit/models/earth/troposphere/FixedTroposphericDelay.java
@@ -147,7 +147,7 @@ public > T pathDelay(final T elevation, final
final T zero = date.getField().getZero();
final T pi = zero.getPi();
// limit the height to 5000 m
- final T h = FastMath.min(FastMath.max(zero, point.getAltitude()), zero.add(5000));
+ final T h = FastMath.min(FastMath.max(zero, point.getAltitude()), zero.newInstance(5000));
// limit the elevation to 0 - π
final T ele = FastMath.min(pi, FastMath.max(zero, elevation));
// mirror elevation at the right angle of π/2
diff --git a/src/main/java/org/orekit/models/earth/troposphere/GlobalMappingFunctionModel.java b/src/main/java/org/orekit/models/earth/troposphere/GlobalMappingFunctionModel.java
index 1c326c24df..be4bcf3d12 100644
--- a/src/main/java/org/orekit/models/earth/troposphere/GlobalMappingFunctionModel.java
+++ b/src/main/java/org/orekit/models/earth/troposphere/GlobalMappingFunctionModel.java
@@ -185,8 +185,8 @@ public > T[] mappingFactors(final T elevation,
final T zero = field.getZero();
// bh and ch constants (Boehm, J et al, 2006) | HYDROSTATIC PART
- final T bh = zero.add(0.0029);
- final T c0h = zero.add(0.062);
+ final T bh = zero.newInstance(0.0029);
+ final T c0h = zero.newInstance(0.062);
final T c10h;
final T c11h;
final T psi;
@@ -197,12 +197,12 @@ public > T[] mappingFactors(final T elevation,
// sin(latitude) > 0 -> northern hemisphere
if (FastMath.sin(latitude.getReal()) > 0) {
- c10h = zero.add(0.001);
- c11h = zero.add(0.005);
+ c10h = zero.newInstance(0.001);
+ c11h = zero.newInstance(0.005);
psi = zero;
} else {
- c10h = zero.add(0.002);
- c11h = zero.add(0.007);
+ c10h = zero.newInstance(0.002);
+ c11h = zero.newInstance(0.007);
psi = zero.getPi();
}
@@ -215,8 +215,8 @@ public > T[] mappingFactors(final T elevation,
final T ch = c11h.divide(2.0).multiply(FastMath.cos(coef).add(1.0)).add(c10h).multiply(FastMath.cos(latitude).negate().add(1.0)).add(c0h);
// bw and cw constants (Boehm, J et al, 2006) | WET PART
- final T bw = zero.add(0.00146);
- final T cw = zero.add(0.04391);
+ final T bw = zero.newInstance(0.00146);
+ final T cw = zero.newInstance(0.04391);
// Compute coefficients ah and aw with spherical harmonics Eq. 3 (Ref 1)
diff --git a/src/main/java/org/orekit/models/earth/troposphere/MendesPavlisModel.java b/src/main/java/org/orekit/models/earth/troposphere/MendesPavlisModel.java
index 4644e256b2..ad5d8625fa 100644
--- a/src/main/java/org/orekit/models/earth/troposphere/MendesPavlisModel.java
+++ b/src/main/java/org/orekit/models/earth/troposphere/MendesPavlisModel.java
@@ -201,7 +201,7 @@ public > T[] computeZenithDelay(final FieldGeo
final T[] delay = MathArrays.buildArray(field, 2);
// Dispertion Equation for the Hydrostatic component
- final T sigma = zero.add(1 / lambda);
+ final T sigma = zero.newInstance(1 / lambda);
final T sigma2 = sigma.multiply(sigma);
final T coef1 = sigma2.add(K_COEFFICIENTS[0]);
final T coef2 = sigma2.negate().add(K_COEFFICIENTS[0]);
diff --git a/src/main/java/org/orekit/models/earth/troposphere/NiellMappingFunctionModel.java b/src/main/java/org/orekit/models/earth/troposphere/NiellMappingFunctionModel.java
index a7392f71e1..0cbae2db97 100644
--- a/src/main/java/org/orekit/models/earth/troposphere/NiellMappingFunctionModel.java
+++ b/src/main/java/org/orekit/models/earth/troposphere/NiellMappingFunctionModel.java
@@ -236,8 +236,8 @@ public > T[] mappingFactors(final T elevation,
function[0] = TroposphericModelUtils.mappingFunction(ah, bh, ch, elevation);
// Wet mapping factor
- function[1] = TroposphericModelUtils.mappingFunction(zero.add(awFunction.value(absLatidude)), zero.add(bwFunction.value(absLatidude)),
- zero.add(cwFunction.value(absLatidude)), elevation);
+ function[1] = TroposphericModelUtils.mappingFunction(zero.newInstance(awFunction.value(absLatidude)), zero.newInstance(bwFunction.value(absLatidude)),
+ zero.newInstance(cwFunction.value(absLatidude)), elevation);
// Apply height correction
final T correction = TroposphericModelUtils.computeHeightCorrection(elevation, point.getAltitude(), field);
diff --git a/src/main/java/org/orekit/models/earth/troposphere/SaastamoinenModel.java b/src/main/java/org/orekit/models/earth/troposphere/SaastamoinenModel.java
index 78c9c99883..ab728dca0a 100644
--- a/src/main/java/org/orekit/models/earth/troposphere/SaastamoinenModel.java
+++ b/src/main/java/org/orekit/models/earth/troposphere/SaastamoinenModel.java
@@ -293,7 +293,7 @@ public > T pathDelay(final T elevation, final
final T zero = field.getZero();
// there are no data in the model for negative altitudes and altitude bigger than 5000 m
// limit the height to a range of [0, 5000] m
- final T fixedHeight = FastMath.min(FastMath.max(zero, point.getAltitude()), zero.add(5000));
+ final T fixedHeight = FastMath.min(FastMath.max(zero, point.getAltitude()), zero.newInstance(5000));
// the corrected temperature using a temperature gradient of -6.5 K/km
final T T = fixedHeight.multiply(6.5e-3).negate().add(t0);
@@ -308,7 +308,7 @@ public > T pathDelay(final T elevation, final
final T e = R.multiply(FastMath.exp(eFunction.value(T)));
// calculate the zenith angle from the elevation
- final T z = FastMath.abs(FastMath.max(elevation, zero.add(lowElevationThreshold)).negate().add(zero.getPi().multiply(0.5)));
+ final T z = FastMath.abs(FastMath.max(elevation, zero.newInstance(lowElevationThreshold)).negate().add(zero.getPi().multiply(0.5)));
// get correction factor
final T deltaR = getDeltaR(fixedHeight, z, field);
@@ -346,7 +346,7 @@ private