ModifiedHopfieldModel.java
/* Copyright 2002-2024 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 org.hipparchus.CalculusFieldElement;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.FieldSinCos;
import org.hipparchus.util.MathUtils;
import org.hipparchus.util.SinCos;
import org.orekit.bodies.FieldGeodeticPoint;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.models.earth.weather.FieldPressureTemperatureHumidity;
import org.orekit.models.earth.weather.PressureTemperatureHumidity;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.FieldAbsoluteDate;
import org.orekit.utils.FieldTrackingCoordinates;
import org.orekit.utils.ParameterDriver;
import org.orekit.utils.TrackingCoordinates;
/** The modified Hopfield model.
* <p>
* This model from Hopfield 1969, 1970, 1972 is described in equations
* 5.105, 5.106, 5.107 and 5.108 in Guochang Xu, GPS - Theory, Algorithms
* and Applications, Springer, 2007.
* </p>
* @author Luc Maisonobe
* @see "Guochang Xu, GPS - Theory, Algorithms and Applications, Springer, 2007"
* @since 12.1
*/
public class ModifiedHopfieldModel implements TroposphericModel {
/** Constant for dry altitude effect. */
private static final double HD0 = 40136.0;
/** Slope for dry altitude effect. */
private static final double HD1 = 148.72;
/** Temperature reference. */
private static final double T0 = 273.16;
/** Constant for wet altitude effect. */
private static final double HW0 = 11000.0;
/** Dry delay factor. */
private static final double ND = 77.64e-6;
/** Wet delay factor, degree 1. */
private static final double NW1 = -12.96e-6;
/** Wet delay factor, degree 2. */
private static final double NW2 = 0.371800;
/** BAse radius. */
private static final double RE = 6378137.0;
/** Create a new Hopfield model.
*/
public ModifiedHopfieldModel() {
// nothing to do
}
/** {@inheritDoc} */
@Override
public TroposphericDelay pathDelay(final TrackingCoordinates trackingCoordinates,
final GeodeticPoint point,
final PressureTemperatureHumidity weather,
final double[] parameters, final AbsoluteDate date) {
// zenith angle
final double zenithAngle = MathUtils.SEMI_PI - trackingCoordinates.getElevation();
// dry component
final double hd = HD0 + HD1 * (weather.getTemperature() - T0);
final double nd = ND * TroposphericModelUtils.HECTO_PASCAL.fromSI(weather.getPressure()) /
weather.getTemperature();
// wet component
final double hw = HW0;
final double nw = (NW1 + NW2 / weather.getTemperature()) / weather.getTemperature();
return new TroposphericDelay(delay(0.0, hd, nd),
delay(0.0, hw, nw),
delay(zenithAngle, hd, nd),
delay(zenithAngle, hw, nw));
}
/** {@inheritDoc}
* <p>
* 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.
* </p>
* <p>
* 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.
* </p>
*/
@Override
public <T extends CalculusFieldElement<T>> FieldTroposphericDelay<T> pathDelay(final FieldTrackingCoordinates<T> trackingCoordinates,
final FieldGeodeticPoint<T> point,
final FieldPressureTemperatureHumidity<T> weather,
final T[] parameters, final FieldAbsoluteDate<T> date) {
// zenith angle
final T zenithAngle = trackingCoordinates.getElevation().negate().add(MathUtils.SEMI_PI);
// dry component
final T hd = weather.getTemperature().subtract(T0).multiply(HD1).add(HD0);
final T nd = TroposphericModelUtils.HECTO_PASCAL.fromSI(weather.getPressure()).
multiply(ND).
divide(weather.getTemperature());
// wet component
final T hw = date.getField().getZero().newInstance(HW0);
final T nw = weather.getTemperature().reciprocal().multiply(NW2).add(NW1).divide(weather.getTemperature());
return new FieldTroposphericDelay<>(delay(date.getField().getZero(), hd, nd),
delay(date.getField().getZero(), hw, nw),
delay(zenithAngle, hd, nd),
delay(zenithAngle, hw, nw));
}
/** {@inheritDoc} */
@Override
public List<ParameterDriver> getParametersDrivers() {
return Collections.emptyList();
}
/** Compute the 9 terms sum delay.
* @param zenithAngle zenith angle
* @param hi altitude effect
* @param ni delay factor
* @return 9 terms sum delay
*/
private double delay(final double zenithAngle, final double hi, final double ni) {
// equation 5.107
final SinCos scZ = FastMath.sinCos(zenithAngle);
final double rePhi = RE + hi;
final double reS = RE * scZ.sin();
final double reC = RE * scZ.cos();
final double ri = FastMath.sqrt(rePhi * rePhi - reS * reS) - reC;
final double ai = -scZ.cos() / hi;
final double bi = -scZ.sin() * scZ.sin() / (2 * hi * RE);
final double ai2 = ai * ai;
final double bi2 = bi * bi;
final double f1i = 1;
final double f2i = 4 * ai;
final double f3i = 6 * ai2 + 4 * bi;
final double f4i = 4 * ai * (ai2 + 3 * bi);
final double f5i = ai2 * ai2 + 12 * ai2 * bi + 6 * bi2;
final double f6i = 4 * ai * bi * (ai2 + 3 * bi);
final double f7i = bi2 * (6 * ai2 + 4 * bi);
final double f8i = 4 * ai * bi * bi2;
final double f9i = bi2 * bi2;
return ni * (ri * (f1i +
ri * (f2i / 2 +
ri * (f3i / 3 +
ri * (f4i / 4 +
ri * (f5i / 5 +
ri * (f6i / 6 +
ri * (f7i / 7 +
ri * (f8i / 8 +
ri * f9i / 9)))))))));
}
/** Compute the 9 terms sum delay.
* @param <T> type of the elements
* @param zenithAngle zenith angle
* @param hi altitude effect
* @param ni delay factor
* @return 9 terms sum delay
*/
private <T extends CalculusFieldElement<T>> T delay(final T zenithAngle, final T hi, final T ni) {
// equation 5.107
final FieldSinCos<T> scZ = FastMath.sinCos(zenithAngle);
final T rePhi = hi.add(RE);
final T reS = scZ.sin().multiply(RE);
final T reC = scZ.cos().multiply(RE);
final T ri = FastMath.sqrt(rePhi.multiply(rePhi).subtract(reS.multiply(reS))).subtract(reC);
final T ai = scZ.cos().negate().divide(hi);
final T bi = scZ.sin().multiply(scZ.sin()).negate().divide(hi.add(hi).multiply(RE));
final T ai2 = ai.multiply(ai);
final T bi2 = bi.multiply(bi);
final T f1i = ai.getField().getOne();
final T f2i = ai.multiply(4);
final T f3i = ai2.multiply(6).add(bi.multiply(4));
final T f4i = ai.multiply(4).multiply(ai2.add(bi.multiply(3)));
final T f5i = ai2.multiply(ai2).add(ai2.multiply(12).multiply(bi)).add(bi2.multiply(6));
final T f6i = ai.multiply(4).multiply(bi).multiply(ai2.add(bi.multiply(3)));
final T f7i = bi2.multiply(ai2.multiply(6).add(bi.multiply(4)));
final T f8i = ai.multiply(4).multiply(bi).multiply(bi2);
final T f9i = bi2.multiply(bi2);
return ni.
multiply(ri.
multiply(f1i.
add(ri.
multiply(f2i.divide(2).
add(ri.
multiply(f3i.divide(3).
add(ri.
multiply(f4i.divide(4).
add(ri.
multiply(f5i.divide(5).
add(ri.
multiply(f6i.divide(6).
add(ri.
multiply(f7i.divide(7).
add(ri.
multiply(f8i.divide(8).
add(ri.multiply(f9i.divide(9)))))))))))))))))));
}
}