NeQuickParameters.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.util.FastMath;
  20. import org.hipparchus.util.SinCos;
  21. import org.orekit.time.DateComponents;
  22. import org.orekit.time.DateTimeComponents;
  23. import org.orekit.time.TimeComponents;

  25. /**
  26.  * This class perfoms the computation of the parameters used by the NeQuick model.
  27.  *
  28.  * @author Bryan Cazabonne
  29.  *
  30.  * @see "European Union (2016). European GNSS (Galileo) Open Service-Ionospheric Correction
  31.  *       Algorithm for Galileo Single Frequency Users. 1.2."
  32.  *
  33.  * @since 10.1
  34.  */
  35. class NeQuickParameters {

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

  40.     /** F2 layer maximum density. */
  41.     private final double nmF2;

  43.     /** F2 layer maximum density height [km]. */
  44.     private final double hmF2;

  46.     /** F1 layer maximum density height [km]. */
  47.     private final double hmF1;

  49.     /** E layer maximum density height [km]. */
  50.     private final double hmE;

  52.     /** F2 layer bottom thickness parameter [km]. */
  53.     private final double b2Bot;

  55.     /** F1 layer top thickness parameter [km]. */
  56.     private final double b1Top;

  58.     /** F1 layer bottom thickness parameter [km]. */
  59.     private final double b1Bot;

  61.     /** E layer top thickness parameter [km]. */
  62.     private final double beTop;

  64.     /** E layer bottom thickness parameter [km]. */
  65.     private final double beBot;

  67.     /** topside thickness parameter [km]. */
  68.     private final double h0;

  70.     /** Layer amplitudes. */
  71.     private final double[] amplitudes;

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

  88.         // MODIP in degrees
  89.         final double modip = computeMODIP(latitude, longitude, modipGrip);
  90.         // Effective ionisation level Az
  91.         final double az = computeAz(modip, alpha);
  92.         // Effective sunspot number (Eq. 19)
  93.         final double azr = FastMath.sqrt(167273.0 + (az - 63.7) * 1123.6) - 408.99;
  94.         // Date and Time components
  95.         final DateComponents date = dateTime.getDate();
  96.         final TimeComponents time = dateTime.getTime();
  97.         // Hours
  98.         final double hours  = time.getSecondsInUTCDay() / 3600.0;
  99.         // Effective solar zenith angle in radians
  100.         final double xeff = computeEffectiveSolarAngle(date.getMonth(), hours, latitude, longitude);

  102.         // E layer maximum density height in km (Eq. 78)
  103.         this.hmE = 120.0;
  104.         // E layer critical frequency in MHz
  105.         final double foE = computefoE(date.getMonth(), az, xeff, latitude);
  106.         // E layer maximum density in 10^11 m-3 (Eq. 36)
  107.         final double nmE = 0.124 * foE * foE;

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

  125.         // F1 layer critical frequency in MHz
  126.         final double foF1 = computefoF1(foE, foF2);
  127.         // F1 layer maximum density in 10^11 m-3
  128.         final double nmF1;
  129.         if (foF1 <= 0.0 && foE > 2.0) {
  130.             final double foEpopf = foE + 0.5;
  131.             nmF1 = 0.124 * foEpopf * foEpopf;
  132.         } else {
  133.             nmF1 = 0.124 * foF1 * foF1;
  134.         }
  135.         // F1 layer maximum density height in km
  136.         this.hmF1 = 0.5 * (hmF2 + hmE);

  138.         // Thickness parameters (Eq. 85 to 89)
  139.         final double a = 0.01 * clipExp(-3.467 + 0.857 * FastMath.log(foF2 * foF2) + 2.02 * FastMath.log(mF2));
  140.         this.b2Bot = 0.385 * nmF2 / a;
  141.         this.b1Top = 0.3 * (hmF2 - hmF1);
  142.         this.b1Bot = 0.5 * (hmF1 - hmE);
  143.         this.beTop = FastMath.max(b1Bot, 7.0);
  144.         this.beBot = 5.0;

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

  149.         // Topside thickness parameter
  150.         this.h0 = computeH0(date.getMonth(), azr);
  151.     }

  153.     /**
  154.      * Get the F2 layer maximum density.
  155.      * @return nmF2
  156.      */
  157.     public double getNmF2() {
  158.         return nmF2;
  159.     }

  161.     /**
  162.      * Get the F2 layer maximum density height.
  163.      * @return hmF2 in km
  164.      */
  165.     public double getHmF2() {
  166.         return hmF2;
  167.     }

  169.     /**
  170.      * Get the F1 layer maximum density height.
  171.      * @return hmF1 in km
  172.      */
  173.     public double getHmF1() {
  174.         return hmF1;
  175.     }

  177.     /**
  178.      * Get the E layer maximum density height.
  179.      * @return hmE in km
  180.      */
  181.     public double getHmE() {
  182.         return hmE;
  183.     }

  185.     /**
  186.      * Get the F2 layer thickness parameter (bottom).
  187.      * @return B2Bot in km
  188.      */
  189.     public double getB2Bot() {
  190.         return b2Bot;
  191.     }

  193.     /**
  194.      * Get the F1 layer thickness parameter (top).
  195.      * @return B1Top in km
  196.      */
  197.     public double getB1Top() {
  198.         return b1Top;
  199.     }

  201.     /**
  202.      * Get the F1 layer thickness parameter (bottom).
  203.      * @return B1Bot in km
  204.      */
  205.     public double getB1Bot() {
  206.         return b1Bot;
  207.     }

  209.     /**
  210.      * Get the E layer thickness parameter (bottom).
  211.      * @return BeBot in km
  212.      */
  213.     public double getBEBot() {
  214.         return beBot;
  215.     }

  217.     /**
  218.      * Get the E layer thickness parameter (top).
  219.      * @return BeTop in km
  220.      */
  221.     public double getBETop() {
  222.         return beTop;
  223.     }

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

  240.     /**
  241.      * Get the topside thickness parameter H0.
  242.      * @return H0 in km
  243.      */
  244.     public double getH0() {
  245.         return h0;
  246.     }

  248.     /**
  249.      * Computes the value of the modified dip latitude (MODIP) for the
  250.      * given latitude and longitude.
  251.      *
  252.      * @param lat receiver latitude, radians
  253.      * @param lon receiver longitude, radians
  254.      * @param stModip modip grid
  255.      * @return the MODIP in degrees
  256.      */
  257.     private double computeMODIP(final double lat, final double lon, final double[][] stModip) {

  259.         // For the MODIP computation, latitude and longitude have to be converted in degrees
  260.         final double latitude = FastMath.toDegrees(lat);
  261.         final double longitude = FastMath.toDegrees(lon);

  263.         // Extreme cases
  264.         if (latitude == 90.0 || latitude == -90.0) {
  265.             return latitude;
  266.         }

  268.         // Auxiliary parameter l (Eq. 6 to 8)
  269.         final int lF = (int) ((longitude + 180) * 0.1);
  270.         int l = lF - 2;
  271.         if (l < -2) {
  272.             l += 36;
  273.         } else if (l > 33) {
  274.             l -= 36;
  275.         }

  277.         // Auxiliary parameter a (Eq. 9 to 11)
  278.         final double a  = 0.2 * (latitude + 90) + 1.0;
  279.         final double aF = FastMath.floor(a);
  280.         // Eq. 10
  281.         final double x = a - aF;
  282.         // Eq. 11
  283.         final int i = (int) aF - 2;

  285.         // zi coefficients (Eq. 12 and 13)
  286.         final double z1 = interpolate(stModip[i + 1][l + 2], stModip[i + 2][l + 2], stModip[i + 3][l + 2], stModip[i + 4][l + 2], x);
  287.         final double z2 = interpolate(stModip[i + 1][l + 3], stModip[i + 2][l + 3], stModip[i + 3][l + 3], stModip[i + 4][l + 3], x);
  288.         final double z3 = interpolate(stModip[i + 1][l + 4], stModip[i + 2][l + 4], stModip[i + 3][l + 4], stModip[i + 4][l + 4], x);
  289.         final double z4 = interpolate(stModip[i + 1][l + 5], stModip[i + 2][l + 5], stModip[i + 3][l + 5], stModip[i + 4][l + 5], x);

  291.         // Auxiliary parameter b (Eq. 14 and 15)
  292.         final double b  = (longitude + 180) * 0.1;
  293.         final double bF = FastMath.floor(b);
  294.         final double y  = b - bF;

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

  299.     }

  301.     /**
  302.      * This method computes the effective ionisation level Az.
  303.      * <p>
  304.      * This parameter is used for the computation of the Total Electron Content (TEC).
  305.      * </p>
  306.      * @param modip modified dip latitude (MODIP) in degrees
  307.      * @param alpha effective ionisation level coefficients
  308.      * @return the ionisation level Az
  309.      */
  310.     private double computeAz(final double modip, final double[] alpha) {
  311.         // Particular condition (Eq. 17)
  312.         if (alpha[0] == 0.0 && alpha[1] == 0.0 && alpha[2] == 0.0) {
  313.             return 63.7;
  314.         }
  315.         // Az = a0 + modip * a1 + modip^2 * a2 (Eq. 18)
  316.         double az = alpha[0] + modip * (alpha[1] + modip * alpha[2]);
  317.         // If Az < 0 -> Az = 0
  318.         az = FastMath.max(0.0, az);
  319.         // If Az > 400 -> Az = 400
  320.         az = FastMath.min(400.0, az);
  321.         return az;
  322.     }

  324.     /**
  325.      * This method computes the effective solar zenith angle.
  326.      * <p>
  327.      * The effective solar zenith angle is compute as a function of the
  328.      * solar zenith angle and the solar zenith angle at day night transition.
  329.      * </p>
  330.      * @param month current month of the year
  331.      * @param hours universal time (hours)
  332.      * @param latitude in radians
  333.      * @param longitude in radians
  334.      * @return the effective solar zenith angle, radians
  335.      */
  336.     private double computeEffectiveSolarAngle(final int month,
  337.                                               final double hours,
  338.                                               final double latitude,
  339.                                               final double longitude) {
  340.         // Local time (Eq.4)
  341.         final double lt = hours + longitude / FastMath.toRadians(15.0);
  342.         // Day of year at the middle of the month (Eq. 20)
  343.         final double dy = 30.5 * month - 15.0;
  344.         // Time (Eq. 21)
  345.         final double t = dy + (18 - hours) / 24;
  346.         // Arguments am and al (Eq. 22 and 23)
  347.         final double am = FastMath.toRadians(0.9856 * t - 3.289);
  348.         final double al = am + FastMath.toRadians(1.916 * FastMath.sin(am) + 0.020 * FastMath.sin(2.0 * am) + 282.634);
  349.         // Sine and cosine of solar declination (Eq. 24 and 25)
  350.         final double sDec = 0.39782 * FastMath.sin(al);
  351.         final double cDec = FastMath.sqrt(1. - sDec * sDec);
  352.         // Solar zenith angle, deg (Eq. 26 and 27)
  353.         final SinCos scLat   = FastMath.sinCos(latitude);
  354.         final double coef    = (FastMath.PI / 12) * (12 - lt);
  355.         final double cZenith = scLat.sin() * sDec + scLat.cos() * cDec * FastMath.cos(coef);
  356.         final double angle   = FastMath.atan2(FastMath.sqrt(1.0 - cZenith * cZenith), cZenith);
  357.         final double x       = FastMath.toDegrees(angle);
  358.         // Effective solar zenith angle (Eq. 28)
  359.         final double xeff = join(90.0 - 0.24 * clipExp(20.0 - 0.2 * x), x, 12.0, x - X0);
  360.         return FastMath.toRadians(xeff);
  361.     }

  363.     /**
  364.      * This method computes the E layer critical frequency at a given location.
  365.      * @param month current month
  366.      * @param az ffective ionisation level
  367.      * @param xeff effective solar zenith angle in radians
  368.      * @param latitude latitude in radians
  369.      * @return the E layer critical frequency at a given location in MHz
  370.      */
  371.     private double computefoE(final int month, final double az,
  372.                               final double xeff, final double latitude) {
  373.         // The latitude has to be converted in degrees
  374.         final double lat = FastMath.toDegrees(latitude);
  375.         // Square root of the effective ionisation level
  376.         final double sqAz = FastMath.sqrt(az);
  377.         // seas parameter (Eq. 30 to 32)
  378.         final int seas;
  379.         if (month == 1 || month == 2 || month == 11 || month == 12) {
  380.             seas = -1;
  381.         } else if (month == 3 || month == 4 || month == 9 || month == 10) {
  382.             seas = 0;
  383.         } else {
  384.             seas = 1;
  385.         }
  386.         // Latitudinal dependence (Eq. 33 and 34)
  387.         final double ee = clipExp(0.3 * lat);
  388.         final double seasp = seas * ((ee - 1.0) / (ee + 1.0));
  389.         // Critical frequency (Eq. 35)
  390.         final double coef = 1.112 - 0.019 * seasp;
  391.         return FastMath.sqrt(coef * coef * sqAz * FastMath.pow(FastMath.cos(xeff), 0.6) + 0.49);

  393.     }

  395.     /**
  396.      * Computes the F2 layer height of maximum electron density.
  397.      * @param foE E layer layer critical frequency in MHz
  398.      * @param foF2 F2 layer layer critical frequency in MHz
  399.      * @param mF2 maximum usable frequency factor
  400.      * @return hmF2 in km
  401.      */
  402.     private double computehmF2(final double foE, final double foF2, final double mF2) {
  403.         // Ratio
  404.         final double fo = foF2 / foE;
  405.         final double ratio = join(fo, 1.75, 20.0, fo - 1.75);

  407.         // deltaM parameter
  408.         double deltaM = -0.012;
  409.         if (foE >= 1e-30) {
  410.             deltaM += 0.253 / (ratio - 1.215);
  411.         }

  413.         // hmF2 Eq. 80
  414.         final double mF2Sq = mF2 * mF2;
  415.         final double temp  = FastMath.sqrt((0.0196 * mF2Sq + 1) / (1.2967 * mF2Sq - 1.0));
  416.         return ((1490.0 * mF2 * temp) / (mF2 + deltaM)) - 176.0;

  418.     }

  420.     /** Compute sines and cosines.
  421.      * @param a argument
  422.      * @param n number of terms
  423.      * @return sin(a), cos(a), sin(2a), cos(2a) … sin(n a), cos(n a) array
  424.      * @since 12.1.3
  425.      */
  426.     private double[] sinCos(final double a, final int n) {

  428.         final SinCos sc0 = FastMath.sinCos(a);
  429.         SinCos sci = sc0;
  430.         final double[] sc = new double[2 * n];
  431.         int isc = 0;
  432.         sc[isc++] = sci.sin();
  433.         sc[isc++] = sci.cos();
  434.         for (int i = 1; i < n; i++) {
  435.             sci = SinCos.sum(sc0, sci);
  436.             sc[isc++] = sci.sin();
  437.             sc[isc++] = sci.cos();
  438.         }

  440.         return sc;

  442.     }

  444.     /**
  445.      * Computes cf2 coefficients.
  446.      * @param flattenF2 F2 coefficients used by the F2 layer (flatten array)
  447.      * @param azr effective sunspot number (Eq. 19)
  448.      * @param scT sines/cosines array of time argument
  449.      * @return the cf2 coefficients array
  450.      */
  451.     private double[] computeCF2(final double[] flattenF2, final double azr, final double[] scT) {

  453.         // interpolation coefficients for effective spot number
  454.         final double azr01 = azr * 0.01;
  455.         final double omazr01 = 1 - azr01;

  457.         // Eq. 44 and Eq. 50 merged into one loop
  458.         final double[] cf2 = new double[76];
  459.         int index = 0;
  460.         for (int i = 0; i < cf2.length; i++) {
  461.             cf2[i] = omazr01 * flattenF2[index     ] + azr01 * flattenF2[index +  1] +
  462.                     (omazr01 * flattenF2[index +  2] + azr01 * flattenF2[index +  3]) * scT[ 0] +
  463.                     (omazr01 * flattenF2[index +  4] + azr01 * flattenF2[index +  5]) * scT[ 1] +
  464.                     (omazr01 * flattenF2[index +  6] + azr01 * flattenF2[index +  7]) * scT[ 2] +
  465.                     (omazr01 * flattenF2[index +  8] + azr01 * flattenF2[index +  9]) * scT[ 3] +
  466.                     (omazr01 * flattenF2[index + 10] + azr01 * flattenF2[index + 11]) * scT[ 4] +
  467.                     (omazr01 * flattenF2[index + 12] + azr01 * flattenF2[index + 13]) * scT[ 5] +
  468.                     (omazr01 * flattenF2[index + 14] + azr01 * flattenF2[index + 15]) * scT[ 6] +
  469.                     (omazr01 * flattenF2[index + 16] + azr01 * flattenF2[index + 17]) * scT[ 7] +
  470.                     (omazr01 * flattenF2[index + 18] + azr01 * flattenF2[index + 19]) * scT[ 8] +
  471.                     (omazr01 * flattenF2[index + 20] + azr01 * flattenF2[index + 21]) * scT[ 9] +
  472.                     (omazr01 * flattenF2[index + 22] + azr01 * flattenF2[index + 23]) * scT[10] +
  473.                     (omazr01 * flattenF2[index + 24] + azr01 * flattenF2[index + 25]) * scT[11];
  474.             index += 26;
  475.         }
  476.         return cf2;
  477.     }

  479.     /**
  480.      * Computes Cm3 coefficients.
  481.      * @param flattenFm3 Fm3 coefficients used by the F2 layer (flatten array)
  482.      * @param azr effective sunspot number (Eq. 19)
  483.      * @param scT sines/cosines array of time argument
  484.      * @return the Cm3 coefficients array
  485.      */
  486.     private double[] computeCm3(final double[] flattenFm3, final double azr, final double[] scT) {

  488.         // interpolation coefficients for effective spot number
  489.         final double azr01 = azr * 0.01;
  490.         final double omazr01 = 1 - azr01;

  492.         // Eq. 44 and Eq. 51 merged into one loop
  493.         final double[] cm3 = new double[49];
  494.         int index = 0;
  495.         for (int i = 0; i < cm3.length; i++) {
  496.             cm3[i] = omazr01 * flattenFm3[index     ] + azr01 * flattenFm3[index +  1] +
  497.                     (omazr01 * flattenFm3[index +  2] + azr01 * flattenFm3[index +  3]) * scT[ 0] +
  498.                     (omazr01 * flattenFm3[index +  4] + azr01 * flattenFm3[index +  5]) * scT[ 1] +
  499.                     (omazr01 * flattenFm3[index +  6] + azr01 * flattenFm3[index +  7]) * scT[ 2] +
  500.                     (omazr01 * flattenFm3[index +  8] + azr01 * flattenFm3[index +  9]) * scT[ 3] +
  501.                     (omazr01 * flattenFm3[index + 10] + azr01 * flattenFm3[index + 11]) * scT[ 4] +
  502.                     (omazr01 * flattenFm3[index + 12] + azr01 * flattenFm3[index + 13]) * scT[ 5] +
  503.                     (omazr01 * flattenFm3[index + 14] + azr01 * flattenFm3[index + 15]) * scT[ 6] +
  504.                     (omazr01 * flattenFm3[index + 16] + azr01 * flattenFm3[index + 17]) * scT[ 7];
  505.             index += 18;
  506.         }
  507.         return cm3;
  508.     }

  510.     /**
  511.      * This method computes the F2 layer critical frequency.
  512.      * @param modip modified DIP latitude, in degrees
  513.      * @param cf2 Fourier time series for foF2
  514.      * @param latitude latitude in radians
  515.      * @param scL sines/cosines array of longitude argument
  516.      * @return the F2 layer critical frequency, MHz
  517.      */
  518.     private double computefoF2(final double modip, final double[] cf2,
  519.                                final double latitude, final double[] scL) {

  521.         // Legendre grades (Eq. 63)
  522.         final int[] q = new int[] {
  523.             12, 12, 9, 5, 2, 1, 1, 1, 1
  524.         };

  526.         double frequency = cf2[0];

  528.         // MODIP coefficients Eq. 57
  529.         final double sinMODIP = FastMath.sin(FastMath.toRadians(modip));
  530.         final double[] m = new double[12];
  531.         m[0] = 1.0;
  532.         for (int i = 1; i < q[0]; i++) {
  533.             m[i] = sinMODIP * m[i - 1];
  534.             frequency += m[i] * cf2[i];
  535.         }

  537.         // latitude and longitude terms
  538.         int index = 12;
  539.         final double cosLat1 = FastMath.cos(latitude);
  540.         double cosLatI = cosLat1;
  541.         for (int i = 1; i < q.length; i++) {
  542.             final double c = cosLatI * scL[2 * i - 1];
  543.             final double s = cosLatI * scL[2 * i - 2];
  544.             for (int j = 0; j < q[i]; j++) {
  545.                 frequency += m[j] * c * cf2[index++];
  546.                 frequency += m[j] * s * cf2[index++];
  547.             }
  548.             cosLatI *= cosLat1; // Eq. 58
  549.         }

  551.         return frequency;

  553.     }

  555.     /**
  556.      * This method computes the Maximum Usable Frequency factor.
  557.      * @param modip modified DIP latitude, in degrees
  558.      * @param cm3 Fourier time series for M(3000)F2
  559.      * @param latitude latitude in radians
  560.      * @param scL sines/cosines array of longitude argument
  561.      * @return the Maximum Usable Frequency factor
  562.      */
  563.     private double computeMF2(final double modip, final double[] cm3,
  564.                               final double latitude, final double[] scL) {

  566.         // Legendre grades (Eq. 71)
  567.         final int[] r = new int[] {
  568.             7, 8, 6, 3, 2, 1, 1
  569.         };

  571.         double m3000 = cm3[0];

  573.         // MODIP coefficients Eq. 57
  574.         final double sinMODIP = FastMath.sin(FastMath.toRadians(modip));
  575.         final double[] m = new double[12];
  576.         m[0] = 1.0;
  577.         for (int i = 1; i < 12; i++) {
  578.             m[i] = sinMODIP * m[i - 1];
  579.             if (i < 7) {
  580.                 m3000 += m[i] * cm3[i];
  581.             }
  582.         }

  584.         // latitude and longitude terms
  585.         int index = 7;
  586.         final double cosLat1 = FastMath.cos(latitude);
  587.         double cosLatI = cosLat1;
  588.         for (int i = 1; i < r.length; i++) {
  589.             final double c = cosLatI * scL[2 * i - 1];
  590.             final double s = cosLatI * scL[2 * i - 2];
  591.             for (int j = 0; j < r[i]; j++) {
  592.                 m3000 += m[j] * c * cm3[index++];
  593.                 m3000 += m[j] * s * cm3[index++];
  594.             }
  595.             cosLatI *= cosLat1; // Eq. 58
  596.         }

  598.         return m3000;

  600.     }

  602.     /**
  603.      * This method computes the F1 layer critical frequency.
  604.      * <p>
  605.      * This computation performs the algorithm exposed in Annex F
  606.      * of the reference document.
  607.      * </p>
  608.      * @param foE the E layer critical frequency, MHz
  609.      * @return the F1 layer critical frequency, MHz
  610.      * @param foF2 the F2 layer critical frequency, MHz
  611.      */
  612.     private double computefoF1(final double foE, final double foF2) {
  613.         final double temp  = join(1.4 * foE, 0.0, 1000.0, foE - 2.0);
  614.         final double temp2 = join(0.0, temp, 1000.0, foE - temp);
  615.         final double value = join(temp2, 0.85 * temp2, 60.0, 0.85 * foF2 - temp2);
  616.         if (value < 1.0E-6) {
  617.             return 0.0;
  618.         } else {
  619.             return value;
  620.         }
  621.     }

  623.     /**
  624.      * This method allows the computation of the F2, F1 and E layer amplitudes.
  625.      * <p>
  626.      * The resulting element is an array having the following form:
  627.      * <ul>
  628.      * <li>double[0] = A1 → F2 layer amplitude
  629.      * <li>double[1] = A2 → F1 layer amplitude
  630.      * <li>double[2] = A3 → E  layer amplitude
  631.      * </ul>
  632.      * </p>
  633.      * @param nmE E layer maximum density in 10^11 m-3
  634.      * @param nmF1 F1 layer maximum density in 10^11 m-3
  635.      * @param foF1 F1 layer critical frequency in MHz
  636.      * @return a three components array containing the layer amplitudes
  637.      */
  638.     private double[] computeLayerAmplitudes(final double nmE, final double nmF1, final double foF1) {
  639.         // Initialize array
  640.         final double[] amplitude = new double[3];

  642.         // F2 layer amplitude (Eq. 90)
  643.         final double a1 = 4.0 * nmF2;
  644.         amplitude[0]   = a1;

  646.         // F1 and E layer amplitudes (Eq. 91 to 98)
  647.         if (foF1 < 0.5) {
  648.             amplitude[1] = 0.0;
  649.             amplitude[2] = 4.0 * (nmE - epst(a1, hmF2, b2Bot, hmE));
  650.         } else {
  651.             double a2a = 0.0;
  652.             double a3a = 4.0 * nmE;
  653.             for (int i = 0; i < 5; i++) {
  654.                 a2a = 4.0 * (nmF1 - epst(a1, hmF2, b2Bot, hmF1) - epst(a3a, hmE, beTop, hmF1));
  655.                 a2a = join(a2a, 0.8 * nmF1, 1.0, a2a - 0.8 * nmF1);
  656.                 a3a = 4.0 * (nmE - epst(a2a, hmF1, b1Bot, hmE) - epst(a1, hmF2, b2Bot, hmE));
  657.             }
  658.             amplitude[1] = a2a;
  659.             amplitude[2] = join(a3a, 0.05, 60.0, a3a - 0.005);
  660.         }

  662.         return amplitude;
  663.     }

  665.     /**
  666.      * This method computes the topside thickness parameter H0.
  667.      *
  668.      * @param month current month
  669.      * @param azr effective sunspot number
  670.      * @return H0 in km
  671.      */
  672.     private double computeH0(final int month, final double azr) {

  674.         // Auxiliary parameter ka (Eq. 99 and 100)
  675.         final double ka;
  676.         if (month > 3 && month < 10) {
  677.             // month = 4,5,6,7,8,9
  678.             ka = 6.705 - 0.014 * azr - 0.008 * hmF2;
  679.         } else {
  680.             // month = 1,2,3,10,11,12
  681.             final double ratio = hmF2 / b2Bot;
  682.             ka = -7.77 + 0.097 * ratio * ratio + 0.153 * nmF2;
  683.         }

  685.         // Auxiliary parameter kb (Eq. 101 and 102)
  686.         double kb = join(ka, 2.0, 1.0, ka - 2.0);
  687.         kb = join(8.0, kb, 1.0, kb - 8.0);

  689.         // Auxiliary parameter Ha (Eq. 103)
  690.         final double hA = kb * b2Bot;

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

  696.         // Topside thickness parameter (Eq. 106)
  697.         return hA / v;

  699.     }

  701.     /**
  702.      * A clipped exponential function.
  703.      * <p>
  704.      * This function, describe in section F.2.12.2 of the reference document, is
  705.      * recommanded for the computation of exponential values.
  706.      * </p>
  707.      * @param power power for exponential function
  708.      * @return clipped exponential value
  709.      */
  710.     private double clipExp(final double power) {
  711.         if (power > 80.0) {
  712.             return 5.5406E34;
  713.         } else if (power < -80) {
  714.             return 1.8049E-35;
  715.         } else {
  716.             return FastMath.exp(power);
  717.         }
  718.     }

  720.     /**
  721.      * This method provides a third order interpolation function
  722.      * as recommended in the reference document (Ref Eq. 128 to Eq. 138)
  723.      *
  724.      * @param z1 z1 coefficient
  725.      * @param z2 z2 coefficient
  726.      * @param z3 z3 coefficient
  727.      * @param z4 z4 coefficient
  728.      * @param x position
  729.      * @return a third order interpolation
  730.      */
  731.     private double interpolate(final double z1, final double z2,
  732.                                final double z3, final double z4,
  733.                                final double x) {

  735.         if (FastMath.abs(2.0 * x) < 1e-10) {
  736.             return z2;
  737.         }

  739.         final double delta = 2.0 * x - 1.0;
  740.         final double g1 = z3 + z2;
  741.         final double g2 = z3 - z2;
  742.         final double g3 = z4 + z1;
  743.         final double g4 = (z4 - z1) / 3.0;
  744.         final double a0 = 9.0 * g1 - g3;
  745.         final double a1 = 9.0 * g2 - g4;
  746.         final double a2 = g3 - g1;
  747.         final double a3 = g4 - g2;
  748.         return 0.0625 * (a0 + delta * (a1 + delta * (a2 + delta * a3)));

  750.     }

  752.     /**
  753.      * Allows smooth joining of functions f1 and f2
  754.      * (i.e. continuous first derivatives) at origin.
  755.      * <p>
  756.      * This function, describe in section F.2.12.1 of the reference document, is
  757.      * recommanded for computational efficiency.
  758.      * </p>
  759.      * @param dF1 first function
  760.      * @param dF2 second function
  761.      * @param dA width of transition region
  762.      * @param dX x value
  763.      * @return the computed value
  764.      */
  765.     private double join(final double dF1, final double dF2,
  766.                         final double dA, final double dX) {
  767.         final double ee = clipExp(dA * dX);
  768.         return (dF1 * ee + dF2) / (ee + 1.0);
  769.     }

  771.     /**
  772.      * The Epstein function.
  773.      * <p>
  774.      * This function, describe in section 2.5.1 of the reference document, is used
  775.      * as a basis analytical function in NeQuick for the construction of the ionospheric layers.
  776.      * </p>
  777.      * @param x x parameter
  778.      * @param y y parameter
  779.      * @param z z parameter
  780.      * @param w w parameter
  781.      * @return value of the epstein function
  782.      */
  783.     private double epst(final double x, final double y,
  784.                         final double z, final double w) {
  785.         final double ex  = clipExp((w - y) / z);
  786.         final double opex = 1.0 + ex;
  787.         return x * ex / (opex * opex);

  789.     }

  791. }
Created with JaCoCo 0.8.12.202403310830