GlobalPressureTemperature.java

/* Copyright 2002-2025 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.weather;

import org.hipparchus.util.FastMath;
import org.hipparchus.util.SinCos;
import org.orekit.annotation.DefaultDataContext;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.data.DataContext;
import org.orekit.models.earth.Geoid;
import org.orekit.models.earth.troposphere.TroposphericModelUtils;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.DateTimeComponents;
import org.orekit.time.TimeScale;
import org.orekit.utils.LegendrePolynomials;

/** The Global Pressure and Temperature model.
 * This model is an empirical model that provides the temperature and the pressure depending
 * the latitude and the longitude of the station.
 * <p>
 * The Global Pressure and Temperature model is based on spherical harmonics up
 * to degree and order of 9. The residual values ​​of this model can reach 20 hPa
 * for pressure and 10 ° C for temperature. They are significant for higher latitudes and
 * small near the equator (Böhm, 2007)
 * </p>
 *
 * @see "J. Böhm, R. Heinkelmann, and H. Schuh (2007),
 * Short Note: A global model of pressure and temperature for geodetic applications. J Geod,
 * doi:10.1007/s00190-007-0135-3."
 *
 * @author Bryan Cazabonne
 * @author Luc Maisonobe
 * @since 12.1
 */
public class GlobalPressureTemperature {

    /** Temperature gradient (°C/m). */
    private static final double TEMPERATURE_GRADIENT = -6.5e-3;

    /** Spherical harmonics degree. */
    private static final int DEGREE = 9;

    /** Spherical harmonics order. */
    private static final int ORDER = 9;

    /** Geoid used to compute the undulations. */
    private final Geoid geoid;

    /** UTC time scale. */
    private final TimeScale utc;

    /** Build a new instance.
     *
     * <p>This method uses the {@link DataContext#getDefault() default data context}.
     *
     * @param geoid level surface of the gravity potential of a body
     * @see #GlobalPressureTemperature(Geoid, TimeScale)
     */
    @DefaultDataContext
    public GlobalPressureTemperature(final Geoid geoid) {
        this(geoid, DataContext.getDefault().getTimeScales().getUTC());
    }

    /** Build a new instance.
     * @param geoid level surface of the gravity potential of a body
     * @param utc UTC time scale.
     */
    public GlobalPressureTemperature(final Geoid geoid, final TimeScale utc) {
        this.geoid = geoid;
        this.utc   = utc;
    }

    /** Provide weather parameters.
     * @param location location at which parameters are requested
     * @param date date at which parameters are requested
     * @return weather parameters
     */
    public PressureTemperature getWeatherParameters(final GeodeticPoint location, final AbsoluteDate date) {

        // Day of year computation
        final DateTimeComponents dtc = date.getComponents(utc);
        final int dofyear = dtc.getDate().getDayOfYear();

        // Reference day: 28 January 1980 (Niell, 1996)
        final int t0 = 28;
        final double coef = ((dofyear + 1 - t0) / 365.25) * 2 * FastMath.PI;
        final double cosCoef = FastMath.cos(coef);

        // Compute Legendre Polynomials Pnm(sin(phi))
        final LegendrePolynomials p = new LegendrePolynomials(DEGREE, ORDER,
                                                              FastMath.sin(location.getLatitude()));

        // Corrected height
        final double correctedheight = FastMath.max(0.0,
                                                    location.getAltitude() - geoid.getUndulation(location.getLatitude(),
                                                                                                 location.getLongitude(),
                                                                                                 date));

        // Eq. 4 (Ref)
        double meanT0      = 0.0;
        double amplitudeT0 = 0.0;
        double meanP0      = 0.0;
        double amplitudeP0 = 0.0;
        final ABCoefficients abCoef = new ABCoefficients();
        int j = 0;
        for (int n = 0; n <= DEGREE; n++) {
            for (int m = 0; m <= n; m++) {
                final SinCos sc = FastMath.sinCos(m * location.getLongitude());
                final double pCosmLambda = p.getPnm(n, m) * sc.cos();
                final double pSinmLambda = p.getPnm(n, m) * sc.sin();

                meanT0      = meanT0 +
                                (abCoef.getAnmTemperatureMean(j) * pCosmLambda + abCoef.getBnmTemperatureMean(j) * pSinmLambda);
                amplitudeT0 = amplitudeT0 +
                                (abCoef.getAnmTemperatureAmpl(j) * pCosmLambda + abCoef.getBnmTemperatureAmpl(j) * pSinmLambda);
                meanP0      = meanP0 +
                                (abCoef.getAnmPressureMean(j) * pCosmLambda + abCoef.getBnmPressureMean(j) * pSinmLambda);
                amplitudeP0 = amplitudeP0 +
                                (abCoef.getAnmPressureAmpl(j) * pCosmLambda + abCoef.getBnmPressureAmpl(j) * pSinmLambda);

                j = j + 1;
            }
        }

        // Eq. 3 (Ref)
        final double temp0 = meanT0 + amplitudeT0 * cosCoef;
        final double pres0 = meanP0 + amplitudeP0 * cosCoef;

        // Compute pressure and temperature Eq. 1 and 2 (Ref)
        final double degrees = temp0 + TEMPERATURE_GRADIENT * correctedheight;
        final double temperature = degrees + 273.15;
        final double pressure    = pres0 * FastMath.pow(1.0 - correctedheight * 0.0000226, 5.225);

        return new PressureTemperature(location.getAltitude(),
                                       TroposphericModelUtils.HECTO_PASCAL.toSI(pressure),
                                       temperature);

    }

    private static class ABCoefficients {

        /** Mean A<sub>nm</sub> coefficients for the pressure. */
        private static final double[] A_PRESSURE_MEAN = {
            1.0108e+03,
            8.4886e+00,
            1.4799e+00,
            -1.3897e+01,
            3.7516e-03,
            -1.4936e-01,
            1.2232e+01,
            -7.6615e-01,
            -6.7699e-02,
            8.1002e-03,
            -1.5874e+01,
            3.6614e-01,
            -6.7807e-02,
            -3.6309e-03,
            5.9966e-04,
            4.8163e+00,
            -3.7363e-01,
            -7.2071e-02,
            1.9998e-03,
            -6.2385e-04,
            -3.7916e-04,
            4.7609e+00,
            -3.9534e-01,
            8.6667e-03,
            1.1569e-02,
            1.1441e-03,
            -1.4193e-04,
            -8.5723e-05,
            6.5008e-01,
            -5.0889e-01,
            -1.5754e-02,
            -2.8305e-03,
            5.7458e-04,
            3.2577e-05,
            -9.6052e-06,
            -2.7974e-06,
            1.3530e+00,
            -2.7271e-01,
            -3.0276e-04,
            3.6286e-03,
            -2.0398e-04,
            1.5846e-05,
            -7.7787e-06,
            1.1210e-06,
            9.9020e-08,
            5.5046e-01,
            -2.7312e-01,
            3.2532e-03,
            -2.4277e-03,
            1.1596e-04,
            2.6421e-07,
            -1.3263e-06,
            2.7322e-07,
            1.4058e-07,
            4.9414e-09
        };

        /** Mean B<sub>nm</sub> coefficients for the pressure. */
        private static final double[] B_PRESSURE_MEAN = {
            0.0000e+00,
            0.0000e+00,
            -1.2878e+00,
            0.0000e+00,
            7.0444e-01,
            3.3222e-01,
            0.0000e+00,
            -2.9636e-01,
            7.2248e-03,
            7.9655e-03,
            0.0000e+00,
            1.0854e+00,
            1.1145e-02,
            -3.6513e-02,
            3.1527e-03,
            0.0000e+00,
            -4.8434e-01,
            5.2023e-02,
            -1.3091e-02,
            1.8515e-03,
            1.5422e-04,
            0.0000e+00,
            6.8298e-01,
            2.5261e-03,
            -9.9703e-04,
            -1.0829e-03,
            +1.7688e-04,
            -3.1418e-05,
            +0.0000e+00,
            -3.7018e-01,
            4.3234e-02,
            7.2559e-03,
            3.1516e-04,
            2.0024e-05,
            -8.0581e-06,
            -2.3653e-06,
            0.0000e+00,
            1.0298e-01,
            -1.5086e-02,
            5.6186e-03,
            3.2613e-05,
            4.0567e-05,
            -1.3925e-06,
            -3.6219e-07,
            -2.0176e-08,
            0.0000e+00,
            -1.8364e-01,
            1.8508e-02,
            7.5016e-04,
            -9.6139e-05,
            -3.1995e-06,
            1.3868e-07,
            -1.9486e-07,
            3.0165e-10,
            -6.4376e-10
        };

        /** Amplitude A<sub>nm</sub> coefficients for the pressure. */
        private static final double[] A_PRESSURE_AMPLITUDE = {
            -1.0444e-01,
            1.6618e-01,
            -6.3974e-02,
            1.0922e+00,
            5.7472e-01,
            -3.0277e-01,
            -3.5087e+00,
            7.1264e-03,
            -1.4030e-01,
            3.7050e-02,
            4.0208e-01,
            -3.0431e-01,
            -1.3292e-01,
            4.6746e-03,
            -1.5902e-04,
            2.8624e+00,
            -3.9315e-01,
            -6.4371e-02,
            1.6444e-02,
            -2.3403e-03,
            4.2127e-05,
            1.9945e+00,
            -6.0907e-01,
            -3.5386e-02,
            -1.0910e-03,
            -1.2799e-04,
            4.0970e-05,
            2.2131e-05,
            -5.3292e-01,
            -2.9765e-01,
            -3.2877e-02,
            1.7691e-03,
            5.9692e-05,
            3.1725e-05,
            2.0741e-05,
            -3.7622e-07,
            2.6372e+00,
            -3.1165e-01,
            1.6439e-02,
            2.1633e-04,
            1.7485e-04,
            2.1587e-05,
            6.1064e-06,
            -1.3755e-08,
            -7.8748e-08,
            -5.9152e-01,
            -1.7676e-01,
            8.1807e-03,
            1.0445e-03,
            2.3432e-04,
            9.3421e-06,
            2.8104e-06,
            -1.5788e-07,
            -3.0648e-08,
            2.6421e-10
        };

        /** Amplitude B<sub>nm</sub> coefficients for the pressure. */
        private static final double[] B_PRESSURE_AMPLITUDE = {
            0.0000e+00,
            0.0000e+00,
            9.3340e-01,
            0.0000e+00,
            8.2346e-01,
            2.2082e-01,
            0.0000e+00,
            9.6177e-01,
            -1.5650e-02,
            1.2708e-03,
            0.0000e+00,
            -3.9913e-01,
            2.8020e-02,
            2.8334e-02,
            8.5980e-04,
            0.0000e+00,
            3.0545e-01,
            -2.1691e-02,
            6.4067e-04,
            -3.6528e-05,
            -1.1166e-04,
            0.0000e+00,
            -7.6974e-02,
            -1.8986e-02,
            +5.6896e-03,
            -2.4159e-04,
            -2.3033e-04,
            -9.6783e-06,
            0.0000e+00,
            -1.0218e-01,
            -1.3916e-02,
            -4.1025e-03,
            -5.1340e-05,
            -7.0114e-05,
            -3.3152e-07,
            1.6901e-06,
            0.0000e+00,
            -1.2422e-02,
            +2.5072e-03,
            +1.1205e-03,
            -1.3034e-04,
            -2.3971e-05,
            -2.6622e-06,
            5.7852e-07,
            4.5847e-08,
            0.0000e+00,
            4.4777e-02,
            -3.0421e-03,
            2.6062e-05,
            -7.2421e-05,
            1.9119e-06,
            3.9236e-07,
            2.2390e-07,
            2.9765e-09,
            -4.6452e-09
        };

        /** Mean A<sub>nm</sub> coefficients for the temperature. */
        private static final double[] A_TEMPERATURE_MEAN = {
            1.6257e+01,
            2.1224e+00,
            9.2569e-01,
            -2.5974e+01,
            1.4510e+00,
            9.2468e-02,
            -5.3192e-01,
            2.1094e-01,
            -6.9210e-02,
            -3.4060e-02,
            -4.6569e+00,
            2.6385e-01,
            -3.6093e-02,
            1.0198e-02,
            -1.8783e-03,
            7.4983e-01,
            1.1741e-01,
            3.9940e-02,
            5.1348e-03,
            5.9111e-03,
            8.6133e-06,
            6.3057e-01,
            1.5203e-01,
            3.9702e-02,
            4.6334e-03,
            2.4406e-04,
            1.5189e-04,
            1.9581e-07,
            5.4414e-01,
            3.5722e-01,
            5.2763e-02,
            4.1147e-03,
            -2.7239e-04,
            -5.9957e-05,
            1.6394e-06,
            -7.3045e-07,
            -2.9394e+00,
            5.5579e-02,
            1.8852e-02,
            3.4272e-03,
            -2.3193e-05,
            -2.9349e-05,
            3.6397e-07,
            2.0490e-06,
            -6.4719e-08,
            -5.2225e-01,
            2.0799e-01,
            1.3477e-03,
            3.1613e-04,
            -2.2285e-04,
            -1.8137e-05,
            -1.5177e-07,
            6.1343e-07,
            7.8566e-08,
            1.0749e-09
        };

        /** Mean B<sub>nm</sub> coefficients for the temperature. */
        private static final double[] B_TEMPERATURE_MEAN = {
            0.0000e+00,
            0.0000e+00,
            1.0210e+00,
            0.0000e+00,
            6.0194e-01,
            1.2292e-01,
            0.0000e+00,
            -4.2184e-01,
            1.8230e-01,
            4.2329e-02,
            0.0000e+00,
            9.3312e-02,
            9.5346e-02,
            -1.9724e-03,
            5.8776e-03,
            0.0000e+00,
            -2.0940e-01,
            3.4199e-02,
            -5.7672e-03,
            -2.1590e-03,
            5.6815e-04,
            0.0000e+00,
            2.2858e-01,
            1.2283e-02,
            -9.3679e-03,
            -1.4233e-03,
            -1.5962e-04,
            4.0160e-05,
            0.0000e+00,
            3.6353e-02,
            -9.4263e-04,
            -3.6762e-03,
            5.8608e-05,
            -2.6391e-05,
            3.2095e-06,
            -1.1605e-06,
            0.0000e+00,
            1.6306e-01,
            1.3293e-02,
            -1.1395e-03,
            5.1097e-05,
            3.3977e-05,
            7.6449e-06,
            -1.7602e-07,
            -7.6558e-08,
            0.0000e+00,
            -4.5415e-02,
            -1.8027e-02,
            3.6561e-04,
            -1.1274e-04,
            1.3047e-05,
            2.0001e-06,
            -1.5152e-07,
            -2.7807e-08,
            7.7491e-09
        };

        /** Amplitude A<sub>nm</sub> coefficients for the temperature. */
        private static final double[] A_TEMPERATURE_AMPLITUDE = {
            -1.8654e+00,
            -9.0041e+00,
            -1.2974e-01,
            -3.6053e+00,
            2.0284e-02,
            2.1872e-01,
            -1.3015e+00,
            4.0355e-01,
            2.2216e-01,
            -4.0605e-03,
            1.9623e+00,
            4.2887e-01,
            2.1437e-01,
            -1.0061e-02,
            -1.1368e-03,
            -6.9235e-02,
            5.6758e-01,
            1.1917e-01,
            -7.0765e-03,
            3.0017e-04,
            3.0601e-04,
            1.6559e+00,
            2.0722e-01,
            6.0013e-02,
            1.7023e-04,
            -9.2424e-04,
            1.1269e-05,
            -6.9911e-06,
            -2.0886e+00,
            -6.7879e-02,
            -8.5922e-04,
            -1.6087e-03,
            -4.5549e-05,
            3.3178e-05,
            -6.1715e-06,
            -1.4446e-06,
            -3.7210e-01,
            1.5775e-01,
            -1.7827e-03,
            -4.4396e-04,
            2.2844e-04,
            -1.1215e-05,
            -2.1120e-06,
            -9.6421e-07,
            -1.4170e-08,
            7.8720e-01,
            -4.4238e-02,
            -1.5120e-03,
            -9.4119e-04,
            4.0645e-06,
            -4.9253e-06,
            -1.8656e-06,
            -4.0736e-07,
            -4.9594e-08,
            1.6134e-09
        };

        /** Amplitude B<sub>nm</sub> coefficients for the temperature. */
        private static final double[] B_TEMPERATURE_AMPLITUDE = {
            0.0000e+00,
            0.0000e+00,
            -8.9895e-01,
            0.0000e+00,
            -1.0790e+00,
            -1.2699e-01,
            0.0000e+00,
            -5.9033e-01,
            3.4865e-02,
            -3.2614e-02,
            0.0000e+00,
            -2.4310e-02,
            1.5607e-02,
            -2.9833e-02,
            -5.9048e-03,
            0.0000e+00,
            2.8383e-01,
            4.0509e-02,
            -1.8834e-02,
            -1.2654e-03,
            -1.3794e-04,
            0.0000e+00,
            1.3306e-01,
            3.4960e-02,
            -3.6799e-03,
            -3.5626e-04,
            1.4814e-04,
            3.7932e-06,
            0.0000e+00,
            2.0801e-01,
            6.5640e-03,
            -3.4893e-03,
            -2.7395e-04,
            7.4296e-05,
            -7.9927e-06,
            -1.0277e-06,
            0.0000e+00,
            3.6515e-02,
            -7.4319e-03,
            -6.2873e-04,
            8.2461e-05,
            3.1095e-05,
            -5.3860e-07,
            -1.2055e-07,
            -1.1517e-07,
            0.0000e+00,
            3.1404e-02,
            1.5580e-02,
            -1.1428e-03,
            3.3529e-05,
            1.0387e-05,
            -1.9378e-06,
            -2.7327e-07,
            7.5833e-09,
            -9.2323e-09
        };

        /** Build a new instance. */
        ABCoefficients() {

        }

        /** Get the value of the mean A<sub>nm</sub> pressure coefficient for the given index.
         * @param index index
         * @return the mean A<sub>nm</sub> pressure coefficient for the given index
         */
        public double getAnmPressureMean(final int index) {
            return A_PRESSURE_MEAN[index];
        }

        /** Get the value of the mean B<sub>nm</sub> pressure coefficient for the given index.
         * @param index index
         * @return the mean B<sub>nm</sub> pressure coefficient for the given index
         */
        public double getBnmPressureMean(final int index) {
            return B_PRESSURE_MEAN[index];
        }

        /** Get the value of the amplitude A<sub>nm</sub> pressure coefficient for the given index.
         * @param index index
         * @return the amplitude A<sub>nm</sub> pressure coefficient for the given index.
         */
        public double getAnmPressureAmpl(final int index) {
            return A_PRESSURE_AMPLITUDE[index];
        }

        /** Get the value of the amplitude B<sub>nm</sub> pressure coefficient for the given index.
         * @param index index
         * @return the amplitude B<sub>nm</sub> pressure coefficient for the given index
         */
        public double getBnmPressureAmpl(final int index) {
            return B_PRESSURE_AMPLITUDE[index];
        }

        /** Get the value of the mean A<sub>nm</sub> temperature coefficient for the given index.
         * @param index index
         * @return the mean A<sub>nm</sub> temperature coefficient for the given index
         */
        public double getAnmTemperatureMean(final int index) {
            return A_TEMPERATURE_MEAN[index];
        }

        /** Get the value of the mean B<sub>nm</sub> temperature coefficient for the given index.
         * @param index index
         * @return the mean B<sub>nm</sub> temperature coefficient for the given index
         */
        public double getBnmTemperatureMean(final int index) {
            return B_TEMPERATURE_MEAN[index];
        }

        /** Get the value of the amplitude A<sub>nm</sub> temperature coefficient for the given index.
         * @param index index
         * @return the amplitude A<sub>nm</sub> temperature coefficient for the given index.
         */
        public double getAnmTemperatureAmpl(final int index) {
            return A_TEMPERATURE_AMPLITUDE[index];
        }

        /** Get the value of the amplitude B<sub>nm</sub> temperature coefficient for the given index.
         * @param index index
         * @return the amplitude B<sub>nm</sub> temperature coefficient for the given index
         */
        public double getBnmTemperatureAmpl(final int index) {
            return B_TEMPERATURE_AMPLITUDE[index];
        }
    }

}