FieldNeQuickParameters.java

  1. /* Copyright 2002-2024 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;

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

  28. /**
  29.  * This class perfoms the computation of the parameters used by the NeQuick model.
  30.  *
  31.  * @author Bryan Cazabonne
  32.  *
  33.  * @see "European Union (2016). European GNSS (Galileo) Open Service-Ionospheric Correction
  34.  *       Algorithm for Galileo Single Frequency Users. 1.2."
  35.  *
  36.  * @since 10.1
  37.  */
  38. class FieldNeQuickParameters <T extends CalculusFieldElement<T>> {

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

  43.     /** F2 layer maximum density. */
  44.     private final T nmF2;

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

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

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

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

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

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

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

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

  70.     /** topside thickness parameter [km]. */
  71.     private final T h0;

  73.     /** Layer amplitudes. */
  74.     private final T[] amplitudes;

  76.     /**
  77.      * Build a new instance.
  78.      * @param dateTime current date time components
  79.      * @param flattenF2 F2 coefficients used by the F2 layer (flatten array)
  80.      * @param flattenFm3 Fm3 coefficients used by the F2 layer (flatten array)
  81.      * @param latitude latitude of a point along the integration path, in radians
  82.      * @param longitude longitude of a point along the integration path, in radians
  83.      * @param alpha effective ionisation level coefficients
  84.      * @param modipGrip modip grid
  85.      */
  86.     FieldNeQuickParameters(final DateTimeComponents dateTime,
  87.                            final double[] flattenF2, final double[] flattenFm3,
  88.                            final T latitude, final T longitude,
  89.                            final double[] alpha, final double[][] modipGrip) {

  91.         // Zero
  92.         final T zero = latitude.getField().getZero();

  94.         // MODIP in degrees
  95.         final T modip = computeMODIP(latitude, longitude, modipGrip);
  96.         // Effective ionisation level Az
  97.         final T az = computeAz(modip, alpha);
  98.         // Effective sunspot number (Eq. 19)
  99.         final T azr = FastMath.sqrt(az.subtract(63.7).multiply(1123.6).add(167273.0)).subtract(408.99);
  100.         // Date and Time components
  101.         final DateComponents date = dateTime.getDate();
  102.         final TimeComponents time = dateTime.getTime();
  103.         // Hours
  104.         final double hours  = time.getSecondsInUTCDay() / 3600.0;
  105.         // Effective solar zenith angle in radians
  106.         final T xeff = computeEffectiveSolarAngle(date.getMonth(), hours, latitude, longitude);

  108.         // E layer maximum density height in km (Eq. 78)
  109.         this.hmE = zero.newInstance(120.0);
  110.         // E layer critical frequency in MHz
  111.         final T foE = computefoE(date.getMonth(), az, xeff, latitude);
  112.         // E layer maximum density in 10^11 m-3 (Eq. 36)
  113.         final T nmE = foE.multiply(foE).multiply(0.124);

  115.         // Time argument (Eq. 49)
  116.         final double t = FastMath.toRadians(15 * hours) - FastMath.PI;
  117.         // Compute Fourier time series for foF2 and M(3000)F2
  118.         final T[] scT = sinCos(zero.newInstance(t), 6);
  119.         final T[] cf2 = computeCF2(flattenF2, azr, scT);
  120.         final T[] cm3 = computeCm3(flattenFm3, azr, scT);
  121.         // F2 layer critical frequency in MHz
  122.         final T[] scL = sinCos(longitude, 8);
  123.         final T foF2 = computefoF2(modip, cf2, latitude, scL);
  124.         // Maximum Usable Frequency factor
  125.         final T mF2  = computeMF2(modip, cm3, latitude, scL);
  126.         // F2 layer maximum density in 10^11 m-3
  127.         this.nmF2 = foF2.multiply(foF2).multiply(0.124);
  128.         // F2 layer maximum density height in km
  129.         this.hmF2 = computehmF2(foE, foF2, mF2);

  131.         // F1 layer critical frequency in MHz
  132.         final T foF1 = computefoF1(foE, foF2);
  133.         // F1 layer maximum density in 10^11 m-3
  134.         final T nmF1;
  135.         if (foF1.getReal() <= 0.0 && foE.getReal() > 2.0) {
  136.             final T foEpopf = foE.add(0.5);
  137.             nmF1 = foEpopf.multiply(foEpopf).multiply(0.124);
  138.         } else {
  139.             nmF1 = foF1.multiply(foF1).multiply(0.124);
  140.         }
  141.         // F1 layer maximum density height in km
  142.         this.hmF1 = hmF2.add(hmE).multiply(0.5);

  144.         // Thickness parameters (Eq. 85 to 89)
  145.         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);
  146.         this.b2Bot = nmF2.divide(a).multiply(0.385);
  147.         this.b1Top = hmF2.subtract(hmF1).multiply(0.3);
  148.         this.b1Bot = hmF1.subtract(hmE).multiply(0.5);
  149.         this.beTop = FastMath.max(b1Bot, zero.newInstance(7.0));
  150.         this.beBot = zero.newInstance(5.0);

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

  155.         // Topside thickness parameter
  156.         this.h0 = computeH0(date.getMonth(), azr);
  157.     }

  159.     /**
  160.      * Get the F2 layer maximum density.
  161.      * @return nmF2
  162.      */
  163.     public T getNmF2() {
  164.         return nmF2;
  165.     }

  167.     /**
  168.      * Get the F2 layer maximum density height.
  169.      * @return hmF2 in km
  170.      */
  171.     public T getHmF2() {
  172.         return hmF2;
  173.     }

  175.     /**
  176.      * Get the F1 layer maximum density height.
  177.      * @return hmF1 in km
  178.      */
  179.     public T getHmF1() {
  180.         return hmF1;
  181.     }

  183.     /**
  184.      * Get the E layer maximum density height.
  185.      * @return hmE in km
  186.      */
  187.     public T getHmE() {
  188.         return hmE;
  189.     }

  191.     /**
  192.      * Get the F2 layer thickness parameter (bottom).
  193.      * @return B2Bot in km
  194.      */
  195.     public T getB2Bot() {
  196.         return b2Bot;
  197.     }

  199.     /**
  200.      * Get the F1 layer thickness parameter (top).
  201.      * @return B1Top in km
  202.      */
  203.     public T getB1Top() {
  204.         return b1Top;
  205.     }

  207.     /**
  208.      * Get the F1 layer thickness parameter (bottom).
  209.      * @return B1Bot in km
  210.      */
  211.     public T getB1Bot() {
  212.         return b1Bot;
  213.     }

  215.     /**
  216.      * Get the E layer thickness parameter (bottom).
  217.      * @return BeBot in km
  218.      */
  219.     public T getBEBot() {
  220.         return beBot;
  221.     }

  223.     /**
  224.      * Get the E layer thickness parameter (top).
  225.      * @return BeTop in km
  226.      */
  227.     public T getBETop() {
  228.         return beTop;
  229.     }

  231.     /**
  232.      * Get the F2, F1 and E layer amplitudes.
  233.      * <p>
  234.      * The resulting element is an array having the following form:
  235.      * <ul>
  236.      * <li>double[0] = A1 → F2 layer amplitude
  237.      * <li>double[1] = A2 → F1 layer amplitude
  238.      * <li>double[2] = A3 → E  layer amplitude
  239.      * </ul>
  240.      * @return layer amplitudes
  241.      */
  242.     public T[] getLayerAmplitudes() {
  243.         return amplitudes.clone();
  244.     }

  246.     /**
  247.      * Get the topside thickness parameter H0.
  248.      * @return H0 in km
  249.      */
  250.     public T getH0() {
  251.         return h0;
  252.     }

  254.     /**
  255.      * Computes the value of the modified dip latitude (MODIP) for the
  256.      * given latitude and longitude.
  257.      *
  258.      * @param lat receiver latitude, radians
  259.      * @param lon receiver longitude, radians
  260.      * @param stModip modip grid
  261.      * @return the MODIP in degrees
  262.      */
  263.     private T computeMODIP(final T lat, final T lon, final double[][] stModip) {

  265.         // Zero
  266.         final T zero = lat.getField().getZero();

  268.         // For the MODIP computation, the latitude and longitude have to be converted in degrees
  269.         final T latitude  = FastMath.toDegrees(lat);
  270.         final T longitude = FastMath.toDegrees(lon);

  272.         // Extreme cases
  273.         if (latitude.getReal() == 90.0 || latitude.getReal() == -90.0) {
  274.             return latitude;
  275.         }

  277.         // Auxiliary parameter l (Eq. 6 to 8)
  278.         final int lF = (int) ((longitude.getReal() + 180) * 0.1);
  279.         int l = lF - 2;
  280.         if (l < -2) {
  281.             l += 36;
  282.         } else if (l > 33) {
  283.             l -= 36;
  284.         }

  286.         // Auxiliary parameter a (Eq. 9 to 11)
  287.         final T a  = latitude.add(90).multiply(0.2).add(1.0);
  288.         final T aF = FastMath.floor(a);
  289.         // Eq. 10
  290.         final T x = a.subtract(aF);
  291.         // Eq. 11
  292.         final int i = (int) aF.getReal() - 2;

  294.         // zi coefficients (Eq. 12 and 13)
  295.         final T z1 = interpolate(zero.add(stModip[i + 1][l + 2]), zero.add(stModip[i + 2][l + 2]),
  296.                                       zero.add(stModip[i + 3][l + 2]), zero.add(stModip[i + 4][l + 2]), x);
  297.         final T z2 = interpolate(zero.add(stModip[i + 1][l + 3]), zero.add(stModip[i + 2][l + 3]),
  298.                                       zero.add(stModip[i + 3][l + 3]), zero.add(stModip[i + 4][l + 3]), x);
  299.         final T z3 = interpolate(zero.add(stModip[i + 1][l + 4]), zero.add(stModip[i + 2][l + 4]),
  300.                                       zero.add(stModip[i + 3][l + 4]), zero.add(stModip[i + 4][l + 4]), x);
  301.         final T z4 = interpolate(zero.add(stModip[i + 1][l + 5]), zero.add(stModip[i + 2][l + 5]),
  302.                                       zero.add(stModip[i + 3][l + 5]), zero.add(stModip[i + 4][l + 5]), x);

  304.         // Auxiliary parameter b (Eq. 14 and 15)
  305.         final T b  = longitude.add(180).multiply(0.1);
  306.         final T bF = FastMath.floor(b);
  307.         final T y  = b.subtract(bF);

  309.         // MODIP (Ref Eq. 16)
  310.         return interpolate(z1, z2, z3, z4, y);

  312.     }

  314.     /**
  315.      * This method computes the effective ionisation level Az.
  316.      * <p>
  317.      * This parameter is used for the computation of the Total Electron Content (TEC).
  318.      * </p>
  319.      * @param modip modified dip latitude (MODIP) in degrees
  320.      * @param alpha effective ionisation level coefficients
  321.      * @return the ionisation level Az
  322.      */
  323.     private T computeAz(final T modip, final double[] alpha) {
  324.         // Field
  325.         final Field<T> field = modip.getField();
  326.         // Zero
  327.         final T zero = field.getZero();
  328.         // Particular condition (Eq. 17)
  329.         if (alpha[0] == 0.0 && alpha[1] == 0.0 && alpha[2] == 0.0) {
  330.             return zero.newInstance(63.7);
  331.         }
  332.         // Az = a0 + modip * a1 + modip^2 * a2 (Eq. 18)
  333.         T az = modip.multiply(alpha[2]).add(alpha[1]).multiply(modip).add(alpha[0]);
  334.         // If Az < 0 -> Az = 0
  335.         az = FastMath.max(zero, az);
  336.         // If Az > 400 -> Az = 400
  337.         az = FastMath.min(zero.newInstance(400.0), az);
  338.         return az;
  339.     }

  341.     /**
  342.      * This method computes the effective solar zenith angle.
  343.      * <p>
  344.      * The effective solar zenith angle is compute as a function of the
  345.      * solar zenith angle and the solar zenith angle at day night transition.
  346.      * </p>
  347.      * @param month current month of the year
  348.      * @param hours universal time (hours)
  349.      * @param latitude in radians
  350.      * @param longitude in radians
  351.      * @return the effective solar zenith angle, radians
  352.      */
  353.     private T computeEffectiveSolarAngle(final int month,
  354.                                          final double hours,
  355.                                          final T latitude,
  356.                                          final T longitude) {
  357.         // Zero
  358.         final T zero = latitude.getField().getZero();
  359.         // Local time (Eq.4)
  360.         final T lt = longitude.divide(FastMath.toRadians(15.0)).add(hours);
  361.         // Day of year at the middle of the month (Eq. 20)
  362.         final double dy = 30.5 * month - 15.0;
  363.         // Time (Eq. 21)
  364.         final double t = dy + (18 - hours) / 24;
  365.         // Arguments am and al (Eq. 22 and 23)
  366.         final double am = FastMath.toRadians(0.9856 * t - 3.289);
  367.         final double al = am + FastMath.toRadians(1.916 * FastMath.sin(am) + 0.020 * FastMath.sin(2.0 * am) + 282.634);
  368.         // Sine and cosine of solar declination (Eq. 24 and 25)
  369.         final double sDec = 0.39782 * FastMath.sin(al);
  370.         final double cDec = FastMath.sqrt(1. - sDec * sDec);
  371.         // Solar zenith angle, deg (Eq. 26 and 27)
  372.         final FieldSinCos<T> scLat   = FastMath.sinCos(latitude);
  373.         final T coef    = lt.negate().add(12.0).multiply(FastMath.PI / 12);
  374.         final T cZenith = scLat.sin().multiply(sDec).add(scLat.cos().multiply(cDec).multiply(FastMath.cos(coef)));
  375.         final T angle   = FastMath.atan2(FastMath.sqrt(cZenith.multiply(cZenith).negate().add(1.0)), cZenith);
  376.         final T x       = FastMath.toDegrees(angle);
  377.         // Effective solar zenith angle (Eq. 28)
  378.         final T xeff = join(clipExp(x.multiply(0.2).negate().add(20.0)).multiply(0.24).negate().add(90.0), x,
  379.                 zero.newInstance(12.0), x.subtract(X0));
  380.         return FastMath.toRadians(xeff);
  381.     }

  383.     /**
  384.      * This method computes the E layer critical frequency at a given location.
  385.      * @param month current month
  386.      * @param az ffective ionisation level
  387.      * @param xeff effective solar zenith angle in radians
  388.      * @param latitude latitude in radians
  389.      * @return the E layer critical frequency at a given location in MHz
  390.      */
  391.     private T computefoE(final int month, final T az,
  392.                          final T xeff, final T latitude) {
  393.         // The latitude has to be converted in degrees
  394.         final T lat = FastMath.toDegrees(latitude);
  395.         // Square root of the effective ionisation level
  396.         final T sqAz = FastMath.sqrt(az);
  397.         // seas parameter (Eq. 30 to 32)
  398.         final int seas;
  399.         if (month == 1 || month == 2 || month == 11 || month == 12) {
  400.             seas = -1;
  401.         } else if (month == 3 || month == 4 || month == 9 || month == 10) {
  402.             seas = 0;
  403.         } else {
  404.             seas = 1;
  405.         }
  406.         // Latitudinal dependence (Eq. 33 and 34)
  407.         final T ee = clipExp(lat.multiply(0.3));
  408.         final T seasp = ee.subtract(1.0).divide(ee.add(1.0)).multiply(seas);
  409.         // Critical frequency (Eq. 35)
  410.         final T coef = seasp.multiply(0.019).negate().add(1.112);
  411.         return FastMath.sqrt(coef .multiply(coef).multiply(sqAz).multiply(FastMath.cos(xeff).pow(0.6)).add(0.49));

  413.     }

  415.     /**
  416.      * Computes the F2 layer height of maximum electron density.
  417.      * @param foE E layer layer critical frequency in MHz
  418.      * @param foF2 F2 layer layer critical frequency in MHz
  419.      * @param mF2 maximum usable frequency factor
  420.      * @return hmF2 in km
  421.      */
  422.     private T computehmF2(final T foE, final T foF2, final T mF2) {
  423.         // Zero
  424.         final T zero = foE.getField().getZero();
  425.         // Ratio
  426.         final T fo = foF2.divide(foE);
  427.         final T ratio = join(fo, zero.newInstance(1.75), zero.newInstance(20.0), fo.subtract(1.75));

  429.         // deltaM parameter
  430.         T deltaM = zero.subtract(0.012);
  431.         if (foE.getReal() >= 1e-30) {
  432.             deltaM = deltaM.add(ratio.subtract(1.215).divide(0.253).reciprocal());
  433.         }

  435.         // hmF2 Eq. 80
  436.         final T mF2Sq = mF2.square();
  437.         final T temp  = FastMath.sqrt(mF2Sq.multiply(0.0196).add(1.0).divide(mF2Sq.multiply(1.2967).subtract(1.0)));
  438.         return mF2.multiply(1490.0).multiply(temp).divide(mF2.add(deltaM)).subtract(176.0);

  440.     }

  442.     /** Compute sines and cosines.
  443.      * @param a argument
  444.      * @param n number of terms
  445.      * @return sin(a), cos(a), sin(2a), cos(2a) … sin(n a), cos(n a) array
  446.      * @since 12.1.3
  447.      */
  448.     private T[] sinCos(final T a, final int n) {

  450.         final FieldSinCos<T> sc0 = FastMath.sinCos(a);
  451.         FieldSinCos<T> sci = sc0;
  452.         final T[] sc = MathArrays.buildArray(a.getField(), 2 * n);
  453.         int isc = 0;
  454.         sc[isc++] = sci.sin();
  455.         sc[isc++] = sci.cos();
  456.         for (int i = 1; i < n; i++) {
  457.             sci = FieldSinCos.sum(sc0, sci);
  458.             sc[isc++] = sci.sin();
  459.             sc[isc++] = sci.cos();
  460.         }

  462.         return sc;

  464.     }

  466.     /**
  467.      * Computes cf2 coefficients.
  468.      * @param flattenF2 F2 coefficients used by the F2 layer (flatten array)
  469.      * @param azr effective sunspot number (Eq. 19)
  470.      * @param scT sines/cosines array of time argument
  471.      * @return the cf2 coefficients array
  472.      */
  473.     private T[] computeCF2(final double[] flattenF2, final T azr, final T[] scT) {

  475.         // interpolation coefficients for effective spot number
  476.         final T azr01 = azr.multiply(0.01);
  477.         final T omazr01 = azr01.negate().add(1);

  479.         // Eq. 44 and Eq. 50 merged into one loop
  480.         final T[] cf2 = MathArrays.buildArray(azr.getField(), 76);
  481.         int index = 0;
  482.         for (int i = 0; i < cf2.length; i++) {
  483.             // CHECKSTYLE: stop Indentation check
  484.             cf2[i] = omazr01.multiply(flattenF2[index     ]).add(azr01.multiply(flattenF2[index +  1])).
  485.                  add(omazr01.multiply(flattenF2[index +  2]).add(azr01.multiply(flattenF2[index +  3])).multiply(scT[ 0])).
  486.                  add(omazr01.multiply(flattenF2[index +  4]).add(azr01.multiply(flattenF2[index +  5])).multiply(scT[ 1])).
  487.                  add(omazr01.multiply(flattenF2[index +  6]).add(azr01.multiply(flattenF2[index +  7])).multiply(scT[ 2])).
  488.                  add(omazr01.multiply(flattenF2[index +  8]).add(azr01.multiply(flattenF2[index +  9])).multiply(scT[ 3])).
  489.                  add(omazr01.multiply(flattenF2[index + 10]).add(azr01.multiply(flattenF2[index + 11])).multiply(scT[ 4])).
  490.                  add(omazr01.multiply(flattenF2[index + 12]).add(azr01.multiply(flattenF2[index + 13])).multiply(scT[ 5])).
  491.                  add(omazr01.multiply(flattenF2[index + 14]).add(azr01.multiply(flattenF2[index + 15])).multiply(scT[ 6])).
  492.                  add(omazr01.multiply(flattenF2[index + 16]).add(azr01.multiply(flattenF2[index + 17])).multiply(scT[ 7])).
  493.                  add(omazr01.multiply(flattenF2[index + 18]).add(azr01.multiply(flattenF2[index + 19])).multiply(scT[ 8])).
  494.                  add(omazr01.multiply(flattenF2[index + 20]).add(azr01.multiply(flattenF2[index + 21])).multiply(scT[ 9])).
  495.                  add(omazr01.multiply(flattenF2[index + 22]).add(azr01.multiply(flattenF2[index + 23])).multiply(scT[10])).
  496.                  add(omazr01.multiply(flattenF2[index + 24]).add(azr01.multiply(flattenF2[index + 25])).multiply(scT[11]));
  497.             index += 26;
  498.             // CHECKSTYLE: resume Indentation check
  499.         }
  500.         return cf2;
  501.     }

  503.     /**
  504.      * Computes Cm3 coefficients.
  505.      * @param flattenFm3 Fm3 coefficients used by the F2 layer (flatten array)
  506.      * @param azr effective sunspot number (Eq. 19)
  507.      * @param scT sines/cosines array of time argument
  508.      * @return the Cm3 coefficients array
  509.      */
  510.     private T[] computeCm3(final double[] flattenFm3, final T azr, final T[] scT) {

  512.         // interpolation coefficients for effective spot number
  513.         final T azr01 = azr.multiply(0.01);
  514.         final T omazr01 = azr01.negate().add(1);

  516.         // Eq. 44 and Eq. 51 merged into one loop
  517.         final T[] cm3 = MathArrays.buildArray(azr.getField(), 49);
  518.         int index = 0;
  519.         for (int i = 0; i < cm3.length; i++) {
  520.             cm3[i] = omazr01.multiply(flattenFm3[index     ]).add(azr01.multiply(flattenFm3[index +  1])).
  521.                  add(omazr01.multiply(flattenFm3[index +  2]).add(azr01.multiply(flattenFm3[index +  3])).multiply(scT[ 0])).
  522.                  add(omazr01.multiply(flattenFm3[index +  4]).add(azr01.multiply(flattenFm3[index +  5])).multiply(scT[ 1])).
  523.                  add(omazr01.multiply(flattenFm3[index +  6]).add(azr01.multiply(flattenFm3[index +  7])).multiply(scT[ 2])).
  524.                  add(omazr01.multiply(flattenFm3[index +  8]).add(azr01.multiply(flattenFm3[index +  9])).multiply(scT[ 3])).
  525.                  add(omazr01.multiply(flattenFm3[index + 10]).add(azr01.multiply(flattenFm3[index + 11])).multiply(scT[ 4])).
  526.                  add(omazr01.multiply(flattenFm3[index + 12]).add(azr01.multiply(flattenFm3[index + 13])).multiply(scT[ 5])).
  527.                  add(omazr01.multiply(flattenFm3[index + 14]).add(azr01.multiply(flattenFm3[index + 15])).multiply(scT[ 6])).
  528.                  add(omazr01.multiply(flattenFm3[index + 16]).add(azr01.multiply(flattenFm3[index + 17])).multiply(scT[ 7]));
  529.             index += 18;
  530.         }
  531.         return cm3;
  532.     }

  534.     /**
  535.      * This method computes the F2 layer critical frequency.
  536.      * @param modip modified DIP latitude, in degrees
  537.      * @param cf2 Fourier time series for foF2
  538.      * @param latitude latitude in radians
  539.      * @param scL sines/cosines array of longitude argument
  540.      * @return the F2 layer critical frequency, MHz
  541.      */
  542.     private T computefoF2(final T modip, final T[] cf2,
  543.                           final T latitude, final T[] scL) {

  545.         // Legendre grades (Eq. 63)
  546.         final int[] q = new int[] {
  547.             12, 12, 9, 5, 2, 1, 1, 1, 1
  548.         };

  550.         T frequency = cf2[0];

  552.         // MODIP coefficients Eq. 57
  553.         final T sinMODIP = FastMath.sin(FastMath.toRadians(modip));
  554.         final T[] m = MathArrays.buildArray(latitude.getField(), 12);
  555.         m[0] = latitude.getField().getOne();
  556.         for (int i = 1; i < q[0]; i++) {
  557.             m[i] = sinMODIP.multiply(m[i - 1]);
  558.             frequency = frequency.add(m[i].multiply(cf2[i]));
  559.         }

  561.         // latitude and longitude terms
  562.         int index = 12;
  563.         final T cosLat1 = FastMath.cos(latitude);
  564.         T cosLatI = cosLat1;
  565.         for (int i = 1; i < q.length; i++) {
  566.             final T c = cosLatI.multiply(scL[2 * i - 1]);
  567.             final T s = cosLatI.multiply(scL[2 * i - 2]);
  568.             for (int j = 0; j < q[i]; j++) {
  569.                 frequency = frequency.add(m[j].multiply(c).multiply(cf2[index++]));
  570.                 frequency = frequency.add(m[j].multiply(s).multiply(cf2[index++]));
  571.             }
  572.             cosLatI = cosLatI.multiply(cosLat1);
  573.         }

  575.         return frequency;

  577.     }

  579.     /**
  580.      * This method computes the Maximum Usable Frequency factor.
  581.      * @param modip modified DIP latitude, in degrees
  582.      * @param cm3 Fourier time series for M(3000)F2
  583.      * @param latitude latitude in radians
  584.      * @param scL sines/cosines array of longitude argument
  585.      * @return the Maximum Usable Frequency factor
  586.      */
  587.     private T computeMF2(final T modip, final T[] cm3,
  588.                          final T latitude, final T[] scL) {

  590.         // Legendre grades (Eq. 71)
  591.         final int[] r = new int[] {
  592.             7, 8, 6, 3, 2, 1, 1
  593.         };

  595.         T m3000 = cm3[0];

  597.         // MODIP coefficients Eq. 57
  598.         final T sinMODIP = FastMath.sin(FastMath.toRadians(modip));
  599.         final T[] m = MathArrays.buildArray(latitude.getField(), 12);
  600.         m[0] = latitude.getField().getOne();
  601.         for (int i = 1; i < 12; i++) {
  602.             m[i] = sinMODIP.multiply(m[i - 1]);
  603.             if (i < 7) {
  604.                 m3000 = m3000.add(m[i].multiply(cm3[i]));
  605.             }
  606.         }

  608.         // latitude and longitude terms
  609.         int index = 7;
  610.         final T cosLat1 = FastMath.cos(latitude);
  611.         T cosLatI = cosLat1;
  612.         for (int i = 1; i < r.length; i++) {
  613.             final T c = cosLatI.multiply(scL[2 * i - 1]);
  614.             final T s = cosLatI.multiply(scL[2 * i - 2]);
  615.             for (int j = 0; j < r[i]; j++) {
  616.                 m3000 = m3000.add(m[j].multiply(c).multiply(cm3[index++]));
  617.                 m3000 = m3000.add(m[j].multiply(s).multiply(cm3[index++]));
  618.             }
  619.             cosLatI = cosLatI.multiply(cosLat1); // Eq. 58
  620.         }

  622.         return m3000;

  624.     }

  626.     /**
  627.      * This method computes the F1 layer critical frequency.
  628.      * <p>
  629.      * This computation performs the algorithm exposed in Annex F
  630.      * of the reference document.
  631.      * </p>
  632.      * @param foE the E layer critical frequency, MHz
  633.      * @return the F1 layer critical frequency, MHz
  634.      * @param foF2 the F2 layer critical frequency, MHz
  635.      */
  636.     private T computefoF1(final T foE, final T foF2) {
  637.         final T zero = foE.getField().getZero();
  638.         final T temp  = join(foE.multiply(1.4), zero, zero.newInstance(1000.0), foE.subtract(2.0));
  639.         final T temp2 = join(zero, temp, zero.newInstance(1000.0), foE.subtract(temp));
  640.         final T value = join(temp2, temp2.multiply(0.85), zero.newInstance(60.0), foF2.multiply(0.85).subtract(temp2));
  641.         if (value.getReal() < 1.0E-6) {
  642.             return zero;
  643.         } else {
  644.             return value;
  645.         }
  646.     }

  648.     /**
  649.      * This method allows the computation of the F2, F1 and E layer amplitudes.
  650.      * <p>
  651.      * The resulting element is an array having the following form:
  652.      * <ul>
  653.      * <li>double[0] = A1 → F2 layer amplitude
  654.      * <li>double[1] = A2 → F1 layer amplitude
  655.      * <li>double[2] = A3 → E  layer amplitude
  656.      * </ul>
  657.      * </p>
  658.      * @param nmE E layer maximum density in 10^11 m-3
  659.      * @param nmF1 F1 layer maximum density in 10^11 m-3
  660.      * @param foF1 F1 layer critical frequency in MHz
  661.      * @return a three components array containing the layer amplitudes
  662.      */
  663.     private T[] computeLayerAmplitudes(final T nmE, final T nmF1, final T foF1) {
  664.         // Zero
  665.         final T zero = nmE.getField().getZero();

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

  670.         // F2 layer amplitude (Eq. 90)
  671.         final T a1 = nmF2.multiply(4.0);
  672.         amplitude[0] = a1;

  674.         // F1 and E layer amplitudes (Eq. 91 to 98)
  675.         if (foF1.getReal() < 0.5) {
  676.             amplitude[1] = zero;
  677.             amplitude[2] = nmE.subtract(epst(a1, hmF2, b2Bot, hmE)).multiply(4.0);
  678.         } else {
  679.             T a2a = zero;
  680.             T a3a = nmE.multiply(4.0);
  681.             for (int i = 0; i < 5; i++) {
  682.                 a2a = nmF1.subtract(epst(a1, hmF2, b2Bot, hmF1)).subtract(epst(a3a, hmE, beTop, hmF1)).multiply(4.0);
  683.                 a2a = join(a2a, nmF1.multiply(0.8), nmE.getField().getOne(), a2a.subtract(nmF1.multiply(0.8)));
  684.                 a3a = nmE.subtract(epst(a2a, hmF1, b1Bot, hmE)).subtract(epst(a1, hmF2, b2Bot, hmE)).multiply(4.0);
  685.             }
  686.             amplitude[1] = a2a;
  687.             amplitude[2] = join(a3a, zero.newInstance(0.05), zero.newInstance(60.0), a3a.subtract(0.005));
  688.         }

  690.         return amplitude;
  691.     }

  693.     /**
  694.      * This method computes the topside thickness parameter H0.
  695.      *
  696.      * @param month current month
  697.      * @param azr effective sunspot number
  698.      * @return H0 in km
  699.      */
  700.     private T computeH0(final int month, final T azr) {

  702.         // One
  703.         final T one = azr.getField().getOne();

  705.         // Auxiliary parameter ka (Eq. 99 and 100)
  706.         final T ka;
  707.         if (month > 3 && month < 10) {
  708.             // month = 4,5,6,7,8,9
  709.             ka = azr.multiply(0.014).add(hmF2.multiply(0.008)).negate().add(6.705);
  710.         } else {
  711.             // month = 1,2,3,10,11,12
  712.             final T ratio = hmF2.divide(b2Bot);
  713.             ka = ratio.multiply(ratio).multiply(0.097).add(nmF2.multiply(0.153)).add(-7.77);
  714.         }

  716.         // Auxiliary parameter kb (Eq. 101 and 102)
  717.         T kb = join(ka, one.newInstance(2.0), one, ka.subtract(2.0));
  718.         kb = join(one.newInstance(8.0), kb, one, kb.subtract(8.0));

  720.         // Auxiliary parameter Ha (Eq. 103)
  721.         final T hA = kb.multiply(b2Bot);

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

  727.         // Topside thickness parameter (Eq. 106)
  728.         return hA.divide(v);

  730.     }

  732.     /**
  733.      * A clipped exponential function.
  734.      * <p>
  735.      * This function, describe in section F.2.12.2 of the reference document, is
  736.      * recommanded for the computation of exponential values.
  737.      * </p>
  738.      * @param power power for exponential function
  739.      * @return clipped exponential value
  740.      */
  741.     private T clipExp(final T power) {
  742.         final T zero = power.getField().getZero();
  743.         if (power.getReal() > 80.0) {
  744.             return zero.newInstance(5.5406E34);
  745.         } else if (power.getReal() < -80) {
  746.             return zero.newInstance(1.8049E-35);
  747.         } else {
  748.             return FastMath.exp(power);
  749.         }
  750.     }

  752.     /**
  753.      * This method provides a third order interpolation function
  754.      * as recommended in the reference document (Ref Eq. 128 to Eq. 138)
  755.      *
  756.      * @param z1 z1 coefficient
  757.      * @param z2 z2 coefficient
  758.      * @param z3 z3 coefficient
  759.      * @param z4 z4 coefficient
  760.      * @param x position
  761.      * @return a third order interpolation
  762.      */
  763.     private T interpolate(final T z1, final T z2,
  764.                           final T z3, final T z4,
  765.                           final T x) {

  767.         if (FastMath.abs(2.0 * x.getReal()) < 1e-10) {
  768.             return z2;
  769.         }

  771.         final T delta = x.multiply(2.0).subtract(1.0);
  772.         final T g1 = z3.add(z2);
  773.         final T g2 = z3.subtract(z2);
  774.         final T g3 = z4.add(z1);
  775.         final T g4 = z4.subtract(z1).divide(3.0);
  776.         final T a0 = g1.multiply(9.0).subtract(g3);
  777.         final T a1 = g2.multiply(9.0).subtract(g4);
  778.         final T a2 = g3.subtract(g1);
  779.         final T a3 = g4.subtract(g2);
  780.         return delta.multiply(a3).add(a2).multiply(delta).add(a1).multiply(delta).add(a0).multiply(0.0625);

  782.     }

  784.     /**
  785.      * Allows smooth joining of functions f1 and f2
  786.      * (i.e. continuous first derivatives) at origin.
  787.      * <p>
  788.      * This function, describe in section F.2.12.1 of the reference document, is
  789.      * recommanded for computational efficiency.
  790.      * </p>
  791.      * @param dF1 first function
  792.      * @param dF2 second function
  793.      * @param dA width of transition region
  794.      * @param dX x value
  795.      * @return the computed value
  796.      */
  797.     private T join(final T dF1, final T dF2,
  798.                    final T dA, final T dX) {
  799.         final T ee = clipExp(dA.multiply(dX));
  800.         return dF1.multiply(ee).add(dF2).divide(ee.add(1.0));
  801.     }

  803.     /**
  804.      * The Epstein function.
  805.      * <p>
  806.      * This function, describe in section 2.5.1 of the reference document, is used
  807.      * as a basis analytical function in NeQuick for the construction of the ionospheric layers.
  808.      * </p>
  809.      * @param x x parameter
  810.      * @param y y parameter
  811.      * @param z z parameter
  812.      * @param w w parameter
  813.      * @return value of the epstein function
  814.      */
  815.     private T epst(final T x, final T y,
  816.                    final T z, final T w) {
  817.         final T ex  = clipExp(w.subtract(y).divide(z));
  818.         final T opex = ex.add(1.0);
  819.         return x.multiply(ex).divide(opex.multiply(opex));

  821.     }

  823. }
Created with JaCoCo 0.8.12.202403310830