NeQuickGalileo.java

  1. /* Copyright 2002-2025 Thales Alenia Space
  2.  * Licensed to CS Communication & Systèmes (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.nequick;

  19. import org.hipparchus.CalculusFieldElement;
  20. import org.hipparchus.util.FastMath;
  21. import org.hipparchus.util.MathArrays;
  22. import org.orekit.annotation.DefaultDataContext;
  23. import org.orekit.bodies.FieldGeodeticPoint;
  24. import org.orekit.bodies.GeodeticPoint;
  25. import org.orekit.data.DataContext;
  26. import org.orekit.data.DataSource;
  27. import org.orekit.errors.OrekitException;
  28. import org.orekit.errors.OrekitMessages;
  29. import org.orekit.time.DateTimeComponents;
  30. import org.orekit.time.TimeScale;

  32. /** Galileo-specific version of NeQuick engine.
  33.  * @since 13.0
  34.  * @author Luc Maisonobe
  35.  */
  36. public class NeQuickGalileo extends NeQuickModel {

  38.     /** Starting number of points for integration. */
  39.     private static final int N_START = 8;

  41.     /** Broadcast ionization engine coefficients. */
  42.     private final double[] alpha;

  44.     /**
  45.      * Build a new instance.
  46.      *
  47.      * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
  48.      *
  49.      * @param alpha effective ionisation level coefficients
  50.      * @see #NeQuickGalileo(double[], TimeScale)
  51.      */
  52.     @DefaultDataContext
  53.     public NeQuickGalileo(final double[] alpha) {
  54.         this(alpha, DataContext.getDefault().getTimeScales().getUTC());
  55.     }

  57.     /**
  58.      * Build a new instance of the Galileo version of the NeQuick-2 model.
  59.      * <p>
  60.      * The Galileo version uses a loose modip grid and 3 broadcast parameters to compute
  61.      * effective ionization level.
  62.      * </p>
  63.      * @param alpha broadcast effective ionisation level coefficients
  64.      * @param utc UTC time scale.
  65.      * @since 10.1
  66.      */
  67.     public NeQuickGalileo(final double[] alpha, final TimeScale utc) {
  68.         super(utc);
  69.         this.alpha = alpha.clone();
  70.     }

  72.     /** Get effective ionisation level coefficients.
  73.      * @return effective ionisation level coefficients
  74.      */
  75.     public double[] getAlpha() {
  76.         return alpha.clone();
  77.     }

  79.     /**
  80.      * Compute effective ionization level.
  81.      *
  82.      * @param modip modified dip latitude at receiver location
  83.      * @return effective ionization level (Az in Nequick Galileo, R12 in original Nequick ITU)
  84.      */
  85.     private double effectiveIonizationLevel(final double modip) {
  86.         // Particular condition (Eq. 17)
  87.         if (alpha[0] == 0.0 && alpha[1] == 0.0 && alpha[2] == 0.0) {
  88.             return 63.7;
  89.         } else {
  90.             // Az = a0 + modip * a1 + modip² * a2 (Eq. 18)
  91.             return FastMath.min(FastMath.max(alpha[0] + modip * (alpha[1] + modip * alpha[2]), 0.0), 400.0);
  92.         }
  93.     }

  95.     /**
  96.      * Compute effective ionization level.
  97.      *
  98.      * @param <T>   type of the field elements
  99.      * @param modip modified dip latitude at receiver location
  100.      * @return effective ionization level (Az in Nequick Galileo, R12 in original Nequick ITU)
  101.      */
  102.     private <T extends CalculusFieldElement<T>> T effectiveIonizationLevel(final T modip) {
  103.         // Particular condition (Eq. 17)
  104.         if (alpha[0] == 0.0 && alpha[1] == 0.0 && alpha[2] == 0.0) {
  105.             return modip.newInstance(63.7);
  106.         } else {
  107.             // Az = a0 + modip * a1 + modip² * a2 (Eq. 18)
  108.             return FastMath.min(FastMath.max(modip.multiply(alpha[2]).add(alpha[1]).multiply(modip).add(alpha[0]),
  109.                                              0.0),
  110.                                 400.0);
  111.         }
  112.     }

  114.     /** {@inheritDoc} */
  115.     @Override
  116.     double stec(final DateTimeComponents dateTime, final Ray ray) {

  118.         // Tolerance for the integration accuracy. Defined inside the reference document, section 2.5.8.1.
  119.         final double h1 = ray.getRecH();
  120.         final double tolerance;
  121.         if (h1 < 1000000.0) {
  122.             tolerance = 0.001;
  123.         } else {
  124.             tolerance = 0.01;
  125.         }

  127.         // Integration
  128.         int n = N_START;
  129.         final Segment seg1 = new Segment(n, ray, ray.getS1(), ray.getS2());
  130.         double gn1 = stecIntegration(dateTime, seg1);
  131.         n *= 2;
  132.         final Segment seg2 = new Segment(n, ray, ray.getS1(), ray.getS2());
  133.         double gn2 = stecIntegration(dateTime, seg2);

  135.         int count = 1;
  136.         while (FastMath.abs(gn2 - gn1) > tolerance * FastMath.abs(gn1) && count < 20) {
  137.             gn1 = gn2;
  138.             n *= 2;
  139.             final Segment seg = new Segment(n, ray, ray.getS1(), ray.getS2());
  140.             gn2 = stecIntegration(dateTime, seg);
  141.             count += 1;
  142.         }

  144.         // If count > 20 the integration did not converge
  145.         if (count == 20) {
  146.             throw new OrekitException(OrekitMessages.STEC_INTEGRATION_DID_NOT_CONVERGE);
  147.         }

  149.         // Eq. 202
  150.         return (gn2 + ((gn2 - gn1) / 15.0)) * 1.0e-16;

  152.     }

  154.     /** {@inheritDoc} */
  155.     @Override
  156.     <T extends CalculusFieldElement<T>> T stec(final DateTimeComponents dateTime, final FieldRay<T> ray) {

  158.         // Tolerance for the integration accuracy. Defined inside the reference document, section 2.5.8.1.
  159.         final T h1 = ray.getRecH();
  160.         final double tolerance;
  161.         if (h1.getReal() < 1000000.0) {
  162.             tolerance = 0.001;
  163.         } else {
  164.             tolerance = 0.01;
  165.         }

  167.         // Integration
  168.         int n = N_START;
  169.         final FieldSegment<T> seg1 = new FieldSegment<>(n, ray, ray.getS1(), ray.getS2());
  170.         T gn1 = stecIntegration(dateTime, seg1);
  171.         n *= 2;
  172.         final FieldSegment<T> seg2 = new FieldSegment<>(n, ray, ray.getS1(), ray.getS2());
  173.         T gn2 = stecIntegration(dateTime, seg2);

  175.         int count = 1;
  176.         while (FastMath.abs(gn2.subtract(gn1)).getReal() > FastMath.abs(gn1).multiply(tolerance).getReal() &&
  177.                count < 20) {
  178.             gn1 = gn2;
  179.             n *= 2;
  180.             final FieldSegment<T> seg = new FieldSegment<>(n, ray, ray.getS1(), ray.getS2());
  181.             gn2 = stecIntegration(dateTime, seg);
  182.             count += 1;
  183.         }

  185.         // If count > 20 the integration did not converge
  186.         if (count == 20) {
  187.             throw new OrekitException(OrekitMessages.STEC_INTEGRATION_DID_NOT_CONVERGE);
  188.         }

  190.         // Eq. 202
  191.         return gn2.add(gn2.subtract(gn1).divide(15.0)).multiply(1.0e-16);

  193.     }

  195.     /**
  196.      * This method performs the STEC integration.
  197.      *
  198.      * @param dateTime current date and time components
  199.      * @param seg      coordinates along the integration path
  200.      * @return result of the integration
  201.      */
  202.     private double stecIntegration(final DateTimeComponents dateTime, final Segment seg) {

  204.         // Compute electron density
  205.         double density = 0.0;
  206.         for (int i = 0; i < seg.getNbPoints(); i++) {
  207.             final GeodeticPoint gp = seg.getPoint(i);
  208.             final double modip = GalileoHolder.INSTANCE.computeMODIP(gp.getLatitude(), gp.getLongitude());
  209.             density += electronDensity(computeFourierTimeSeries(dateTime, effectiveIonizationLevel(modip)),
  210.                                        gp.getLatitude(), gp.getLongitude(), gp.getAltitude());
  211.         }

  213.         return 0.5 * seg.getInterval() * density;
  214.     }

  216.     /**
  217.      * This method performs the STEC integration.
  218.      *
  219.      * @param <T>      type of the elements
  220.      * @param dateTime current date and time components
  221.      * @param seg      coordinates along the integration path
  222.      * @return result of the integration
  223.      */
  224.     private <T extends CalculusFieldElement<T>> T stecIntegration(final DateTimeComponents dateTime,
  225.                                                               final FieldSegment<T> seg) {
  226.         // Compute electron density
  227.         T density = seg.getInterval().getField().getZero();
  228.         for (int i = 0; i < seg.getNbPoints(); i++) {
  229.             final FieldGeodeticPoint<T> gp = seg.getPoint(i);
  230.             final T modip = GalileoHolder.INSTANCE.computeMODIP(gp.getLatitude(), gp.getLongitude());
  231.             density = density.add(electronDensity(computeFourierTimeSeries(dateTime, effectiveIonizationLevel(modip)),
  232.                                                   gp.getLatitude(), gp.getLongitude(), gp.getAltitude()));
  233.         }

  235.         return seg.getInterval().multiply(density).multiply(0.5);
  236.     }

  238.     /** {@inheritDoc} */
  239.     @Override
  240.     protected double computeMODIP(final double latitude, final double longitude) {
  241.         return GalileoHolder.INSTANCE.computeMODIP(latitude, longitude);
  242.     }

  244.     /** {@inheritDoc} */
  245.     @Override
  246.     protected <T extends CalculusFieldElement<T>> T computeMODIP(final T latitude, final T longitude) {
  247.         return GalileoHolder.INSTANCE.computeMODIP(latitude, longitude);
  248.     }

  250.     /** {@inheritDoc} */
  251.     @Override
  252.     boolean applyChapmanParameters(final double hInKm) {
  253.         return hInKm < 100.0;
  254.     }

  256.     /** {@inheritDoc} */
  257.     @Override
  258.     double[] computeExponentialArguments(final double h, final NeQuickParameters parameters) {

  260.         // If h < 100km we use h = 100km as recommended in the reference document
  261.         // for equations 111 to 113
  262.         final double clippedH = FastMath.max(h, 100.0);

  264.         // Eq. 111 to 113
  265.         final double   exp       = clipExp(10.0 / (1.0 + FastMath.abs(clippedH - parameters.getHmF2())));
  266.         final double[] arguments = new double[3];
  267.         arguments[0] =  (clippedH - parameters.getHmF2()) / parameters.getB2Bot();
  268.         arguments[1] = ((clippedH - parameters.getHmF1()) / parameters.getBF1(h)) * exp;
  269.         arguments[2] = ((clippedH - parameters.getHmE())  / parameters.getBE(h))  * exp;
  270.         return arguments;

  272.     }

  274.     /** {@inheritDoc} */
  275.     @Override
  276.     <T extends CalculusFieldElement<T>> T[] computeExponentialArguments(final T h,
  277.                                                                         final FieldNeQuickParameters<T> parameters) {
  278.         // If h < 100km we use h = 100km as recommended in the reference document
  279.         // for equations 111 to 113
  280.         final T clippedH = FastMath.max(h, 100);

  282.         // Eq. 111 to 113
  283.         final T   exp   = clipExp(FastMath.abs(clippedH.subtract(parameters.getHmF2())).add(1.0).reciprocal().multiply(10.0));
  284.         final T[] arguments = MathArrays.buildArray(h.getField(), 3);
  285.         arguments[0] = clippedH.subtract(parameters.getHmF2()).divide(parameters.getB2Bot());
  286.         arguments[1] = clippedH.subtract(parameters.getHmF1()).divide(parameters.getBF1(h)).multiply(exp);
  287.         arguments[2] = clippedH.subtract(parameters.getHmE()).divide(parameters.getBE(h)).multiply(exp);
  288.         return arguments;

  290.     }

  292.     /** {@inheritDoc} */
  293.     @Override
  294.     double computeH0(final NeQuickParameters parameters) {

  296.         final int month = parameters.getDateTime().getDate().getMonth();

  298.         // Auxiliary parameter ka (Eq. 99 and 100)
  299.         final double ka;
  300.         if (month > 3 && month < 10) {
  301.             // month = 4,5,6,7,8,9
  302.             ka = 6.705 - 0.014 * parameters.getAzr() - 0.008 * parameters.getHmF2();
  303.         } else {
  304.             // month = 1,2,3,10,11,12
  305.             final double ratio = parameters.getHmF2() / parameters.getB2Bot();
  306.             ka = -7.77 + 0.097 * ratio * ratio + 0.153 * parameters.getNmF2();
  307.         }

  309.         // Auxiliary parameter kb (Eq. 101 and 102)
  310.         double kb = parameters.join(ka, 2.0, 1.0, ka - 2.0);
  311.         kb = parameters.join(8.0, kb, 1.0, kb - 8.0);

  313.         // Auxiliary parameter Ha (Eq. 103)
  314.         final double hA = kb * parameters.getB2Bot();

  316.         // Auxiliary parameters x and v (Eq. 104 and 105)
  317.         final double x = 0.01 * (hA - 150.0);
  318.         final double v = (0.041163 * x - 0.183981) * x + 1.424472;

  320.         // Topside thickness parameter (Eq. 106)
  321.         return hA / v;

  323.     }

  325.     /** {@inheritDoc} */
  326.     @Override
  327.     <T extends CalculusFieldElement<T>> T computeH0(final FieldNeQuickParameters<T> parameters) {

  329.         final int month = parameters.getDateTime().getDate().getMonth();

  331.         // One
  332.         final T one = parameters.getAzr().getField().getOne();

  334.         // Auxiliary parameter ka (Eq. 99 and 100)
  335.         final T ka;
  336.         if (month > 3 && month < 10) {
  337.             // month = 4,5,6,7,8,9
  338.             ka = parameters.getAzr().multiply(0.014).add(parameters.getHmF2().multiply(0.008)).negate().add(6.705);
  339.         } else {
  340.             // month = 1,2,3,10,11,12
  341.             final T ratio = parameters.getHmF2().divide(parameters.getB2Bot());
  342.             ka = ratio.multiply(ratio).multiply(0.097).add(parameters.getNmF2().multiply(0.153)).add(-7.77);
  343.         }

  345.         // Auxiliary parameter kb (Eq. 101 and 102)
  346.         T kb = parameters.join(ka, one.newInstance(2.0), one, ka.subtract(2.0));
  347.         kb = parameters.join(one.newInstance(8.0), kb, one, kb.subtract(8.0));

  349.         // Auxiliary parameter Ha (Eq. 103)
  350.         final T hA = kb.multiply(parameters.getB2Bot());

  352.         // Auxiliary parameters x and v (Eq. 104 and 105)
  353.         final T x = hA.subtract(150.0).multiply(0.01);
  354.         final T v = x.multiply(0.041163).subtract(0.183981).multiply(x).add(1.424472);

  356.         // Topside thickness parameter (Eq. 106)
  357.         return hA.divide(v);

  359.     }

  361.     /** Holder for the Galileo-specific modip singleton.
  362.      * <p>
  363.      * We use the initialization on demand holder idiom to store the singleton,
  364.      * as it is both thread-safe, efficient (no synchronization) and works with
  365.      * all versions of java.
  366.      * </p>
  367.      */
  368.     private static class GalileoHolder {

  370.         /** Resource for modip grid. */
  371.         private static final String MODIP_GRID = "/assets/org/orekit/nequick/modipNeQG_wrapped.asc";

  373.         /** Unique instance. */
  374.         private static final ModipGrid INSTANCE =
  375.             new ModipGrid(36, 36,
  376.                           new DataSource(MODIP_GRID, () -> GalileoHolder.class.getResourceAsStream(MODIP_GRID)),
  377.                           true);

  379.         /** Private constructor.
  380.          * <p>This class is a utility class, it should neither have a public
  381.          * nor a default constructor. This private constructor prevents
  382.          * the compiler from generating one automatically.</p>
  383.          */
  384.         private GalileoHolder() {
  385.         }

  387.     }

  389. }
Created with JaCoCo 0.8.12.202403310830