KlobucharIonoModel.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.Field;
  21. import org.hipparchus.CalculusFieldElement;
  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.annotation.DefaultDataContext;
  29. import org.orekit.bodies.FieldGeodeticPoint;
  30. import org.orekit.bodies.GeodeticPoint;
  31. import org.orekit.data.DataContext;
  32. import org.orekit.frames.FieldStaticTransform;
  33. import org.orekit.frames.StaticTransform;
  34. import org.orekit.frames.TopocentricFrame;
  35. import org.orekit.propagation.FieldSpacecraftState;
  36. import org.orekit.propagation.SpacecraftState;
  37. import org.orekit.time.AbsoluteDate;
  38. import org.orekit.time.DateTimeComponents;
  39. import org.orekit.time.FieldAbsoluteDate;
  40. import org.orekit.time.TimeScale;
  41. import org.orekit.utils.Constants;
  42. import org.orekit.utils.ParameterDriver;

  43. /**
  44.  * Klobuchar ionospheric delay model.
  45.  * Klobuchar ionospheric delay model is designed as a GNSS correction model.
  46.  * The parameters for the model are provided by the GPS satellites in their broadcast
  47.  * messsage.
  48.  * This model is based on the assumption the electron content is concentrated
  49.  * in 350 km layer.
  50.  *
  51.  * The delay refers to L1 (1575.42 MHz).
  52.  * If the delay is sought for L2 (1227.60 MHz), multiply the result by 1.65 (Klobuchar, 1996).
  53.  * More generally, since ionospheric delay is inversely proportional to the square of the signal
  54.  * frequency f, to adapt this model to other GNSS frequencies f, multiply by (L1 / f)^2.
  55.  *
  56.  * References:
  57.  *     ICD-GPS-200, Rev. C, (1997), pp. 125-128
  58.  *     Klobuchar, J.A., Ionospheric time-delay algorithm for single-frequency GPS users,
  59.  *         IEEE Transactions on Aerospace and Electronic Systems, Vol. 23, No. 3, May 1987
  60.  *     Klobuchar, J.A., "Ionospheric Effects on GPS", Global Positioning System: Theory and
  61.  *         Applications, 1996, pp.513-514, Parkinson, Spilker.
  62.  *
  63.  * @author Joris Olympio
  64.  * @since 7.1
  65.  *
  66.  */
  67. public class KlobucharIonoModel implements IonosphericModel, IonosphericDelayModel {

  68.     /** The 4 coefficients of a cubic equation representing the amplitude of the vertical delay. Units are sec/semi-circle^(i-1) for the i-th coefficient, i=1, 2, 3, 4. */
  69.     private final double[] alpha;

  70.     /** The 4 coefficients of a cubic equation representing the period of the model. Units are sec/semi-circle^(i-1) for the i-th coefficient, i=1, 2, 3, 4. */
  71.     private final double[] beta;

  72.     /** GPS time scale. */
  73.     private final TimeScale gps;

  74.     /** Create a new Klobuchar ionospheric delay model, when a single frequency system is used.
  75.      * This model accounts for at least 50 percent of RMS error due to ionospheric propagation effect (ICD-GPS-200)
  76.      *
  77.      * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
  78.      *
  79.      * @param alpha coefficients of a cubic equation representing the amplitude of the vertical delay.
  80.      * @param beta coefficients of a cubic equation representing the period of the model.
  81.      * @see #KlobucharIonoModel(double[], double[], TimeScale)
  82.      */
  83.     @DefaultDataContext
  84.     public KlobucharIonoModel(final double[] alpha, final double[] beta) {
  85.         this(alpha, beta, DataContext.getDefault().getTimeScales().getGPS());
  86.     }

  87.     /**
  88.      * Create a new Klobuchar ionospheric delay model, when a single frequency system is
  89.      * used. This model accounts for at least 50 percent of RMS error due to ionospheric
  90.      * propagation effect (ICD-GPS-200)
  91.      *
  92.      * @param alpha coefficients of a cubic equation representing the amplitude of the
  93.      *              vertical delay.
  94.      * @param beta  coefficients of a cubic equation representing the period of the
  95.      *              model.
  96.      * @param gps   GPS time scale.
  97.      * @since 10.1
  98.      */
  99.     public KlobucharIonoModel(final double[] alpha,
  100.                               final double[] beta,
  101.                               final TimeScale gps) {
  102.         this.alpha = alpha.clone();
  103.         this.beta  = beta.clone();
  104.         this.gps = gps;
  105.     }

  106.     /**
  107.      * Calculates the ionospheric path delay for the signal path from a ground
  108.      * station to a satellite.
  109.      * <p>
  110.      * The path delay is computed for any elevation angle.
  111.      * </p>
  112.      * @param date current date
  113.      * @param geo geodetic point of receiver/station
  114.      * @param elevation elevation of the satellite in radians
  115.      * @param azimuth azimuth of the satellite in radians
  116.      * @param frequency frequency of the signal in Hz
  117.      * @param parameters ionospheric model parameters
  118.      * @return the path delay due to the ionosphere in m
  119.      */
  120.     public double pathDelay(final AbsoluteDate date, final GeodeticPoint geo,
  121.                             final double elevation, final double azimuth, final double frequency,
  122.                             final double[] parameters) {

  123.         // Sine and cosine of the azimuth
  124.         final SinCos sc = FastMath.sinCos(azimuth);

  125.         // degrees to semicircles
  126.         final double rad2semi = 1. / FastMath.PI;
  127.         final double semi2rad = FastMath.PI;

  128.         // Earth Centered angle
  129.         final double psi = 0.0137 / (elevation / FastMath.PI + 0.11) - 0.022;

  130.         // Subionospheric latitude: the latitude of the IPP (Ionospheric Pierce Point)
  131.         // in [-0.416, 0.416], semicircle
  132.         final double latIono = FastMath.min(
  133.                                       FastMath.max(geo.getLatitude() * rad2semi + psi * sc.cos(), -0.416),
  134.                                       0.416);

  135.         // Subionospheric longitude: the longitude of the IPP
  136.         // in semicircle
  137.         final double lonIono = geo.getLongitude() * rad2semi + (psi * sc.sin() / FastMath.cos(latIono * semi2rad));

  138.         // Geomagnetic latitude, semicircle
  139.         final double latGeom = latIono + 0.064 * FastMath.cos((lonIono - 1.617) * semi2rad);

  140.         // day of week and tow (sec)
  141.         // Note: Sunday=0, Monday=1, Tuesday=2, Wednesday=3, Thursday=4, Friday=5, Saturday=6
  142.         final DateTimeComponents dtc = date.getComponents(gps);
  143.         final int dofweek = dtc.getDate().getDayOfWeek();
  144.         final double secday = dtc.getTime().getSecondsInLocalDay();
  145.         final double tow = dofweek * 86400. + secday;

  146.         final double t = 43200. * lonIono + tow;
  147.         final double tsec = t - FastMath.floor(t / 86400.) * 86400; // Seconds of day

  148.         // Slant factor, semicircle
  149.         final double slantFactor = 1.0 + 16.0 * FastMath.pow(0.53 - elevation / FastMath.PI, 3);

  150.         // Period of model, seconds
  151.         final double period = FastMath.max(72000., beta[0] + (beta[1]  + (beta[2] + beta[3] * latGeom) * latGeom) * latGeom);

  152.         // Phase of the model, radians
  153.         // (Max at 14.00 = 50400 sec local time)
  154.         final double x = 2.0 * FastMath.PI * (tsec - 50400.0) / period;

  155.         // Amplitude of the model, seconds
  156.         final double amplitude = FastMath.max(0, alpha[0] + (alpha[1]  + (alpha[2] + alpha[3] * latGeom) * latGeom) * latGeom);

  157.         // Ionospheric correction (L1)
  158.         double ionoTimeDelayL1 = slantFactor * (5. * 1e-9);
  159.         if (FastMath.abs(x) < 1.570) {
  160.             ionoTimeDelayL1 += slantFactor * (amplitude * (1.0 - FastMath.pow(x, 2) / 2.0 + FastMath.pow(x, 4) / 24.0));
  161.         }

  162.         // Ionospheric delay for the L1 frequency, in meters, with slant correction.
  163.         final double ratio = FastMath.pow(1575.42e6 / frequency, 2);
  164.         return ratio * Constants.SPEED_OF_LIGHT * ionoTimeDelayL1;
  165.     }

  166.     /** {@inheritDoc} */
  167.     @Override
  168.     public double pathDelay(final SpacecraftState state, final TopocentricFrame baseFrame,
  169.                             final double frequency, final double[] parameters) {
  170.         return pathDelay(state, baseFrame, state.getDate(), frequency, parameters);
  171.     }

  172.     /** {@inheritDoc} */
  173.     @Override
  174.     public double pathDelay(final SpacecraftState state,
  175.                             final TopocentricFrame baseFrame, final AbsoluteDate receptionDate,
  176.                             final double frequency, final double[] parameters) {

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

  182.         // Elevation in radians
  183.         final double   elevation = position.getDelta();

  184.         // Only consider measures above the horizon
  185.         if (elevation > 0.0) {
  186.             // Date
  187.             final AbsoluteDate date = state.getDate();
  188.             // Geodetic point
  189.             final GeodeticPoint geo = baseFrame.getPoint();
  190.             // Azimuth angle in radians
  191.             double azimuth = FastMath.atan2(position.getX(), position.getY());
  192.             if (azimuth < 0.) {
  193.                 azimuth += MathUtils.TWO_PI;
  194.             }
  195.             // Delay
  196.             return pathDelay(date, geo, elevation, azimuth, frequency, parameters);
  197.         }

  198.         return 0.0;
  199.     }

  200.     /**
  201.      * Calculates the ionospheric path delay for the signal path from a ground
  202.      * station to a satellite.
  203.      * <p>
  204.      * The path delay is computed for any elevation angle.
  205.      * </p>
  206.      * @param <T> type of the elements
  207.      * @param date current date
  208.      * @param geo geodetic point of receiver/station
  209.      * @param elevation elevation of the satellite in radians
  210.      * @param azimuth azimuth of the satellite in radians
  211.      * @param frequency frequency of the signal in Hz
  212.      * @param parameters ionospheric model parameters
  213.      * @return the path delay due to the ionosphere in m
  214.      */
  215.     public <T extends CalculusFieldElement<T>> T pathDelay(final FieldAbsoluteDate<T> date, final FieldGeodeticPoint<T> geo,
  216.                                                        final T elevation, final T azimuth, final double frequency,
  217.                                                        final T[] parameters) {

  218.         // Sine and cosine of the azimuth
  219.         final FieldSinCos<T> sc = FastMath.sinCos(azimuth);

  220.         // Field
  221.         final Field<T> field = date.getField();
  222.         final T zero = field.getZero();
  223.         final T one  = field.getOne();

  224.         // degrees to semicircles
  225.         final T pi       = one.getPi();
  226.         final T rad2semi = pi.reciprocal();

  227.         // Earth Centered angle
  228.         final T psi = elevation.divide(pi).add(0.11).divide(0.0137).reciprocal().subtract(0.022);

  229.         // Subionospheric latitude: the latitude of the IPP (Ionospheric Pierce Point)
  230.         // in [-0.416, 0.416], semicircle
  231.         final T latIono = FastMath.min(
  232.                                       FastMath.max(geo.getLatitude().multiply(rad2semi).add(psi.multiply(sc.cos())), zero.subtract(0.416)),
  233.                                       zero.newInstance(0.416));

  234.         // Subionospheric longitude: the longitude of the IPP
  235.         // in semicircle
  236.         final T lonIono = geo.getLongitude().multiply(rad2semi).add(psi.multiply(sc.sin()).divide(FastMath.cos(latIono.multiply(pi))));

  237.         // Geomagnetic latitude, semicircle
  238.         final T latGeom = latIono.add(FastMath.cos(lonIono.subtract(1.617).multiply(pi)).multiply(0.064));

  239.         // day of week and tow (sec)
  240.         // Note: Sunday=0, Monday=1, Tuesday=2, Wednesday=3, Thursday=4, Friday=5, Saturday=6
  241.         final DateTimeComponents dtc = date.getComponents(gps);
  242.         final int dofweek = dtc.getDate().getDayOfWeek();
  243.         final double secday = dtc.getTime().getSecondsInLocalDay();
  244.         final double tow = dofweek * 86400. + secday;

  245.         final T t = lonIono.multiply(43200.).add(tow);
  246.         final T tsec = t.subtract(FastMath.floor(t.divide(86400.)).multiply(86400.)); // Seconds of day

  247.         // Slant factor, semicircle
  248.         final T slantFactor = FastMath.pow(elevation.divide(pi).negate().add(0.53), 3).multiply(16.0).add(one);

  249.         // Period of model, seconds
  250.         final T period = FastMath.max(zero.newInstance(72000.), latGeom.multiply(latGeom.multiply(latGeom.multiply(beta[3]).add(beta[2])).add(beta[1])).add(beta[0]));

  251.         // Phase of the model, radians
  252.         // (Max at 14.00 = 50400 sec local time)
  253.         final T x = tsec.subtract(50400.0).multiply(pi.multiply(2.0)).divide(period);

  254.         // Amplitude of the model, seconds
  255.         final T amplitude = FastMath.max(zero, latGeom.multiply(latGeom.multiply(latGeom.multiply(alpha[3]).add(alpha[2])).add(alpha[1])).add(alpha[0]));

  256.         // Ionospheric correction (L1)
  257.         T ionoTimeDelayL1 = slantFactor.multiply(5. * 1e-9);
  258.         if (FastMath.abs(x.getReal()) < 1.570) {
  259.             ionoTimeDelayL1 = ionoTimeDelayL1.add(slantFactor.multiply(amplitude.multiply(one.subtract(FastMath.pow(x, 2).multiply(0.5)).add(FastMath.pow(x, 4).divide(24.0)))));
  260.         }

  261.         // Ionospheric delay for the L1 frequency, in meters, with slant correction.
  262.         final double ratio = FastMath.pow(1575.42e6 / frequency, 2);
  263.         return ionoTimeDelayL1.multiply(Constants.SPEED_OF_LIGHT).multiply(ratio);
  264.     }

  265.     /** {@inheritDoc} */
  266.     @Override
  267.     public <T extends CalculusFieldElement<T>> T pathDelay(final FieldSpacecraftState<T> state, final TopocentricFrame baseFrame,
  268.                                                        final double frequency, final T[] parameters) {
  269.         return pathDelay(state, baseFrame, state.getDate(), frequency, parameters);
  270.     }

  271.     /** {@inheritDoc} */
  272.     @Override
  273.     public <T extends CalculusFieldElement<T>> T pathDelay(final FieldSpacecraftState<T> state,
  274.                                                            final TopocentricFrame baseFrame, final FieldAbsoluteDate<T> receptionDate,
  275.                                                            final double frequency, final T[] parameters) {

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

  281.         // Elevation in radians
  282.         final T elevation = position.getDelta();

  283.         if (elevation.getReal() > 0.0) {
  284.             // Date
  285.             final FieldAbsoluteDate<T> date = state.getDate();
  286.             // Geodetic point
  287.             final FieldGeodeticPoint<T> geo = baseFrame.getPoint(date.getField());
  288.             // Azimuth angle in radians
  289.             T azimuth = FastMath.atan2(position.getX(), position.getY());
  290.             if (azimuth.getReal() < 0.) {
  291.                 azimuth = azimuth.add(MathUtils.TWO_PI);
  292.             }
  293.             // Delay
  294.             return pathDelay(date, geo, elevation, azimuth, frequency, parameters);
  295.         }

  296.         return elevation.getField().getZero();
  297.     }

  298.     @Override
  299.     public List<ParameterDriver> getParametersDrivers() {
  300.         return Collections.emptyList();
  301.     }
  302. }