FieldNeQuickParameters.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.nequick;

  19. import org.hipparchus.CalculusFieldElement;
  20. import org.hipparchus.util.FastMath;
  21. import org.hipparchus.util.FieldSinCos;
  22. import org.hipparchus.util.MathArrays;
  23. import org.orekit.time.DateComponents;
  24. import org.orekit.time.DateTimeComponents;
  25. import org.orekit.time.TimeComponents;

  27. /**
  28.  * This class performs the computation of the parameters used by the NeQuick model.
  29.  *
  30.  * @author Bryan Cazabonne
  31.  *
  32.  * @see "European Union (2016). European GNSS (Galileo) Open Service-Ionospheric Correction
  33.  *       Algorithm for Galileo Single Frequency Users. 1.2."
  34.  * @see <a href="https://www.itu.int/rec/R-REC-P.531/en">ITU-R P.531</a>
  35.  *
  36.  * @since 10.1
  37.  * @param <T> type of the field elements
  38.  */
  39. class FieldNeQuickParameters <T extends CalculusFieldElement<T>> {

  41.     /** Solar zenith angle at day night transition, degrees. */
  42.     private static final double X0 = 86.23292796211615;

  44.     /** Current date time components.
  45.      * @since 13.0
  46.      */
  47.     private final DateTimeComponents dateTime;

  49.     /** Effective sunspot number.
  50.      * @since 13.0
  51.      */
  52.     private final T azr;

  54.     /** F2 layer critical frequency.
  55.      * @since 13.0
  56.      */
  57.     private final T foF2;

  59.     /** F2 layer maximum density. */
  60.     private final T nmF2;

  62.     /** F2 layer maximum density height [km]. */
  63.     private final T hmF2;

  65.     /** F1 layer maximum density height [km]. */
  66.     private final T hmF1;

  68.     /** E layer maximum density height [km]. */
  69.     private final T hmE;

  71.     /** F2 layer bottom thickness parameter [km]. */
  72.     private final T b2Bot;

  74.     /** F1 layer top thickness parameter [km]. */
  75.     private final T b1Top;

  77.     /** F1 layer bottom thickness parameter [km]. */
  78.     private final T b1Bot;

  80.     /** E layer top thickness parameter [km]. */
  81.     private final T beTop;

  83.     /** E layer bottom thickness parameter [km]. */
  84.     private final T beBot;

  86.     /** Layer amplitudes. */
  87.     private final T[] amplitudes;

  89.     /**
  90.      * Build a new instance.
  91.      * @param dateTime current date time components
  92.      * @param flattenF2 F2 coefficients used by the F2 layer (flatten array)
  93.      * @param flattenFm3 Fm3 coefficients used by the F2 layer (flatten array)
  94.      * @param latitude latitude of a point along the integration path, in radians
  95.      * @param longitude longitude of a point along the integration path, in radians
  96.      * @param az effective ionisation level
  97.      * @param modip modip
  98.      */
  99.     FieldNeQuickParameters(final DateTimeComponents dateTime, final double[] flattenF2,
  100.                            final double[] flattenFm3, final T latitude, final T longitude, final T az,
  101.                            final T modip) {
  102.         this(new FieldFourierTimeSeries<>(dateTime, az, flattenF2, flattenFm3),
  103.              latitude, longitude, modip);
  104.     }

  106.     /**
  107.      * Build a new instance.
  108.      * @param fourierTimeSeries Fourier time series for foF2 and M(3000)F2 layer
  109.      * @param latitude latitude of a point along the integration path, in radians
  110.      * @param longitude longitude of a point along the integration path, in radians
  111.      * @param modip modip
  112.      */
  113.     FieldNeQuickParameters(final FieldFourierTimeSeries<T> fourierTimeSeries,
  114.                            final T latitude, final T longitude, final T modip) {

  116.         // Zero
  117.         final T zero = latitude.getField().getZero();

  119.         this.dateTime = fourierTimeSeries.getDateTime();
  120.         this.azr      = fourierTimeSeries.getAzr();

  122.         // Date and Time components
  123.         final DateComponents date = dateTime.getDate();
  124.         final TimeComponents time = dateTime.getTime();
  125.         // Hours
  126.         final double hours  = time.getSecondsInUTCDay() / 3600.0;
  127.         // Effective solar zenith angle in radians
  128.         final T xeff = computeEffectiveSolarAngle(date.getMonth(), hours, latitude, longitude);

  130.         // E layer maximum density height in km (Eq. 78)
  131.         this.hmE = zero.newInstance(120.0);
  132.         // E layer critical frequency in MHz
  133.         final T foE = computefoE(date.getMonth(), fourierTimeSeries.getAz(), xeff, latitude);
  134.         // E layer maximum density in 10^11 m-3 (Eq. 36)
  135.         final T nmE = foE.multiply(foE).multiply(0.124);

  137.         // F2 layer critical frequency in MHz
  138.         final T[] scL = FieldFourierTimeSeries.sinCos(longitude, 8);
  139.         this.foF2 = computefoF2(modip, fourierTimeSeries.getCf2Reference(), latitude, scL);
  140.         // Maximum Usable Frequency factor
  141.         final T mF2  = computeMF2(modip, fourierTimeSeries.getCm3Reference(), latitude, scL);
  142.         // F2 layer maximum density in 10^11 m-3
  143.         this.nmF2 = foF2.multiply(foF2).multiply(0.124);
  144.         // F2 layer maximum density height in km
  145.         this.hmF2 = computehmF2(foE, mF2);

  147.         // F1 layer critical frequency in MHz
  148.         final T foF1 = computefoF1(foE);
  149.         // F1 layer maximum density in 10^11 m-3
  150.         final T nmF1;
  151.         if (foF1.getReal() <= 0.0 && foE.getReal() > 2.0) {
  152.             final T foEpopf = foE.add(0.5);
  153.             nmF1 = foEpopf.multiply(foEpopf).multiply(0.124);
  154.         } else {
  155.             nmF1 = foF1.multiply(foF1).multiply(0.124);
  156.         }
  157.         // F1 layer maximum density height in km
  158.         this.hmF1 = hmF2.add(hmE).multiply(0.5);

  160.         // Thickness parameters (Eq. 85 to 89)
  161.         final T a = clipExp(FastMath.log(foF2.multiply(foF2)).multiply(0.857).add(FastMath.log(mF2).multiply(2.02)).add(-3.467)).multiply(0.01);
  162.         this.b2Bot = nmF2.divide(a).multiply(0.385);
  163.         this.b1Top = hmF2.subtract(hmF1).multiply(0.3);
  164.         this.b1Bot = hmF1.subtract(hmE).multiply(0.5);
  165.         this.beTop = FastMath.max(b1Bot, zero.newInstance(7.0));
  166.         this.beBot = zero.newInstance(5.0);

  168.         // Layer amplitude coefficients
  169.         this.amplitudes = computeLayerAmplitudes(nmE, nmF1, foF1);

  171.     }

  173.     /**
  174.      * Get current date time components.
  175.      * @return current date time components
  176.      * @since 13.0
  177.      */
  178.     public DateTimeComponents getDateTime() {
  179.         return dateTime;
  180.     }

  182.     /**
  183.      * Get effective sunspot number.
  184.      * @return effective sunspot number
  185.      * @since 13.0
  186.      */
  187.     public T getAzr() {
  188.         return azr;
  189.     }

  191.     /**
  192.      * Get F2 layer critical frequency.
  193.      * @return F2 layer critical frequency
  194.      * @since 13.0
  195.      */
  196.     public T getFoF2() {
  197.         return foF2;
  198.     }

  200.     /**
  201.      * Get the F2 layer maximum density.
  202.      * @return nmF2
  203.      */
  204.     public T getNmF2() {
  205.         return nmF2;
  206.     }

  208.     /**
  209.      * Get the F2 layer maximum density height.
  210.      * @return hmF2 in km
  211.      */
  212.     public T getHmF2() {
  213.         return hmF2;
  214.     }

  216.     /**
  217.      * Get the F1 layer maximum density height.
  218.      * @return hmF1 in km
  219.      */
  220.     public T getHmF1() {
  221.         return hmF1;
  222.     }

  224.     /**
  225.      * Get the E layer maximum density height.
  226.      * @return hmE in km
  227.      */
  228.     public T getHmE() {
  229.         return hmE;
  230.     }

  232.     /**
  233.      * Get the F2 layer thickness parameter (bottom).
  234.      * @return B2Bot in km
  235.      */
  236.     public T getB2Bot() {
  237.         return b2Bot;
  238.     }

  240.     /**
  241.      * Get the F1 layer thickness parameter.
  242.      * @param h current height (km)
  243.      * @return B1 in km
  244.      * @since 13.0
  245.      */
  246.     public T getBF1(final T h) {
  247.         // Eq. 110
  248.         return (h.getReal() > hmF1.getReal()) ? b1Top : b1Bot;
  249.     }

  251.     /**
  252.      * Get the E layer thickness parameter.
  253.      * @param h current height (km)
  254.      * @return Be in km
  255.      * @since 13.0
  256.      */
  257.     public T getBE(final T h) {
  258.         // Eq. 109
  259.         return (h.getReal() > hmE.getReal()) ? beTop : beBot;
  260.     }

  262.     /**
  263.      * Get the F2, F1 and E layer amplitudes.
  264.      * <p>
  265.      * The resulting element is an array having the following form:
  266.      * <ul>
  267.      * <li>double[0] = A1 → F2 layer amplitude
  268.      * <li>double[1] = A2 → F1 layer amplitude
  269.      * <li>double[2] = A3 → E  layer amplitude
  270.      * </ul>
  271.      * @return layer amplitudes
  272.      */
  273.     public T[] getLayerAmplitudes() {
  274.         return amplitudes.clone();
  275.     }

  277.     /**
  278.      * This method computes the effective solar zenith angle.
  279.      * <p>
  280.      * The effective solar zenith angle is compute as a function of the
  281.      * solar zenith angle and the solar zenith angle at day night transition.
  282.      * </p>
  283.      * @param month current month of the year
  284.      * @param hours universal time (hours)
  285.      * @param latitude in radians
  286.      * @param longitude in radians
  287.      * @return the effective solar zenith angle, radians
  288.      */
  289.     private T computeEffectiveSolarAngle(final int month,
  290.                                          final double hours,
  291.                                          final T latitude,
  292.                                          final T longitude) {
  293.         // Zero
  294.         final T zero = latitude.getField().getZero();
  295.         // Local time (Eq.4)
  296.         final T lt = longitude.divide(FastMath.toRadians(15.0)).add(hours);
  297.         // Day of year at the middle of the month (Eq. 20)
  298.         final double dy = 30.5 * month - 15.0;
  299.         // Time (Eq. 21)
  300.         final double t = dy + (18 - hours) / 24;
  301.         // Arguments am and al (Eq. 22 and 23)
  302.         final double am = FastMath.toRadians(0.9856 * t - 3.289);
  303.         final double al = am + FastMath.toRadians(1.916 * FastMath.sin(am) + 0.020 * FastMath.sin(2.0 * am) + 282.634);
  304.         // Sine and cosine of solar declination (Eq. 24 and 25)
  305.         final double sDec = 0.39782 * FastMath.sin(al);
  306.         final double cDec = FastMath.sqrt(1. - sDec * sDec);
  307.         // Solar zenith angle, deg (Eq. 26 and 27)
  308.         final FieldSinCos<T> scLat   = FastMath.sinCos(latitude);
  309.         final T coef    = lt.negate().add(12.0).multiply(FastMath.PI / 12);
  310.         final T cZenith = scLat.sin().multiply(sDec).add(scLat.cos().multiply(cDec).multiply(FastMath.cos(coef)));
  311.         final T angle   = FastMath.atan2(FastMath.sqrt(cZenith.multiply(cZenith).negate().add(1.0)), cZenith);
  312.         final T x       = FastMath.toDegrees(angle);
  313.         // Effective solar zenith angle (Eq. 28)
  314.         final T xeff = join(clipExp(x.multiply(0.2).negate().add(20.0)).multiply(0.24).negate().add(90.0), x,
  315.                 zero.newInstance(12.0), x.subtract(X0));
  316.         return FastMath.toRadians(xeff);
  317.     }

  319.     /**
  320.      * This method computes the E layer critical frequency at a given location.
  321.      * @param month current month
  322.      * @param az ffective ionisation level
  323.      * @param xeff effective solar zenith angle in radians
  324.      * @param latitude latitude in radians
  325.      * @return the E layer critical frequency at a given location in MHz
  326.      */
  327.     private T computefoE(final int month, final T az,
  328.                          final T xeff, final T latitude) {
  329.         // The latitude has to be converted in degrees
  330.         final T lat = FastMath.toDegrees(latitude);
  331.         // Square root of the effective ionisation level
  332.         final T sqAz = FastMath.sqrt(az);
  333.         // seas parameter (Eq. 30 to 32)
  334.         final int seas;
  335.         if (month == 1 || month == 2 || month == 11 || month == 12) {
  336.             seas = -1;
  337.         } else if (month == 3 || month == 4 || month == 9 || month == 10) {
  338.             seas = 0;
  339.         } else {
  340.             seas = 1;
  341.         }
  342.         // Latitudinal dependence (Eq. 33 and 34)
  343.         final T ee = clipExp(lat.multiply(0.3));
  344.         final T seasp = ee.subtract(1.0).divide(ee.add(1.0)).multiply(seas);
  345.         // Critical frequency (Eq. 35)
  346.         final T coef = seasp.multiply(0.019).negate().add(1.112);
  347.         return FastMath.sqrt(coef .multiply(coef).multiply(sqAz).multiply(FastMath.cos(xeff).pow(0.6)).add(0.49));

  349.     }

  351.     /**
  352.      * Computes the F2 layer height of maximum electron density.
  353.      * @param foE E layer layer critical frequency in MHz
  354.      * @param mF2 maximum usable frequency factor
  355.      * @return hmF2 in km
  356.      */
  357.     private T computehmF2(final T foE, final T mF2) {
  358.         // Zero
  359.         final T zero = foE.getField().getZero();
  360.         // Ratio
  361.         final T fo = foF2.divide(foE);
  362.         final T ratio = join(fo, zero.newInstance(1.75), zero.newInstance(20.0), fo.subtract(1.75));

  364.         // deltaM parameter
  365.         T deltaM = zero.subtract(0.012);
  366.         if (foE.getReal() >= 1e-30) {
  367.             deltaM = deltaM.add(ratio.subtract(1.215).divide(0.253).reciprocal());
  368.         }

  370.         // hmF2 Eq. 80
  371.         final T mF2Sq = mF2.square();
  372.         final T temp  = FastMath.sqrt(mF2Sq.multiply(0.0196).add(1.0).divide(mF2Sq.multiply(1.2967).subtract(1.0)));
  373.         return mF2.multiply(1490.0).multiply(temp).divide(mF2.add(deltaM)).subtract(176.0);

  375.     }

  377.     /**
  378.      * This method computes the F2 layer critical frequency.
  379.      * @param modip modified DIP latitude, in degrees
  380.      * @param cf2 Fourier time series for foF2
  381.      * @param latitude latitude in radians
  382.      * @param scL sines/cosines array of longitude argument
  383.      * @return the F2 layer critical frequency, MHz
  384.      */
  385.     private T computefoF2(final T modip, final T[] cf2,
  386.                           final T latitude, final T[] scL) {

  388.         // Legendre grades (Eq. 63)
  389.         final int[] q = new int[] {
  390.             12, 12, 9, 5, 2, 1, 1, 1, 1
  391.         };

  393.         T frequency = cf2[0];

  395.         // ModipGrid coefficients Eq. 57
  396.         final T sinMODIP = FastMath.sin(FastMath.toRadians(modip));
  397.         final T[] m = MathArrays.buildArray(latitude.getField(), 12);
  398.         m[0] = latitude.getField().getOne();
  399.         for (int i = 1; i < q[0]; i++) {
  400.             m[i] = sinMODIP.multiply(m[i - 1]);
  401.             frequency = frequency.add(m[i].multiply(cf2[i]));
  402.         }

  404.         // latitude and longitude terms
  405.         int index = 12;
  406.         final T cosLat1 = FastMath.cos(latitude);
  407.         T cosLatI = cosLat1;
  408.         for (int i = 1; i < q.length; i++) {
  409.             final T c = cosLatI.multiply(scL[2 * i - 1]);
  410.             final T s = cosLatI.multiply(scL[2 * i - 2]);
  411.             for (int j = 0; j < q[i]; j++) {
  412.                 frequency = frequency.add(m[j].multiply(c).multiply(cf2[index++]));
  413.                 frequency = frequency.add(m[j].multiply(s).multiply(cf2[index++]));
  414.             }
  415.             cosLatI = cosLatI.multiply(cosLat1);
  416.         }

  418.         return frequency;

  420.     }

  422.     /**
  423.      * This method computes the Maximum Usable Frequency factor.
  424.      * @param modip modified DIP latitude, in degrees
  425.      * @param cm3 Fourier time series for M(3000)F2
  426.      * @param latitude latitude in radians
  427.      * @param scL sines/cosines array of longitude argument
  428.      * @return the Maximum Usable Frequency factor
  429.      */
  430.     private T computeMF2(final T modip, final T[] cm3,
  431.                          final T latitude, final T[] scL) {

  433.         // Legendre grades (Eq. 71)
  434.         final int[] r = new int[] {
  435.             7, 8, 6, 3, 2, 1, 1
  436.         };

  438.         T m3000 = cm3[0];

  440.         // ModipGrid coefficients Eq. 57
  441.         final T sinMODIP = FastMath.sin(FastMath.toRadians(modip));
  442.         final T[] m = MathArrays.buildArray(latitude.getField(), 12);
  443.         m[0] = latitude.getField().getOne();
  444.         for (int i = 1; i < 12; i++) {
  445.             m[i] = sinMODIP.multiply(m[i - 1]);
  446.             if (i < 7) {
  447.                 m3000 = m3000.add(m[i].multiply(cm3[i]));
  448.             }
  449.         }

  451.         // latitude and longitude terms
  452.         int index = 7;
  453.         final T cosLat1 = FastMath.cos(latitude);
  454.         T cosLatI = cosLat1;
  455.         for (int i = 1; i < r.length; i++) {
  456.             final T c = cosLatI.multiply(scL[2 * i - 1]);
  457.             final T s = cosLatI.multiply(scL[2 * i - 2]);
  458.             for (int j = 0; j < r[i]; j++) {
  459.                 m3000 = m3000.add(m[j].multiply(c).multiply(cm3[index++]));
  460.                 m3000 = m3000.add(m[j].multiply(s).multiply(cm3[index++]));
  461.             }
  462.             cosLatI = cosLatI.multiply(cosLat1); // Eq. 58
  463.         }

  465.         return m3000;

  467.     }

  469.     /**
  470.      * This method computes the F1 layer critical frequency.
  471.      * <p>
  472.      * This computation performs the algorithm exposed in Annex F
  473.      * of the reference document.
  474.      * </p>
  475.      * @param foE the E layer critical frequency, MHz
  476.      * @return the F1 layer critical frequency, MHz
  477.      */
  478.     private T computefoF1(final T foE) {
  479.         final T zero = foE.getField().getZero();
  480.         final T temp  = join(foE.multiply(1.4), zero, zero.newInstance(1000.0), foE.subtract(2.0));
  481.         final T temp2 = join(zero, temp, zero.newInstance(1000.0), foE.subtract(temp));
  482.         final T value = join(temp2, temp2.multiply(0.85), zero.newInstance(60.0), foF2.multiply(0.85).subtract(temp2));
  483.         if (value.getReal() < 1.0E-6) {
  484.             return zero;
  485.         } else {
  486.             return value;
  487.         }
  488.     }

  490.     /**
  491.      * This method allows the computation of the F2, F1 and E layer amplitudes.
  492.      * <p>
  493.      * The resulting element is an array having the following form:
  494.      * <ul>
  495.      * <li>double[0] = A1 → F2 layer amplitude
  496.      * <li>double[1] = A2 → F1 layer amplitude
  497.      * <li>double[2] = A3 → E  layer amplitude
  498.      * </ul>
  499.      * </p>
  500.      * @param nmE E layer maximum density in 10^11 m-3
  501.      * @param nmF1 F1 layer maximum density in 10^11 m-3
  502.      * @param foF1 F1 layer critical frequency in MHz
  503.      * @return a three components array containing the layer amplitudes
  504.      */
  505.     private T[] computeLayerAmplitudes(final T nmE, final T nmF1, final T foF1) {
  506.         // Zero
  507.         final T zero = nmE.getField().getZero();

  509.         // Initialize array
  510.         final T[] amplitude = MathArrays.buildArray(nmE.getField(), 3);

  512.         // F2 layer amplitude (Eq. 90)
  513.         final T a1 = nmF2.multiply(4.0);
  514.         amplitude[0] = a1;

  516.         // F1 and E layer amplitudes (Eq. 91 to 98)
  517.         if (foF1.getReal() < 0.5) {
  518.             amplitude[1] = zero;
  519.             amplitude[2] = nmE.subtract(epst(a1, hmF2, b2Bot, hmE)).multiply(4.0);
  520.         } else {
  521.             T a2a = zero;
  522.             T a3a = nmE.multiply(4.0);
  523.             for (int i = 0; i < 5; i++) {
  524.                 a2a = nmF1.subtract(epst(a1, hmF2, b2Bot, hmF1)).subtract(epst(a3a, hmE, beTop, hmF1)).multiply(4.0);
  525.                 a2a = join(a2a, nmF1.multiply(0.8), nmE.getField().getOne(), a2a.subtract(nmF1.multiply(0.8)));
  526.                 a3a = nmE.subtract(epst(a2a, hmF1, b1Bot, hmE)).subtract(epst(a1, hmF2, b2Bot, hmE)).multiply(4.0);
  527.             }
  528.             amplitude[1] = a2a;
  529.             amplitude[2] = join(a3a, zero.newInstance(0.05), zero.newInstance(60.0), a3a.subtract(0.005));
  530.         }

  532.         return amplitude;
  533.     }

  535.     /**
  536.      * A clipped exponential function.
  537.      * <p>
  538.      * This function, describe in section F.2.12.2 of the reference document, is
  539.      * recommanded for the computation of exponential values.
  540.      * </p>
  541.      * @param power power for exponential function
  542.      * @return clipped exponential value
  543.      */
  544.     private T clipExp(final T power) {
  545.         final T zero = power.getField().getZero();
  546.         if (power.getReal() > 80.0) {
  547.             return zero.newInstance(5.5406E34);
  548.         } else if (power.getReal() < -80) {
  549.             return zero.newInstance(1.8049E-35);
  550.         } else {
  551.             return FastMath.exp(power);
  552.         }
  553.     }

  555.     /**
  556.      * Allows smooth joining of functions f1 and f2
  557.      * (i.e. continuous first derivatives) at origin.
  558.      * <p>
  559.      * This function, describe in section F.2.12.1 of the reference document, is
  560.      * recommanded for computational efficiency.
  561.      * </p>
  562.      * @param dF1 first function
  563.      * @param dF2 second function
  564.      * @param dA width of transition region
  565.      * @param dX x value
  566.      * @return the computed value
  567.      */
  568.     T join(final T dF1, final T dF2, final T dA, final T dX) {
  569.         final T ee = clipExp(dA.multiply(dX));
  570.         return dF1.multiply(ee).add(dF2).divide(ee.add(1.0));
  571.     }

  573.     /**
  574.      * The Epstein function.
  575.      * <p>
  576.      * This function, describe in section 2.5.1 of the reference document, is used
  577.      * as a basis analytical function in NeQuick for the construction of the ionospheric layers.
  578.      * </p>
  579.      * @param x x parameter
  580.      * @param y y parameter
  581.      * @param z z parameter
  582.      * @param w w parameter
  583.      * @return value of the epstein function
  584.      */
  585.     private T epst(final T x, final T y,
  586.                    final T z, final T w) {
  587.         final T ex  = clipExp(w.subtract(y).divide(z));
  588.         final T opex = ex.add(1.0);
  589.         return x.multiply(ex).divide(opex.multiply(opex));

  591.     }

  593. }
Created with JaCoCo 0.8.12.202403310830