GlobalMappingFunctionModel.java
/* Copyright 2002-2023 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.models.earth.troposphere;
import org.hipparchus.Field;
import org.hipparchus.CalculusFieldElement;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.FieldSinCos;
import org.hipparchus.util.MathArrays;
import org.hipparchus.util.SinCos;
import org.orekit.annotation.DefaultDataContext;
import org.orekit.bodies.FieldGeodeticPoint;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.data.DataContext;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.DateTimeComponents;
import org.orekit.time.FieldAbsoluteDate;
import org.orekit.time.TimeScale;
import org.orekit.utils.FieldLegendrePolynomials;
import org.orekit.utils.LegendrePolynomials;
/** The Global Mapping Function model for radio techniques.
* This model is an empirical mapping function. It only needs the
* values of the station latitude, longitude, height and the
* date for the computations.
* <p>
* The Global Mapping Function is based on spherical harmonics up
* to degree and order of 9. It was developed to be consistent
* with the {@link ViennaOneModel Vienna1} mapping function model.
* </p>
*
* @see "Boehm, J., A.E. Niell, P. Tregoning, H. Schuh (2006),
* Global Mapping Functions (GMF): A new empirical mapping function based
* on numerical weather model data, Geoph. Res. Letters, Vol. 33, L07304,
* doi:10.1029/2005GL025545."
*
* @see "Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
* IERS Technical Note No. 36, BKG (2010)"
*
* @author Bryan Cazabonne
*
*/
public class GlobalMappingFunctionModel implements MappingFunction {
/** Multiplication factor for mapping function coefficients. */
private static final double FACTOR = 1.0e-5;
/** UTC time scale. */
private final TimeScale utc;
/** Build a new instance.
*
* <p>This constructor uses the {@link DataContext#getDefault() default data context}.
*
* @see #GlobalMappingFunctionModel(TimeScale)
*/
@DefaultDataContext
public GlobalMappingFunctionModel() {
this(DataContext.getDefault().getTimeScales().getUTC());
}
/** Build a new instance.
* @param utc UTC time scale.
* @since 10.1
*/
public GlobalMappingFunctionModel(final TimeScale utc) {
this.utc = utc;
}
/** {@inheritDoc} */
@Override
public double[] mappingFactors(final double elevation, final GeodeticPoint point,
final AbsoluteDate date) {
// Day of year computation
final DateTimeComponents dtc = date.getComponents(utc);
final int dofyear = dtc.getDate().getDayOfYear();
// bh and ch constants (Boehm, J et al, 2006) | HYDROSTATIC PART
final double bh = 0.0029;
final double c0h = 0.062;
final double c10h;
final double c11h;
final double psi;
// Latitude and longitude of the station
final double latitude = point.getLatitude();
final double longitude = point.getLongitude();
if (FastMath.sin(latitude) > 0) {
// northern hemisphere case
c10h = 0.001;
c11h = 0.005;
psi = 0.0;
} else {
// southern hemisphere case
c10h = 0.002;
c11h = 0.007;
psi = FastMath.PI;
}
double t0 = 28;
if (latitude < 0) {
// southern hemisphere: t0 = 28 + an integer half of year
t0 += 183;
}
final double coef = ((dofyear + 1 - t0) / 365.25) * 2 * FastMath.PI + psi;
final double ch = c0h + ((FastMath.cos(coef) + 1) * (c11h / 2.0) + c10h) * (1.0 - FastMath.cos(latitude));
// bw and cw constants (Boehm, J et al, 2006) | WET PART
final double bw = 0.00146;
final double cw = 0.04391;
// Compute coefficients ah and aw with spherical harmonics Eq. 3 (Ref 1)
// Compute Legendre Polynomials Pnm(sin(phi))
final int degree = 9;
final int order = 9;
final LegendrePolynomials p = new LegendrePolynomials(degree, order, FastMath.sin(latitude));
double a0Hydro = 0.;
double amplHydro = 0.;
double a0Wet = 0.;
double amplWet = 0.;
final ABCoefficients abCoef = new ABCoefficients();
int j = 0;
for (int n = 0; n <= 9; n++) {
for (int m = 0; m <= n; m++) {
// Sine and cosine of m * longitude
final SinCos sc = FastMath.sinCos(m * longitude);
// Compute coefficients
a0Hydro = a0Hydro + (abCoef.getAHMean(j) * p.getPnm(n, m) * sc.cos() +
abCoef.getBHMean(j) * p.getPnm(n, m) * sc.sin()) * FACTOR;
a0Wet = a0Wet + (abCoef.getAWMean(j) * p.getPnm(n, m) * sc.cos() +
abCoef.getBWMean(j) * p.getPnm(n, m) * sc.sin()) * FACTOR;
amplHydro = amplHydro + (abCoef.getAHAmplitude(j) * p.getPnm(n, m) * sc.cos() +
abCoef.getBHAmplitude(j) * p.getPnm(n, m) * sc.sin()) * FACTOR;
amplWet = amplWet + (abCoef.getAWAmplitude(j) * p.getPnm(n, m) * sc.cos() +
abCoef.getBWAmplitude(j) * p.getPnm(n, m) * sc.sin()) * FACTOR;
j = j + 1;
}
}
// Eq. 2 (Ref 1)
final double ah = a0Hydro + amplHydro * FastMath.cos(coef - psi);
final double aw = a0Wet + amplWet * FastMath.cos(coef - psi);
final double[] function = new double[2];
function[0] = TroposphericModelUtils.mappingFunction(ah, bh, ch, elevation);
function[1] = TroposphericModelUtils.mappingFunction(aw, bw, cw, elevation);
// Apply height correction
final double correction = TroposphericModelUtils.computeHeightCorrection(elevation, point.getAltitude());
function[0] = function[0] + correction;
return function;
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T[] mappingFactors(final T elevation, final FieldGeodeticPoint<T> point,
final FieldAbsoluteDate<T> date) {
// Day of year computation
final DateTimeComponents dtc = date.getComponents(utc);
final int dofyear = dtc.getDate().getDayOfYear();
final Field<T> field = date.getField();
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 c10h;
final T c11h;
final T psi;
// Latitude and longitude of the station
final T latitude = point.getLatitude();
final T longitude = point.getLongitude();
// sin(latitude) > 0 -> northern hemisphere
if (FastMath.sin(latitude.getReal()) > 0) {
c10h = zero.add(0.001);
c11h = zero.add(0.005);
psi = zero;
} else {
c10h = zero.add(0.002);
c11h = zero.add(0.007);
psi = zero.getPi();
}
double t0 = 28;
if (latitude.getReal() < 0) {
// southern hemisphere: t0 = 28 + an integer half of year
t0 += 183;
}
final T coef = psi.add(zero.getPi().multiply(2.0).multiply((dofyear + 1 - t0) / 365.25));
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);
// Compute coefficients ah and aw with spherical harmonics Eq. 3 (Ref 1)
// Compute Legendre Polynomials Pnm(sin(phi))
final int degree = 9;
final int order = 9;
final FieldLegendrePolynomials<T> p = new FieldLegendrePolynomials<>(degree, order, FastMath.sin(latitude));
T a0Hydro = zero;
T amplHydro = zero;
T a0Wet = zero;
T amplWet = zero;
final ABCoefficients abCoef = new ABCoefficients();
int j = 0;
for (int n = 0; n <= 9; n++) {
for (int m = 0; m <= n; m++) {
// Sine and cosine of m * longitude
final FieldSinCos<T> sc = FastMath.sinCos(longitude.multiply(m));
// Compute coefficients
a0Hydro = a0Hydro.add(p.getPnm(n, m).multiply(abCoef.getAHMean(j)).multiply(sc.cos()).
add(p.getPnm(n, m).multiply(abCoef.getBHMean(j)).multiply(sc.sin())).
multiply(FACTOR));
a0Wet = a0Wet.add(p.getPnm(n, m).multiply(abCoef.getAWMean(j)).multiply(sc.cos()).
add(p.getPnm(n, m).multiply(abCoef.getBWMean(j)).multiply(sc.sin())).
multiply(FACTOR));
amplHydro = amplHydro.add(p.getPnm(n, m).multiply(abCoef.getAHAmplitude(j)).multiply(sc.cos()).
add(p.getPnm(n, m).multiply(abCoef.getBHAmplitude(j)).multiply(sc.sin())).
multiply(FACTOR));
amplWet = amplWet.add(p.getPnm(n, m).multiply(abCoef.getAWAmplitude(j)).multiply(sc.cos()).
add(p.getPnm(n, m).multiply(abCoef.getBWAmplitude(j)).multiply(sc.sin())).
multiply(FACTOR));
j = j + 1;
}
}
// Eq. 2 (Ref 1)
final T ah = a0Hydro.add(amplHydro.multiply(FastMath.cos(coef.subtract(psi))));
final T aw = a0Wet.add(amplWet.multiply(FastMath.cos(coef.subtract(psi))));
final T[] function = MathArrays.buildArray(field, 2);
function[0] = TroposphericModelUtils.mappingFunction(ah, bh, ch, elevation);
function[1] = TroposphericModelUtils.mappingFunction(aw, bw, cw, elevation);
// Apply height correction
final T correction = TroposphericModelUtils.computeHeightCorrection(elevation, point.getAltitude(), field);
function[0] = function[0].add(correction);
return function;
}
private static class ABCoefficients {
/** Mean hydrostatic coefficients a.*/
private static final double[] AH_MEAN = {
1.2517e02,
8.503e-01,
6.936e-02,
-6.760e+00,
1.771e-01,
1.130e-02,
5.963e-01,
1.808e-02,
2.801e-03,
-1.414e-03,
-1.212e+00,
9.300e-02,
3.683e-03,
1.095e-03,
4.671e-05,
3.959e-01,
-3.867e-02,
5.413e-03,
-5.289e-04,
3.229e-04,
2.067e-05,
3.000e-01,
2.031e-02,
5.900e-03,
4.573e-04,
-7.619e-05,
2.327e-06,
3.845e-06,
1.182e-01,
1.158e-02,
5.445e-03,
6.219e-05,
4.204e-06,
-2.093e-06,
1.540e-07,
-4.280e-08,
-4.751e-01,
-3.490e-02,
1.758e-03,
4.019e-04,
-2.799e-06,
-1.287e-06,
5.468e-07,
7.580e-08,
-6.300e-09,
-1.160e-01,
8.301e-03,
8.771e-04,
9.955e-05,
-1.718e-06,
-2.012e-06,
1.170e-08,
1.790e-08,
-1.300e-09,
1.000e-10
};
/** Mean hydrostatic coefficients b.*/
private static final double[] BH_MEAN = {
0.000e+00,
0.000e+00,
3.249e-02,
0.000e+00,
3.324e-02,
1.850e-02,
0.000e+00,
-1.115e-01,
2.519e-02,
4.923e-03,
0.000e+00,
2.737e-02,
1.595e-02,
-7.332e-04,
1.933e-04,
0.000e+00,
-4.796e-02,
6.381e-03,
-1.599e-04,
-3.685e-04,
1.815e-05,
0.000e+00,
7.033e-02,
2.426e-03,
-1.111e-03,
-1.357e-04,
-7.828e-06,
2.547e-06,
0.000e+00,
5.779e-03,
3.133e-03,
-5.312e-04,
-2.028e-05,
2.323e-07,
-9.100e-08,
-1.650e-08,
0.000e+00,
3.688e-02,
-8.638e-04,
-8.514e-05,
-2.828e-05,
5.403e-07,
4.390e-07,
1.350e-08,
1.800e-09,
0.000e+00,
-2.736e-02,
-2.977e-04,
8.113e-05,
2.329e-07,
8.451e-07,
4.490e-08,
-8.100e-09,
-1.500e-09,
2.000e-10
};
/** Amplitude for hydrostatic coefficients a.*/
private static final double[] AH_AMPL = {
-2.738e-01,
-2.837e+00,
1.298e-02,
-3.588e-01,
2.413e-02,
3.427e-02,
-7.624e-01,
7.272e-02,
2.160e-02,
-3.385e-03,
4.424e-01,
3.722e-02,
2.195e-02,
-1.503e-03,
2.426e-04,
3.013e-01,
5.762e-02,
1.019e-02,
-4.476e-04,
6.790e-05,
3.227e-05,
3.123e-01,
-3.535e-02,
4.840e-03,
3.025e-06,
-4.363e-05,
2.854e-07,
-1.286e-06,
-6.725e-01,
-3.730e-02,
8.964e-04,
1.399e-04,
-3.990e-06,
7.431e-06,
-2.796e-07,
-1.601e-07,
4.068e-02,
-1.352e-02,
7.282e-04,
9.594e-05,
2.070e-06,
-9.620e-08,
-2.742e-07,
-6.370e-08,
-6.300e-09,
8.625e-02,
-5.971e-03,
4.705e-04,
2.335e-05,
4.226e-06,
2.475e-07,
-8.850e-08,
-3.600e-08,
-2.900e-09,
0.000e+00
};
/** Amplitude for hydrostatic coefficients b.*/
private static final double[] BH_AMPL = {
0.000e+00,
0.000e+00,
-1.136e-01,
0.000e+00,
-1.868e-01,
-1.399e-02,
0.000e+00,
-1.043e-01,
1.175e-02,
-2.240e-03,
0.000e+00,
-3.222e-02,
1.333e-02,
-2.647e-03,
-2.316e-05,
0.000e+00,
5.339e-02,
1.107e-02,
-3.116e-03,
-1.079e-04,
-1.299e-05,
0.000e+00,
4.861e-03,
8.891e-03,
-6.448e-04,
-1.279e-05,
6.358e-06,
-1.417e-07,
0.000e+00,
3.041e-02,
1.150e-03,
-8.743e-04,
-2.781e-05,
6.367e-07,
-1.140e-08,
-4.200e-08,
0.000e+00,
-2.982e-02,
-3.000e-03,
1.394e-05,
-3.290e-05,
-1.705e-07,
7.440e-08,
2.720e-08,
-6.600e-09,
0.000e+00,
1.236e-02,
-9.981e-04,
-3.792e-05,
-1.355e-05,
1.162e-06,
-1.789e-07,
1.470e-08,
-2.400e-09,
-4.000e-10
};
/** Mean wet coefficients a.*/
private static final double[] AW_MEAN = {
5.640e+01,
1.555e+00,
-1.011e+00,
-3.975e+00,
3.171e-02,
1.065e-01,
6.175e-01,
1.376e-01,
4.229e-02,
3.028e-03,
1.688e+00,
-1.692e-01,
5.478e-02,
2.473e-02,
6.059e-04,
2.278e+00,
6.614e-03,
-3.505e-04,
-6.697e-03,
8.402e-04,
7.033e-04,
-3.236e+00,
2.184e-01,
-4.611e-02,
-1.613e-02,
-1.604e-03,
5.420e-05,
7.922e-05,
-2.711e-01,
-4.406e-01,
-3.376e-02,
-2.801e-03,
-4.090e-04,
-2.056e-05,
6.894e-06,
2.317e-06,
1.941e+00,
-2.562e-01,
1.598e-02,
5.449e-03,
3.544e-04,
1.148e-05,
7.503e-06,
-5.667e-07,
-3.660e-08,
8.683e-01,
-5.931e-02,
-1.864e-03,
-1.277e-04,
2.029e-04,
1.269e-05,
1.629e-06,
9.660e-08,
-1.015e-07,
-5.000e-10
};
/** Mean wet coefficients b.*/
private static final double[] BW_MEAN = {
0.000e+00,
0.000e+00,
2.592e-01,
0.000e+00,
2.974e-02,
-5.471e-01,
0.000e+00,
-5.926e-01,
-1.030e-01,
-1.567e-02,
0.000e+00,
1.710e-01,
9.025e-02,
2.689e-02,
2.243e-03,
0.000e+00,
3.439e-01,
2.402e-02,
5.410e-03,
1.601e-03,
9.669e-05,
0.000e+00,
9.502e-02,
-3.063e-02,
-1.055e-03,
-1.067e-04,
-1.130e-04,
2.124e-05,
0.000e+00,
-3.129e-01,
8.463e-03,
2.253e-04,
7.413e-05,
-9.376e-05,
-1.606e-06,
2.060e-06,
0.000e+00,
2.739e-01,
1.167e-03,
-2.246e-05,
-1.287e-04,
-2.438e-05,
-7.561e-07,
1.158e-06,
4.950e-08,
0.000e+00,
-1.344e-01,
5.342e-03,
3.775e-04,
-6.756e-05,
-1.686e-06,
-1.184e-06,
2.768e-07,
2.730e-08,
5.700e-09
};
/** Amplitude for wet coefficients a.*/
private static final double[] AW_AMPL = {
1.023e-01,
-2.695e+00,
3.417e-01,
-1.405e-01,
3.175e-01,
2.116e-01,
3.536e+00,
-1.505e-01,
-1.660e-02,
2.967e-02,
3.819e-01,
-1.695e-01,
-7.444e-02,
7.409e-03,
-6.262e-03,
-1.836e+00,
-1.759e-02,
-6.256e-02,
-2.371e-03,
7.947e-04,
1.501e-04,
-8.603e-01,
-1.360e-01,
-3.629e-02,
-3.706e-03,
-2.976e-04,
1.857e-05,
3.021e-05,
2.248e+00,
-1.178e-01,
1.255e-02,
1.134e-03,
-2.161e-04,
-5.817e-06,
8.836e-07,
-1.769e-07,
7.313e-01,
-1.188e-01,
1.145e-02,
1.011e-03,
1.083e-04,
2.570e-06,
-2.140e-06,
-5.710e-08,
2.000e-08,
-1.632e+00,
-6.948e-03,
-3.893e-03,
8.592e-04,
7.577e-05,
4.539e-06,
-3.852e-07,
-2.213e-07,
-1.370e-08,
5.800e-09
};
/** Amplitude for wet coefficients b.*/
private static final double[] BW_AMPL = {
0.000e+00,
0.000e+00,
-8.865e-02,
0.000e+00,
-4.309e-01,
6.340e-02,
0.000e+00,
1.162e-01,
6.176e-02,
-4.234e-03,
0.000e+00,
2.530e-01,
4.017e-02,
-6.204e-03,
4.977e-03,
0.000e+00,
-1.737e-01,
-5.638e-03,
1.488e-04,
4.857e-04,
-1.809e-04,
0.000e+00,
-1.514e-01,
-1.685e-02,
5.333e-03,
-7.611e-05,
2.394e-05,
8.195e-06,
0.000e+00,
9.326e-02,
-1.275e-02,
-3.071e-04,
5.374e-05,
-3.391e-05,
-7.436e-06,
6.747e-07,
0.000e+00,
-8.637e-02,
-3.807e-03,
-6.833e-04,
-3.861e-05,
-2.268e-05,
1.454e-06,
3.860e-07,
-1.068e-07,
0.000e+00,
-2.658e-02,
-1.947e-03,
7.131e-04,
-3.506e-05,
1.885e-07,
5.792e-07,
3.990e-08,
2.000e-08,
-5.700e-09
};
/** Build a new instance. */
ABCoefficients() {
}
/** Get the value of the mean hydrostatique coefficient ah for the given index.
* @param index index
* @return the mean hydrostatique coefficient ah for the given index
*/
public double getAHMean(final int index) {
return AH_MEAN[index];
}
/** Get the value of the mean hydrostatique coefficient bh for the given index.
* @param index index
* @return the mean hydrostatique coefficient bh for the given index
*/
public double getBHMean(final int index) {
return BH_MEAN[index];
}
/** Get the value of the mean wet coefficient aw for the given index.
* @param index index
* @return the mean wet coefficient aw for the given index
*/
public double getAWMean(final int index) {
return AW_MEAN[index];
}
/** Get the value of the mean wet coefficient bw for the given index.
* @param index index
* @return the mean wet coefficient bw for the given index
*/
public double getBWMean(final int index) {
return BW_MEAN[index];
}
/** Get the value of the amplitude of the hydrostatique coefficient ah for the given index.
* @param index index
* @return the amplitude of the hydrostatique coefficient ah for the given index
*/
public double getAHAmplitude(final int index) {
return AH_AMPL[index];
}
/** Get the value of the amplitude of the hydrostatique coefficient bh for the given index.
* @param index index
* @return the amplitude of the hydrostatique coefficient bh for the given index
*/
public double getBHAmplitude(final int index) {
return BH_AMPL[index];
}
/** Get the value of the amplitude of the wet coefficient aw for the given index.
* @param index index
* @return the amplitude of the wet coefficient aw for the given index
*/
public double getAWAmplitude(final int index) {
return AW_AMPL[index];
}
/** Get the value of the amplitude of the wet coefficient bw for the given index.
* @param index index
* @return the amplitude of the wet coefficient bw for the given index
*/
public double getBWAmplitude(final int index) {
return BW_AMPL[index];
}
}
}