SsrVtecIonosphericModel.java

  1. /* Copyright 2002-2025 CS GROUP
  2.  * Licensed to CS GROUP (CS) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * CS licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *   http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.orekit.models.earth.ionosphere;

  18. import java.util.Collections;
  19. import java.util.List;

  20. import org.hipparchus.CalculusFieldElement;
  21. import org.hipparchus.Field;
  22. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  23. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  24. import org.hipparchus.util.FastMath;
  25. import org.hipparchus.util.FieldSinCos;
  26. import org.hipparchus.util.MathUtils;
  27. import org.hipparchus.util.SinCos;
  28. import org.orekit.bodies.FieldGeodeticPoint;
  29. import org.orekit.bodies.GeodeticPoint;
  30. import org.orekit.frames.FieldStaticTransform;
  31. import org.orekit.frames.StaticTransform;
  32. import org.orekit.frames.TopocentricFrame;
  33. import org.orekit.gnss.metric.messages.ssr.subtype.SsrIm201;
  34. import org.orekit.gnss.metric.messages.ssr.subtype.SsrIm201Data;
  35. import org.orekit.gnss.metric.messages.ssr.subtype.SsrIm201Header;
  36. import org.orekit.propagation.FieldSpacecraftState;
  37. import org.orekit.propagation.SpacecraftState;
  38. import org.orekit.time.AbsoluteDate;
  39. import org.orekit.time.FieldAbsoluteDate;
  40. import org.orekit.utils.Constants;
  41. import org.orekit.utils.FieldLegendrePolynomials;
  42. import org.orekit.utils.LegendrePolynomials;
  43. import org.orekit.utils.ParameterDriver;

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

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

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

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

  65.     /**
  66.      * Constructor.
  67.      * @param vtecMessage SSR Ionosphere VTEC Spherical Harmonics Message.
  68.      */
  69.     public SsrVtecIonosphericModel(final SsrIm201 vtecMessage) {
  70.         this.vtecMessage = vtecMessage;
  71.     }

  72.     /** {@inheritDoc} */
  73.     @Override
  74.     public double pathDelay(final SpacecraftState state, final TopocentricFrame baseFrame,
  75.                             final double frequency, final double[] parameters) {
  76.         return pathDelay(state, baseFrame, state.getDate(), frequency, parameters);
  77.     }

  78.     /** {@inheritDoc} */
  79.     @Override
  80.     public double pathDelay(final SpacecraftState state,
  81.                             final TopocentricFrame baseFrame, final AbsoluteDate receptionDate,
  82.                             final double frequency, final double[] parameters) {

  83.         // we use transform from base frame to inert frame and invert it
  84.         // because base frame could be peered with inertial frame (hence improving performances with caching)
  85.         // but the reverse is almost never used
  86.         final StaticTransform base2Inert = baseFrame.getStaticTransformTo(state.getFrame(), receptionDate);
  87.         final Vector3D        position   = base2Inert.getInverse().transformPosition(state.getPosition());

  88.         // Elevation in radians
  89.         final double   elevation = position.getDelta();

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

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

  97.             // Initialize slant TEC
  98.             double stec = 0.0;

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

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

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

  107.         }

  108.         // Delay is equal to 0.0
  109.         return 0.0;

  110.     }

  111.     /** {@inheritDoc} */
  112.     @Override
  113.     public <T extends CalculusFieldElement<T>> T pathDelay(final FieldSpacecraftState<T> state, final TopocentricFrame baseFrame,
  114.                                                        final double frequency, final T[] parameters) {
  115.         return pathDelay(state, baseFrame, state.getDate(), frequency, parameters);
  116.     }

  117.     /** {@inheritDoc} */
  118.     @Override
  119.     public <T extends CalculusFieldElement<T>> T pathDelay(final FieldSpacecraftState<T> state,
  120.                                                            final TopocentricFrame baseFrame, final FieldAbsoluteDate<T> receptionDate,
  121.                                                            final double frequency, final T[] parameters) {

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

  124.         // we use transform from base frame to inert frame and invert it
  125.         // because base frame could be peered with inertial frame (hence improving performances with caching)
  126.         // but the reverse is almost never used
  127.         final FieldStaticTransform<T> base2Inert = baseFrame.getStaticTransformTo(state.getFrame(), receptionDate);
  128.         final FieldVector3D<T>        position   = base2Inert.getInverse().transformPosition(state.getPosition());

  129.         // Elevation in radians
  130.         final T                elevation = position.getDelta();

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

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

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

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

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

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

  148.         }

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

  151.     }

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

  157.     /**
  158.      * Calculates the slant TEC for a given ionospheric layer.
  159.      * @param im201Data ionospheric data for the current layer
  160.      * @param im201Header container for data contained in the header
  161.      * @param elevation satellite elevation angle [rad]
  162.      * @param azimuth satellite azimuth angle [rad]
  163.      * @param point geodetic point
  164.      * @return the slant TEC for the current ionospheric layer
  165.      */
  166.     private static double stecIonosphericLayer(final SsrIm201Data im201Data, final SsrIm201Header im201Header,
  167.                                                final double elevation, final double azimuth,
  168.                                                final GeodeticPoint point) {

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

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

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

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

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

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

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

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

  195.     }

  196.     /**
  197.      * Calculates the slant TEC for a given ionospheric layer.
  198.      * @param im201Data ionospheric data for the current layer
  199.      * @param im201Header container for data contained in the header
  200.      * @param elevation satellite elevation angle [rad]
  201.      * @param azimuth satellite azimuth angle [rad]
  202.      * @param point geodetic point
  203.      * @param <T> type of the elements
  204.      * @return the slant TEC for the current ionospheric layer
  205.      */
  206.     private static <T extends CalculusFieldElement<T>> T stecIonosphericLayer(final SsrIm201Data im201Data, final SsrIm201Header im201Header,
  207.                                                                           final T elevation, final T azimuth,
  208.                                                                           final FieldGeodeticPoint<T> point) {

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

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

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

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

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

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

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

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

  235.     }

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

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

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

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

  286.     /**
  287.      * Calculates the longitude of the pierce point in the spherical Earth model.
  288.      * @param scA sine and cosine of the azimuth angle
  289.      * @param phiPP the latitude of the pierce point in the spherical Earth model in radians
  290.      * @param psiPP the spherical Earth’s central angle in radians
  291.      * @param phiR the geocentric latitude of the station in radians
  292.      * @param lambdaR the geocentric longitude of the station
  293.      * @return the longitude of the pierce point in the spherical Earth model in radians
  294.      */
  295.     private static double calculatePiercePointLongitude(final SinCos scA,
  296.                                                         final double phiPP, final double psiPP,
  297.                                                         final double phiR, final double lambdaR) {

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

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

  302.     }

  303.     /**
  304.      * Calculates the longitude of the pierce point in the spherical Earth model.
  305.      * @param scA sine and cosine of the azimuth angle
  306.      * @param phiPP the latitude of the pierce point in the spherical Earth model in radians
  307.      * @param psiPP the spherical Earth’s central angle in radians
  308.      * @param phiR the geocentric latitude of the station in radians
  309.      * @param lambdaR the geocentric longitude of the station
  310.      * @param <T> type of the elements
  311.      * @return the longitude of the pierce point in the spherical Earth model in radians
  312.      */
  313.     private static <T extends CalculusFieldElement<T>> T calculatePiercePointLongitude(final FieldSinCos<T> scA,
  314.                                                                                    final T phiPP, final T psiPP,
  315.                                                                                    final T phiR, final T lambdaR) {

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

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

  321.     }

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

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

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

  356.         // Initialize VTEC value
  357.         double vtec = 0.0;

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

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

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

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

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

  369.             }

  370.         }

  371.         // Return the VTEC
  372.         return vtec;

  373.     }

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

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

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

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

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

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

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

  401.             }

  402.         }

  403.         // Return the VTEC
  404.         return vtec;

  405.     }

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

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

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

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

  429.     }

  430. }