SsrVtecIonosphericModel.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.ionosphere;

import java.util.Collections;
import java.util.List;

import org.hipparchus.CalculusFieldElement;
import org.hipparchus.Field;
import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
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.frames.TopocentricFrame;
import org.orekit.gnss.metric.messages.ssr.subtype.SsrIm201;
import org.orekit.gnss.metric.messages.ssr.subtype.SsrIm201Data;
import org.orekit.gnss.metric.messages.ssr.subtype.SsrIm201Header;
import org.orekit.propagation.FieldSpacecraftState;
import org.orekit.propagation.SpacecraftState;
import org.orekit.utils.Constants;
import org.orekit.utils.FieldLegendrePolynomials;
import org.orekit.utils.LegendrePolynomials;
import org.orekit.utils.ParameterDriver;

/**
 * Ionospheric model based on SSR IM201 message.
 * <p>
 * Within this message, the ionospheric VTEC is provided
 * using spherical harmonic expansions. For a given ionospheric
 * layer, the slant TEC value is calculated using the satellite
 * elevation and the height of the corresponding layer. The
 * total slant TEC is computed by the sum of the individual slant
 * TEC for each layer.
 * </p>
 * @author Bryan Cazabonne
 * @since 11.0
 * @see "IGS State Space Representation (SSR) Format, Version 1.00, October 2020."
 */
public class SsrVtecIonosphericModel implements IonosphericModel {

    /** Earth radius in meters (see reference). */
    private static final double EARTH_RADIUS = 6370000.0;

    /** Multiplication factor for path delay computation. */
    private static final double FACTOR = 40.3e16;

    /** SSR Ionosphere VTEC Spherical Harmonics Message.. */
    private final transient SsrIm201 vtecMessage;

    /**
     * Constructor.
     * @param vtecMessage SSR Ionosphere VTEC Spherical Harmonics Message.
     */
    public SsrVtecIonosphericModel(final SsrIm201 vtecMessage) {
        this.vtecMessage = vtecMessage;
    }

    /** {@inheritDoc} */
    @Override
    public double pathDelay(final SpacecraftState state, final TopocentricFrame baseFrame,
                            final double frequency, final double[] parameters) {

        // Elevation in radians
        final Vector3D position  = state.getPosition(baseFrame);
        final double   elevation = position.getDelta();

        // Only consider measures above the horizon
        if (elevation > 0.0) {

            // Azimuth angle in radians
            double azimuth = FastMath.atan2(position.getX(), position.getY());
            if (azimuth < 0.) {
                azimuth += MathUtils.TWO_PI;
            }

            // Initialize slant TEC
            double stec = 0.0;

            // Message header
            final SsrIm201Header header = vtecMessage.getHeader();

            // Loop on ionospheric layers
            for (final SsrIm201Data data : vtecMessage.getData()) {
                stec += stecIonosphericLayer(data, header, elevation, azimuth, baseFrame.getPoint());
            }

            // Return the path delay
            return FACTOR * stec / (frequency * frequency);

        }

        // Delay is equal to 0.0
        return 0.0;

    }

    /** {@inheritDoc} */
    @Override
    public <T extends CalculusFieldElement<T>> T pathDelay(final FieldSpacecraftState<T> state, final TopocentricFrame baseFrame,
                                                       final double frequency, final T[] parameters) {

        // Field
        final Field<T> field = state.getDate().getField();

        // Elevation in radians
        final FieldVector3D<T> position  = state.getPosition(baseFrame);
        final T                elevation = position.getDelta();

        // Only consider measures above the horizon
        if (elevation.getReal() > 0.0) {

            // Azimuth angle in radians
            T azimuth = FastMath.atan2(position.getX(), position.getY());
            if (azimuth.getReal() < 0.) {
                azimuth = azimuth.add(MathUtils.TWO_PI);
            }

            // Initialize slant TEC
            T stec = field.getZero();

            // Message header
            final SsrIm201Header header = vtecMessage.getHeader();

            // Loop on ionospheric layers
            for (SsrIm201Data data : vtecMessage.getData()) {
                stec = stec.add(stecIonosphericLayer(data, header, elevation, azimuth, baseFrame.getPoint(field)));
            }

            // Return the path delay
            return stec.multiply(FACTOR).divide(frequency * frequency);

        }

        // Delay is equal to 0.0
        return field.getZero();

    }

    /** {@inheritDoc} */
    @Override
    public List<ParameterDriver> getParametersDrivers() {
        return Collections.emptyList();
    }

    /**
     * Calculates the slant TEC for a given ionospheric layer.
     * @param im201Data ionospheric data for the current layer
     * @param im201Header container for data contained in the header
     * @param elevation satellite elevation angle [rad]
     * @param azimuth satellite azimuth angle [rad]
     * @param point geodetic point
     * @return the slant TEC for the current ionospheric layer
     */
    private static double stecIonosphericLayer(final SsrIm201Data im201Data, final SsrIm201Header im201Header,
                                               final double elevation, final double azimuth,
                                               final GeodeticPoint point) {

        // Geodetic point data
        final double phiR    = point.getLatitude();
        final double lambdaR = point.getLongitude();
        final double hR      = point.getAltitude();

        // Data contained in the message
        final double     hI     = im201Data.getHeightIonosphericLayer();
        final int        degree = im201Data.getSphericalHarmonicsDegree();
        final int        order  = im201Data.getSphericalHarmonicsOrder();
        final double[][] cnm    = im201Data.getCnm();
        final double[][] snm    = im201Data.getSnm();

        // Spherical Earth's central angle
        final double psiPP = calculatePsi(hR, hI, elevation);

        // Sine and cosine of useful angles
        final SinCos scA    = FastMath.sinCos(azimuth);
        final SinCos scPhiR = FastMath.sinCos(phiR);
        final SinCos scPsi  = FastMath.sinCos(psiPP);

        // Pierce point latitude and longitude
        final double phiPP    = calculatePiercePointLatitude(scPhiR, scPsi, scA);
        final double lambdaPP = calculatePiercePointLongitude(scA, phiPP, psiPP, phiR, lambdaR);

        // Mean sun fixed longitude (modulo 2pi)
        final double lambdaS = calculateSunLongitude(im201Header, lambdaPP);

        // VTEC
        // According to the documentation, negative VTEC values must be ignored and shall be replaced by 0.0
        final double vtec = FastMath.max(0.0, calculateVTEC(degree, order, cnm, snm, phiPP, lambdaS));

        // Return STEC for the current ionospheric layer
        return vtec / FastMath.sin(elevation + psiPP);

    }

    /**
     * Calculates the slant TEC for a given ionospheric layer.
     * @param im201Data ionospheric data for the current layer
     * @param im201Header container for data contained in the header
     * @param elevation satellite elevation angle [rad]
     * @param azimuth satellite azimuth angle [rad]
     * @param point geodetic point
     * @param <T> type of the elements
     * @return the slant TEC for the current ionospheric layer
     */
    private static <T extends CalculusFieldElement<T>> T stecIonosphericLayer(final SsrIm201Data im201Data, final SsrIm201Header im201Header,
                                                                          final T elevation, final T azimuth,
                                                                          final FieldGeodeticPoint<T> point) {

        // Geodetic point data
        final T phiR    = point.getLatitude();
        final T lambdaR = point.getLongitude();
        final T hR      = point.getAltitude();

        // Data contained in the message
        final double     hI     = im201Data.getHeightIonosphericLayer();
        final int        degree = im201Data.getSphericalHarmonicsDegree();
        final int        order  = im201Data.getSphericalHarmonicsOrder();
        final double[][] cnm    = im201Data.getCnm();
        final double[][] snm    = im201Data.getSnm();

        // Spherical Earth's central angle
        final T psiPP = calculatePsi(hR, hI, elevation);

        // Sine and cosine of useful angles
        final FieldSinCos<T> scA    = FastMath.sinCos(azimuth);
        final FieldSinCos<T> scPhiR = FastMath.sinCos(phiR);
        final FieldSinCos<T> scPsi  = FastMath.sinCos(psiPP);

        // Pierce point latitude and longitude
        final T phiPP    = calculatePiercePointLatitude(scPhiR, scPsi, scA);
        final T lambdaPP = calculatePiercePointLongitude(scA, phiPP, psiPP, phiR, lambdaR);

        // Mean sun fixed longitude (modulo 2pi)
        final T lambdaS = calculateSunLongitude(im201Header, lambdaPP);

        // VTEC
        // According to the documentation, negative VTEC values must be ignored and shall be replaced by 0.0
        final T vtec = FastMath.max(phiR.getField().getZero(), calculateVTEC(degree, order, cnm, snm, phiPP, lambdaS));

        // Return STEC for the current ionospheric layer
        return vtec.divide(FastMath.sin(elevation.add(psiPP)));

    }

    /**
     * Calculates the spherical Earth’s central angle between station position and
     * the projection of the pierce point to the spherical Earth’s surface.
     * @param hR height of station position in meters
     * @param hI height of ionospheric layer in meters
     * @param elevation satellite elevation angle in radians
     * @return the spherical Earth’s central angle in radians
     */
    private static double calculatePsi(final double hR, final double hI,
                                       final double elevation) {
        final double ratio = (EARTH_RADIUS + hR) / (EARTH_RADIUS + hI);
        return MathUtils.SEMI_PI - elevation - FastMath.asin(ratio * FastMath.cos(elevation));
    }

    /**
     * Calculates the spherical Earth’s central angle between station position and
     * the projection of the pierce point to the spherical Earth’s surface.
     * @param hR height of station position in meters
     * @param hI height of ionospheric layer in meters
     * @param elevation satellite elevation angle in radians
     * @param <T> type of the elements
     * @return the spherical Earth’s central angle in radians
     */
    private static <T extends CalculusFieldElement<T>> T calculatePsi(final T hR, final double hI,
                                                                  final T elevation) {
        final T ratio = hR.add(EARTH_RADIUS).divide(EARTH_RADIUS + hI);
        return hR.getPi().multiply(0.5).subtract(elevation).subtract(FastMath.asin(ratio.multiply(FastMath.cos(elevation))));
    }

    /**
     * Calculates the latitude of the pierce point in the spherical Earth model.
     * @param scPhiR sine and cosine of the geocentric latitude of the station
     * @param scPsi sine and cosine of the spherical Earth's central angle
     * @param scA sine and cosine of the azimuth angle
     * @return the latitude of the pierce point in the spherical Earth model in radians
     */
    private static double calculatePiercePointLatitude(final SinCos scPhiR, final SinCos scPsi, final SinCos scA) {
        return FastMath.asin(scPhiR.sin() * scPsi.cos() + scPhiR.cos() * scPsi.sin() * scA.cos());
    }

    /**
     * Calculates the latitude of the pierce point in the spherical Earth model.
     * @param scPhiR sine and cosine of the geocentric latitude of the station
     * @param scPsi sine and cosine of the spherical Earth's central angle
     * @param scA sine and cosine of the azimuth angle
     * @param <T> type of the elements
     * @return the latitude of the pierce point in the spherical Earth model in radians
     */
    private static <T extends CalculusFieldElement<T>> T calculatePiercePointLatitude(final FieldSinCos<T> scPhiR,
                                                                                  final FieldSinCos<T> scPsi,
                                                                                  final FieldSinCos<T> scA) {
        return FastMath.asin(scPhiR.sin().multiply(scPsi.cos()).add(scPhiR.cos().multiply(scPsi.sin()).multiply(scA.cos())));
    }

    /**
     * Calculates the longitude of the pierce point in the spherical Earth model.
     * @param scA sine and cosine of the azimuth angle
     * @param phiPP the latitude of the pierce point in the spherical Earth model in radians
     * @param psiPP the spherical Earth’s central angle in radians
     * @param phiR the geocentric latitude of the station in radians
     * @param lambdaR the geocentric longitude of the station
     * @return the longitude of the pierce point in the spherical Earth model in radians
     */
    private static double calculatePiercePointLongitude(final SinCos scA,
                                                        final double phiPP, final double psiPP,
                                                        final double phiR, final double lambdaR) {

        // arcSin(sin(PsiPP) * sin(Azimuth) / cos(PhiPP))
        final double arcSin = FastMath.asin(FastMath.sin(psiPP) * scA.sin() / FastMath.cos(phiPP));

        // Return
        return verifyCondition(scA.cos(), psiPP, phiR) ? lambdaR + FastMath.PI - arcSin : lambdaR + arcSin;

    }

    /**
     * Calculates the longitude of the pierce point in the spherical Earth model.
     * @param scA sine and cosine of the azimuth angle
     * @param phiPP the latitude of the pierce point in the spherical Earth model in radians
     * @param psiPP the spherical Earth’s central angle in radians
     * @param phiR the geocentric latitude of the station in radians
     * @param lambdaR the geocentric longitude of the station
     * @param <T> type of the elements
     * @return the longitude of the pierce point in the spherical Earth model in radians
     */
    private static <T extends CalculusFieldElement<T>> T calculatePiercePointLongitude(final FieldSinCos<T> scA,
                                                                                   final T phiPP, final T psiPP,
                                                                                   final T phiR, final T lambdaR) {

        // arcSin(sin(PsiPP) * sin(Azimuth) / cos(PhiPP))
        final T arcSin = FastMath.asin(FastMath.sin(psiPP).multiply(scA.sin()).divide(FastMath.cos(phiPP)));

        // Return
        return verifyCondition(scA.cos().getReal(), psiPP.getReal(), phiR.getReal()) ?
                                               lambdaR.add(arcSin.getPi()).subtract(arcSin) : lambdaR.add(arcSin);

    }

    /**
     * Calculate the mean sun fixed longitude phase.
     * @param im201Header header of the IM201 message
     * @param lambdaPP the longitude of the pierce point in the spherical Earth model in radians
     * @return the mean sun fixed longitude phase in radians
     */
    private static double calculateSunLongitude(final SsrIm201Header im201Header, final double lambdaPP) {
        final double t = getTime(im201Header);
        return MathUtils.normalizeAngle(lambdaPP + (t - 50400.0) * FastMath.PI / 43200.0, FastMath.PI);
    }

    /**
     * Calculate the mean sun fixed longitude phase.
     * @param im201Header header of the IM201 message
     * @param lambdaPP the longitude of the pierce point in the spherical Earth model in radians
     * @param <T> type of the elements
     * @return the mean sun fixed longitude phase in radians
     */
    private static <T extends CalculusFieldElement<T>> T calculateSunLongitude(final SsrIm201Header im201Header, final T lambdaPP) {
        final double t = getTime(im201Header);
        return MathUtils.normalizeAngle(lambdaPP.add(lambdaPP.getPi().multiply(t - 50400.0).divide(43200.0)), lambdaPP.getPi());
    }

    /**
     * Calculate the VTEC contribution for a given ionospheric layer.
     * @param degree degree of spherical expansion
     * @param order order of spherical expansion
     * @param cnm cosine coefficients for the layer in TECU
     * @param snm sine coefficients for the layer in TECU
     * @param phiPP geocentric latitude of ionospheric pierce point for the layer in radians
     * @param lambdaS mean sun fixed and phase shifted longitude of ionospheric pierce point
     * @return the VTEC contribution for the current ionospheric layer in TECU
     */
    private static double calculateVTEC(final int degree, final int order,
                                        final double[][] cnm, final double[][] snm,
                                        final double phiPP, final double lambdaS) {

        // Initialize VTEC value
        double vtec = 0.0;

        // Compute Legendre Polynomials Pnm(sin(phiPP))
        final LegendrePolynomials p = new LegendrePolynomials(degree, order, FastMath.sin(phiPP));

        // Calculate VTEC
        for (int n = 0; n <= degree; n++) {

            for (int m = 0; m <= FastMath.min(n, order); m++) {

                // Legendre coefficients
                final SinCos sc = FastMath.sinCos(m * lambdaS);
                final double pCosmLambda = p.getPnm(n, m) * sc.cos();
                final double pSinmLambda = p.getPnm(n, m) * sc.sin();

                // Update VTEC value
                vtec += cnm[n][m] * pCosmLambda + snm[n][m] * pSinmLambda;

            }

        }

        // Return the VTEC
        return vtec;

    }

    /**
     * Calculate the VTEC contribution for a given ionospheric layer.
     * @param degree degree of spherical expansion
     * @param order order of spherical expansion
     * @param cnm cosine coefficients for the layer in TECU
     * @param snm sine coefficients for the layer in TECU
     * @param phiPP geocentric latitude of ionospheric pierce point for the layer in radians
     * @param lambdaS mean sun fixed and phase shifted longitude of ionospheric pierce point
     * @param <T> type of the elements
     * @return the VTEC contribution for the current ionospheric layer in TECU
     */
    private static <T extends CalculusFieldElement<T>> T calculateVTEC(final int degree, final int order,
                                                                   final double[][] cnm, final double[][] snm,
                                                                   final T phiPP, final T lambdaS) {

        // Initialize VTEC value
        T vtec = phiPP.getField().getZero();

        // Compute Legendre Polynomials Pnm(sin(phiPP))
        final FieldLegendrePolynomials<T> p = new FieldLegendrePolynomials<>(degree, order, FastMath.sin(phiPP));

        // Calculate VTEC
        for (int n = 0; n <= degree; n++) {

            for (int m = 0; m <= FastMath.min(n, order); m++) {

                // Legendre coefficients
                final FieldSinCos<T> sc = FastMath.sinCos(lambdaS.multiply(m));
                final T pCosmLambda = p.getPnm(n, m).multiply(sc.cos());
                final T pSinmLambda = p.getPnm(n, m).multiply(sc.sin());

                // Update VTEC value
                vtec = vtec.add(pCosmLambda.multiply(cnm[n][m]).add(pSinmLambda.multiply(snm[n][m])));

            }

        }

        // Return the VTEC
        return vtec;

    }

    /**
     * Get the SSR epoch time of computation modulo 86400 seconds.
     * @param im201Header header data
     * @return the SSR epoch time of computation modulo 86400 seconds
     */
    private static double getTime(final SsrIm201Header im201Header) {
        final double ssrEpochTime = im201Header.getSsrEpoch1s();
        return ssrEpochTime - FastMath.floor(ssrEpochTime / Constants.JULIAN_DAY) * Constants.JULIAN_DAY;
    }

    /**
     * Verify the condition for the calculation of the pierce point longitude.
     * @param scACos cosine of the azimuth angle
     * @param psiPP the spherical Earth’s central angle in radians
     * @param phiR the geocentric latitude of the station in radians
     * @return true if the condition is respected
     */
    private static boolean verifyCondition(final double scACos, final double psiPP,
                                           final double phiR) {

        // tan(PsiPP) * cos(Azimuth)
        final double tanPsiCosA = FastMath.tan(psiPP) * scACos;

        // Verify condition
        return phiR >= 0 && tanPsiCosA > FastMath.tan(MathUtils.SEMI_PI - phiR) ||
                        phiR < 0 && -tanPsiCosA > FastMath.tan(MathUtils.SEMI_PI + phiR);

    }

}