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.             final double az = effectiveIonizationLevel(modip);
  210.             density += electronDensity(dateTime, az, 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.             final T az = effectiveIonizationLevel(modip);
  232.             density = density.add(electronDensity(dateTime, az,
  233.                                                   gp.getLatitude(), gp.getLongitude(), gp.getAltitude()));
  234.         }

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

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

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

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

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

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

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

  273.     }

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

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

  291.     }

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

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

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

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

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

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

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

  324.     }

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

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

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

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

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

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

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

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

  360.     }

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

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

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

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

  388.     }

  390. }
Created with JaCoCo 0.8.12.202403310830