NeQuickItu.java

  1. /* Copyright 2022-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.bodies.FieldGeodeticPoint;
  23. import org.orekit.bodies.GeodeticPoint;
  24. import org.orekit.data.DataSource;
  25. import org.orekit.time.DateTimeComponents;
  26. import org.orekit.time.TimeScale;
  27. import org.orekit.utils.units.Unit;

  29. /** Original model from Aeronomy and Radiopropagation Laboratory
  30.  * of the Abdus Salam International Centre for Theoretical Physics Trieste, Italy.
  31.  * <p>
  32.  * None of the code from Abdus Salam International Centre for Theoretical Physics Trieste
  33.  * has been used, the models have been reimplemented from scratch by the Orekit team.
  34.  * </p>
  35.  * @since 13.0
  36.  * @author Luc Maisonobe
  37.  */
  38. public class NeQuickItu extends NeQuickModel {

  40.     /** One thousand kilometer height. */
  41.     private static final double H_1000 = 1000000.0;

  43.     /** Two thousands kilometer height. */
  44.     private static final double H_2000 = 2000000.0;

  46.     /** H0 (km). */
  47.     private static final double H0 = 90.0;

  49.     /** HD (km). */
  50.     private static final double HD = 5.0;

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

  55.     /** Max number of points for integration. */
  56.     private static final int N_STOP = 1024;

  58.     /** Small convergence criterion. */
  59.     private static final double EPS_SMALL = 1.0e-3;

  61.     /** Medium convergence criterion. */
  62.     private static final double EPS_MEDIUM = 1.0e-2;

  64.     /** Solar flux. */
  65.     private final double f107;

  67.     /** Build a new instance.
  68.      * @param f107 solar flux
  69.      * @param utc UTC time scale
  70.      */
  71.     public NeQuickItu(final double f107, final TimeScale utc) {
  72.         super(utc);
  73.         this.f107 = f107;
  74.     }

  76.     /** Get solar flux.
  77.      * @return solar flux
  78.      */
  79.     public double getF107() {
  80.         return f107;
  81.     }

  83.     /** {@inheritDoc} */
  84.     @Override
  85.     double stec(final DateTimeComponents dateTime, final Ray ray) {
  86.         if (ray.getSatH() <= H_2000) {
  87.             if (ray.getRecH() >= H_1000) {
  88.                 // only one integration interval
  89.                 return stecIntegration(dateTime, EPS_SMALL, ray, ray.getS1(), ray.getS2());
  90.             } else {
  91.                 // two integration intervals, below and above 1000km
  92.                 final double h1000 = NeQuickModel.RE + H_1000;
  93.                 final double s1000 = FastMath.sqrt(h1000 * h1000 - ray.getRadius() * ray.getRadius());
  94.                 return stecIntegration(dateTime, EPS_SMALL, ray, ray.getS1(), s1000) +
  95.                        stecIntegration(dateTime, EPS_MEDIUM, ray, s1000, ray.getS2());
  96.             }
  97.         } else {
  98.             if (ray.getRecH() >= H_2000) {
  99.                 // only one integration interval
  100.                 return stecIntegration(dateTime, EPS_SMALL, ray, ray.getS1(), ray.getS2());
  101.             } else {
  102.                 final double h2000 = NeQuickModel.RE + H_2000;
  103.                 final double s2000 = FastMath.sqrt(h2000 * h2000 - ray.getRadius() * ray.getRadius());
  104.                 if (ray.getRecH() >= H_1000) {
  105.                     // two integration intervals, below and above 2000km
  106.                     return stecIntegration(dateTime, EPS_SMALL, ray, ray.getS1(), s2000) +
  107.                            stecIntegration(dateTime, EPS_SMALL, ray, s2000, ray.getS2());
  108.                 } else {
  109.                     // three integration intervals, below 1000km, between 1000km and 2000km, and above 2000km
  110.                     final double h1000 = NeQuickModel.RE + H_1000;
  111.                     final double s1000 = FastMath.sqrt(h1000 * h1000 - ray.getRadius() * ray.getRadius());
  112.                     return stecIntegration(dateTime, EPS_SMALL, ray, ray.getS1(), s1000) +
  113.                            stecIntegration(dateTime, EPS_MEDIUM, ray, s1000, s2000) +
  114.                            stecIntegration(dateTime, EPS_MEDIUM, ray, s2000, ray.getS2());
  115.                 }
  116.             }
  117.         }
  118.     }

  120.     /** {@inheritDoc} */
  121.     @Override
  122.     <T extends CalculusFieldElement<T>> T stec(final DateTimeComponents dateTime, final FieldRay<T> ray) {
  123.         if (ray.getSatH().getReal() <= H_2000) {
  124.             if (ray.getRecH().getReal() >= H_1000) {
  125.                 // only one integration interval
  126.                 return stecIntegration(dateTime, EPS_SMALL, ray, ray.getS1(), ray.getS2());
  127.             } else {
  128.                 // two integration intervals, below and above 1000km
  129.                 final double h1000 = NeQuickModel.RE + H_1000;
  130.                 final T s1000 = FastMath.sqrt(ray.getRadius().multiply(ray.getRadius()).negate().add(h1000 * h1000));
  131.                 return stecIntegration(dateTime, EPS_SMALL, ray, ray.getS1(), s1000).
  132.                        add(stecIntegration(dateTime, EPS_MEDIUM, ray, s1000, ray.getS2()));
  133.             }
  134.         } else {
  135.             if (ray.getRecH().getReal() >= H_2000) {
  136.                 // only one integration interval
  137.                 return stecIntegration(dateTime, EPS_SMALL, ray, ray.getS1(), ray.getS2());
  138.             } else {
  139.                 final double h2000 = NeQuickModel.RE + H_2000;
  140.                 final T s2000 = FastMath.sqrt(ray.getRadius().multiply(ray.getRadius()).negate().add(h2000 * h2000));
  141.                 if (ray.getRecH().getReal() >= H_1000) {
  142.                     // two integration intervals, below and above 2000km
  143.                     return stecIntegration(dateTime, EPS_SMALL, ray, ray.getS1(), s2000).
  144.                            add(stecIntegration(dateTime, EPS_SMALL, ray, s2000, ray.getS2()));
  145.                 } else {
  146.                     // three integration intervals, below 1000km, between 1000km and 2000km, and above 2000km
  147.                     final double h1000 = NeQuickModel.RE + H_1000;
  148.                     final T s1000 = FastMath.sqrt(ray.getRadius().multiply(ray.getRadius()).negate().add(h1000 * h1000));
  149.                     return stecIntegration(dateTime, EPS_SMALL, ray, ray.getS1(), s1000).
  150.                            add(stecIntegration(dateTime, EPS_MEDIUM, ray, s1000, s2000)).
  151.                            add(stecIntegration(dateTime, EPS_MEDIUM, ray, s2000, ray.getS2()));
  152.                 }
  153.             }
  154.         }
  155.     }

  157.     /**
  158.      * This method performs the STEC integration.
  159.      *
  160.      * @param dateTime current date and time components
  161.      * @param eps convergence criterion
  162.      * @param ray ray-perigee parameters
  163.      * @param s1  lower boundary of integration
  164.      * @param s2  upper boundary for integration
  165.      * @return result of the integration
  166.      */
  167.     private double stecIntegration(final DateTimeComponents dateTime, final double eps, final Ray ray, final double s1,
  168.                                    final double s2) {

  170.         double gInt1 = Double.NaN;
  171.         double gInt2 = Double.NaN;

  173.         for (int n = N_START; n <= N_STOP; n = 2 * n) {

  175.             // integrate with n intervals (2n points)
  176.             final Segment segment = new Segment(n, ray, s1, s2);
  177.             double sum = 0;
  178.             for (int i = 0; i < segment.getNbPoints(); ++i) {
  179.                 final GeodeticPoint gp = segment.getPoint(i);
  180.                 final double ed = electronDensity(dateTime, f107,
  181.                                                   gp.getLatitude(), gp.getLongitude(), gp.getAltitude());
  182.                 sum += ed;
  183.             }

  185.             gInt1 = gInt2;
  186.             gInt2 = sum * 0.5 * segment.getInterval();
  187.             if (FastMath.abs(gInt1 - gInt2) <= FastMath.abs(gInt1 * eps)) {
  188.                 // convergence reached
  189.                 break;
  190.             }

  192.         }

  194.         return Unit.TOTAL_ELECTRON_CONTENT_UNIT.fromSI(gInt2 + (gInt2 - gInt1) / 15.0);

  196.     }

  198.     /**
  199.      * This method performs the STEC integration.
  200.      *
  201.      * @param <T> type of the field elements
  202.      * @param dateTime current date and time components
  203.      * @param eps convergence criterion
  204.      * @param ray ray-perigee parameters
  205.      * @param s1  lower boundary of integration
  206.      * @param s2  upper boundary for integration
  207.      * @return result of the integration
  208.      */
  209.     private <T extends CalculusFieldElement<T>> T stecIntegration(final DateTimeComponents dateTime,
  210.                                                                   final double eps,
  211.                                                                   final FieldRay<T> ray, final T s1, final T s2) {

  213.         T gInt1 = s1.newInstance(Double.NaN);
  214.         T gInt2 = s1.newInstance(Double.NaN);
  215.         final T f107T = s1.newInstance(f107);

  217.         for (int n = N_START; n <= N_STOP; n = 2 * n) {

  219.             // integrate with n intervals (2n points)
  220.             final FieldSegment<T> segment = new FieldSegment<>(n, ray, s1, s2);
  221.             T sum = s1.getField().getZero();
  222.             for (int i = 0; i < segment.getNbPoints(); ++i) {
  223.                 final FieldGeodeticPoint<T> gp = segment.getPoint(i);
  224.                 final T ed =  electronDensity(dateTime, f107T,
  225.                                               gp.getLatitude(), gp.getLongitude(), gp.getAltitude());
  226.                 sum = sum.add(ed);
  227.             }

  229.             gInt1 = gInt2;
  230.             gInt2 = sum.multiply(0.5).multiply(segment.getInterval());
  231.             if (FastMath.abs(gInt1.subtract(gInt2).getReal()) <= FastMath.abs(gInt1.getReal() * eps)) {
  232.                 // convergence reached
  233.                 break;
  234.             }

  236.         }

  238.         return Unit.TOTAL_ELECTRON_CONTENT_UNIT.fromSI(gInt2.add(gInt2.subtract(gInt1).divide(15.0)));

  240.     }

  242.     /** {@inheritDoc} */
  243.     @Override
  244.     protected double computeMODIP(final double latitude, final double longitude) {
  245.         return ItuHolder.INSTANCE.computeMODIP(latitude, longitude);
  246.     }

  248.     /** {@inheritDoc} */
  249.     @Override
  250.     protected <T extends CalculusFieldElement<T>> T computeMODIP(final T latitude, final T longitude) {
  251.         return ItuHolder.INSTANCE.computeMODIP(latitude, longitude);
  252.     }

  254.     /** {@inheritDoc} */
  255.     @Override
  256.     boolean applyChapmanParameters(final double hInKm) {
  257.         return false;
  258.     }

  260.     /** {@inheritDoc} */
  261.     @Override
  262.     double[] computeExponentialArguments(final double h, final NeQuickParameters parameters) {
  263.         final double   exp   = clipExp(10.0 / (1.0 + FastMath.abs(h - parameters.getHmF2())));
  264.         final double[] arguments = new double[3];
  265.         arguments[0] = fixLowHArg( (h - parameters.getHmF2()) / parameters.getB2Bot(), h);
  266.         arguments[1] = fixLowHArg(((h - parameters.getHmF1()) / parameters.getBF1(h)) * exp, h);
  267.         arguments[2] = fixLowHArg(((h - parameters.getHmE())  / parameters.getBE(h))  * exp, h);
  268.         return arguments;
  269.     }

  271.     /** {@inheritDoc} */
  272.     @Override
  273.     <T extends CalculusFieldElement<T>> T[] computeExponentialArguments(final T h,
  274.                                                                         final FieldNeQuickParameters<T> parameters) {
  275.         final T   exp   = clipExp(FastMath.abs(h.subtract(parameters.getHmF2())).add(1.0).reciprocal().multiply(10.0));
  276.         final T[] arguments = MathArrays.buildArray(h.getField(), 3);
  277.         arguments[0] = fixLowHArg(h.subtract(parameters.getHmF2()).divide(parameters.getB2Bot()), h);
  278.         arguments[1] = fixLowHArg(h.subtract(parameters.getHmF1()).divide(parameters.getBF1(h)).multiply(exp), h);
  279.         arguments[2] = fixLowHArg(h.subtract(parameters.getHmE()).divide(parameters.getBE(h)).multiply(exp), h);
  280.         return arguments;
  281.     }

  283.     /**
  284.      * Fix arguments for low altitudes.
  285.      * @param arg argument of the exponential
  286.      * @param h height in km
  287.      * @return fixed argument
  288.      * @since 13.0
  289.      */
  290.     private double fixLowHArg(final double arg, final double h) {
  291.         return h < H0 ? arg * (HD + H0 - h) / HD : arg;
  292.     }

  294.     /**
  295.      * Fix arguments for low altitudes.
  296.      * @param <T> type of the field elements
  297.      * @param arg argument of the exponential
  298.      * @param h height in km
  299.      * @return fixed argument
  300.      * @since 13.0
  301.      */
  302.     private <T extends CalculusFieldElement<T>> T fixLowHArg(final T arg, final T h) {
  303.         return h.getReal() < H0 ? arg.multiply(h.negate().add(HD + H0)).divide(HD) : arg;
  304.     }

  306.     /** {@inheritDoc} */
  307.     @Override
  308.     double computeH0(final NeQuickParameters parameters) {
  309.         final double b2k = -0.0538  * parameters.getFoF2() -
  310.                             0.00664 * parameters.getHmF2() +
  311.                             0.113   * parameters.getHmF2() / parameters.getB2Bot() +
  312.                             0.00257 * parameters.getAzr() +
  313.                             3.22;
  314.         return parameters.getB2Bot() * parameters.join(b2k, 1.0, 2.0, b2k - 1.0);
  315.     }

  317.     /** {@inheritDoc} */
  318.     @Override
  319.     <T extends CalculusFieldElement<T>> T computeH0(final FieldNeQuickParameters<T> parameters) {
  320.         final T b2k = parameters.getFoF2().multiply(-0.0538).
  321.                       subtract(parameters.getHmF2().multiply(0.00664)).
  322.                       add(parameters.getHmF2().multiply(0.113).divide(parameters.getB2Bot())).
  323.                       add(parameters.getAzr().multiply(0.00257)).
  324.                       add (3.22);
  325.         return parameters.getB2Bot().
  326.                multiply(parameters.join(b2k, b2k.newInstance(1.0), b2k.newInstance(2.0), b2k.subtract(1.0)));
  327.     }

  329.     /** Holder for the ITU modip singleton.
  330.      * <p>
  331.      * We use the initialization on demand holder idiom to store the singleton,
  332.      * as it is both thread-safe, efficient (no synchronization) and works with
  333.      * all versions of java.
  334.      * </p>
  335.      */
  336.     private static class ItuHolder {

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

  341.         /** Unique instance. */
  342.         private static final ModipGrid INSTANCE =
  343.             new ModipGrid(180, 180,
  344.                           new DataSource(MODIP_GRID, () -> ItuHolder.class.getResourceAsStream(MODIP_GRID)),
  345.                           false);

  347.         /** Private constructor.
  348.          * <p>This class is a utility class, it should neither have a public
  349.          * nor a default constructor. This private constructor prevents
  350.          * the compiler from generating one automatically.</p>
  351.          */
  352.         private ItuHolder() {
  353.         }

  355.     }

  357. }
Created with JaCoCo 0.8.12.202403310830