NeQuickParameters.java
/* Copyright 2002-2025 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.models.earth.ionosphere.nequick;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.SinCos;
import org.orekit.time.DateComponents;
import org.orekit.time.DateTimeComponents;
import org.orekit.time.TimeComponents;
/**
* This class performs the computation of the parameters used by the NeQuick model.
*
* @author Bryan Cazabonne
*
* @see "European Union (2016). European GNSS (Galileo) Open Service-Ionospheric Correction
* Algorithm for Galileo Single Frequency Users. 1.2."
* @see <a href="https://www.itu.int/rec/R-REC-P.531/en">ITU-R P.531</a>
*
* @since 10.1
*/
public class NeQuickParameters {
/** Solar zenith angle at day night transition, degrees. */
private static final double X0 = 86.23292796211615;
/** Current date time components.
* @since 13.0
*/
private final DateTimeComponents dateTime;
/** Effective sunspot number.
* @since 13.0
*/
private final double azr;
/** F2 layer critical frequency.
* @since 13.0
*/
private final double foF2;
/** F2 layer maximum density. */
private final double nmF2;
/** F2 layer maximum density height [km]. */
private final double hmF2;
/** F1 layer maximum density height [km]. */
private final double hmF1;
/** E layer maximum density height [km]. */
private final double hmE;
/** F2 layer bottom thickness parameter [km]. */
private final double b2Bot;
/** F1 layer top thickness parameter [km]. */
private final double b1Top;
/** F1 layer bottom thickness parameter [km]. */
private final double b1Bot;
/** E layer top thickness parameter [km]. */
private final double beTop;
/** E layer bottom thickness parameter [km]. */
private final double beBot;
/** Layer amplitudes. */
private final double[] amplitudes;
/**
* Build a new instance.
* @param dateTime current date time components
* @param flattenF2 F2 coefficients used by the F2 layer (flatten array)
* @param flattenFm3 Fm3 coefficients used by the M(3000)F2 layer (flatten array)
* @param latitude latitude of a point along the integration path, in radians
* @param longitude longitude of a point along the integration path, in radians
* @param az effective ionisation level
* @param modip modip
*/
public NeQuickParameters(final DateTimeComponents dateTime,
final double[] flattenF2, final double[] flattenFm3,
final double latitude, final double longitude, final double az, final double modip) {
this(new FourierTimeSeries(dateTime, az, flattenF2, flattenFm3),
latitude, longitude, modip);
}
/**
* Build a new instance.
* @param fourierTimeSeries Fourier time series for foF2 and M(3000)F2 layers
* @param latitude latitude of a point along the integration path, in radians
* @param longitude longitude of a point along the integration path, in radians
* @param modip modip
* @since 13.0.1
*/
public NeQuickParameters(final FourierTimeSeries fourierTimeSeries,
final double latitude, final double longitude, final double modip) {
this.dateTime = fourierTimeSeries.getDateTime();
this.azr = fourierTimeSeries.getAzr();
// Date and Time components
final DateComponents date = dateTime.getDate();
final TimeComponents time = dateTime.getTime();
// Hours
final double hours = time.getSecondsInUTCDay() / 3600.0;
// Effective solar zenith angle in radians
final double xeff = computeEffectiveSolarAngle(date.getMonth(), hours, latitude, longitude);
// E layer maximum density height in km (Eq. 78)
this.hmE = 120.0;
// E layer critical frequency in MHz
final double foE = computefoE(date.getMonth(), fourierTimeSeries.getAz(), xeff, latitude);
// E layer maximum density in 10^11 m⁻³ (Eq. 36)
final double nmE = 0.124 * foE * foE;
// F2 layer critical frequency in MHz
final double[] scL = FourierTimeSeries.sinCos(longitude, 8);
this.foF2 = computefoF2(modip, fourierTimeSeries.getCf2Reference(), latitude, scL);
// Maximum Usable Frequency factor
final double mF2 = computeMF2(modip, fourierTimeSeries.getCm3Reference(), latitude, scL);
// F2 layer maximum density in 10^11 m⁻³
this.nmF2 = 0.124 * foF2 * foF2;
// F2 layer maximum density height in km
this.hmF2 = computehmF2(foE, mF2);
// F1 layer critical frequency in MHz
final double foF1 = computefoF1(foE);
// F1 layer maximum density in 10^11 m⁻³
final double nmF1;
if (foF1 <= 0.0 && foE > 2.0) {
final double foEpopf = foE + 0.5;
nmF1 = 0.124 * foEpopf * foEpopf;
} else {
nmF1 = 0.124 * foF1 * foF1;
}
// F1 layer maximum density height in km
this.hmF1 = 0.5 * (hmF2 + hmE);
// Thickness parameters (Eq. 85 to 89)
final double a = 0.01 * clipExp(-3.467 + 0.857 * FastMath.log(foF2 * foF2) + 2.02 * FastMath.log(mF2));
this.b2Bot = 0.385 * nmF2 / a;
this.b1Top = 0.3 * (hmF2 - hmF1);
this.b1Bot = 0.5 * (hmF1 - hmE);
this.beTop = FastMath.max(b1Bot, 7.0);
this.beBot = 5.0;
// Layer amplitude coefficients
this.amplitudes = computeLayerAmplitudes(nmE, nmF1, foF1);
}
/**
* Get current date time components.
* @return current date time components
* @since 13.0
*/
public DateTimeComponents getDateTime() {
return dateTime;
}
/**
* Get effective sunspot number.
* @return effective sunspot number
* @since 13.0
*/
public double getAzr() {
return azr;
}
/**
* Get F2 layer critical frequency.
* @return F2 layer critical frequency
* @since 13.0
*/
public double getFoF2() {
return foF2;
}
/**
* Get the F2 layer maximum density.
* @return nmF2
*/
public double getNmF2() {
return nmF2;
}
/**
* Get the F2 layer maximum density height.
* @return hmF2 in km
*/
public double getHmF2() {
return hmF2;
}
/**
* Get the F1 layer maximum density height.
* @return hmF1 in km
*/
public double getHmF1() {
return hmF1;
}
/**
* Get the E layer maximum density height.
* @return hmE in km
*/
public double getHmE() {
return hmE;
}
/**
* Get the F2 layer thickness parameter (bottom).
* @return B2Bot in km
*/
public double getB2Bot() {
return b2Bot;
}
/**
* Get the F1 layer thickness parameter.
* @param h current height (km)
* @return B1 in km
* @since 13.0
*/
public double getBF1(final double h) {
// Eq. 110
return (h > hmF1) ? b1Top : b1Bot;
}
/**
* Get the E layer thickness parameter.
* @param h current height (km)
* @return Be in km
* @since 13.0
*/
public double getBE(final double h) {
// Eq. 109
return (h > hmE) ? beTop : beBot;
}
/**
* Get the F2, F1 and E layer amplitudes.
* <p>
* The resulting element is an array having the following form:
* <ul>
* <li>double[0] = A1 → F2 layer amplitude
* <li>double[1] = A2 → F1 layer amplitude
* <li>double[2] = A3 → E layer amplitude
* </ul>
* @return layer amplitudes
*/
public double[] getLayerAmplitudes() {
return amplitudes.clone();
}
/**
* This method computes the effective solar zenith angle.
* <p>
* The effective solar zenith angle is compute as a function of the
* solar zenith angle and the solar zenith angle at day night transition.
* </p>
* @param month current month of the year
* @param hours universal time (hours)
* @param latitude in radians
* @param longitude in radians
* @return the effective solar zenith angle, radians
*/
private double computeEffectiveSolarAngle(final int month,
final double hours,
final double latitude,
final double longitude) {
// Local time (Eq.4)
final double lt = hours + longitude / FastMath.toRadians(15.0);
// Day of year at the middle of the month (Eq. 20)
final double dy = 30.5 * month - 15.0;
// Time (Eq. 21)
final double t = dy + (18 - hours) / 24;
// Arguments am and al (Eq. 22 and 23)
final double am = FastMath.toRadians(0.9856 * t - 3.289);
final double al = am + FastMath.toRadians(1.916 * FastMath.sin(am) + 0.020 * FastMath.sin(2.0 * am) + 282.634);
// Sine and cosine of solar declination (Eq. 24 and 25)
final double sDec = 0.39782 * FastMath.sin(al);
final double cDec = FastMath.sqrt(1. - sDec * sDec);
// Solar zenith angle, deg (Eq. 26 and 27)
final SinCos scLat = FastMath.sinCos(latitude);
final double coef = (FastMath.PI / 12) * (12 - lt);
final double cZenith = scLat.sin() * sDec + scLat.cos() * cDec * FastMath.cos(coef);
final double angle = FastMath.atan2(FastMath.sqrt(1.0 - cZenith * cZenith), cZenith);
final double x = FastMath.toDegrees(angle);
// Effective solar zenith angle (Eq. 28)
final double xeff = join(90.0 - 0.24 * clipExp(20.0 - 0.2 * x), x, 12.0, x - X0);
return FastMath.toRadians(xeff);
}
/**
* This method computes the E layer critical frequency at a given location.
* @param month current month
* @param az ffective ionisation level
* @param xeff effective solar zenith angle in radians
* @param latitude latitude in radians
* @return the E layer critical frequency at a given location in MHz
*/
private double computefoE(final int month, final double az,
final double xeff, final double latitude) {
// The latitude has to be converted in degrees
final double lat = FastMath.toDegrees(latitude);
// Square root of the effective ionisation level
final double sqAz = FastMath.sqrt(az);
// seas parameter (Eq. 30 to 32)
final int seas;
if (month == 1 || month == 2 || month == 11 || month == 12) {
seas = -1;
} else if (month == 3 || month == 4 || month == 9 || month == 10) {
seas = 0;
} else {
seas = 1;
}
// Latitudinal dependence (Eq. 33 and 34)
final double ee = clipExp(0.3 * lat);
final double seasp = seas * ((ee - 1.0) / (ee + 1.0));
// Critical frequency (Eq. 35)
final double coef = 1.112 - 0.019 * seasp;
return FastMath.sqrt(coef * coef * sqAz * FastMath.pow(FastMath.cos(xeff), 0.6) + 0.49);
}
/**
* Computes the F2 layer height of maximum electron density.
* @param foE E layer layer critical frequency in MHz
* @param mF2 maximum usable frequency factor
* @return hmF2 in km
*/
private double computehmF2(final double foE, final double mF2) {
// Ratio
final double fo = foF2 / foE;
final double ratio = join(fo, 1.75, 20.0, fo - 1.75);
// deltaM parameter
double deltaM = -0.012;
if (foE >= 1e-30) {
deltaM += 0.253 / (ratio - 1.215);
}
// hmF2 Eq. 80
final double mF2Sq = mF2 * mF2;
final double temp = FastMath.sqrt((0.0196 * mF2Sq + 1) / (1.2967 * mF2Sq - 1.0));
return ((1490.0 * mF2 * temp) / (mF2 + deltaM)) - 176.0;
}
/**
* This method computes the F2 layer critical frequency.
* @param modip modified DIP latitude, in degrees
* @param cf2 Fourier time series for foF2
* @param latitude latitude in radians
* @param scL sines/cosines array of longitude argument
* @return the F2 layer critical frequency, MHz
*/
private double computefoF2(final double modip, final double[] cf2,
final double latitude, final double[] scL) {
// Legendre grades (Eq. 63)
final int[] q = new int[] {
12, 12, 9, 5, 2, 1, 1, 1, 1
};
double frequency = cf2[0];
// ModipGrid coefficients Eq. 57
final double sinMODIP = FastMath.sin(FastMath.toRadians(modip));
final double[] m = new double[12];
m[0] = 1.0;
for (int i = 1; i < q[0]; i++) {
m[i] = sinMODIP * m[i - 1];
frequency += m[i] * cf2[i];
}
// latitude and longitude terms
int index = 12;
final double cosLat1 = FastMath.cos(latitude);
double cosLatI = cosLat1;
for (int i = 1; i < q.length; i++) {
final double c = cosLatI * scL[2 * i - 1];
final double s = cosLatI * scL[2 * i - 2];
for (int j = 0; j < q[i]; j++) {
frequency += m[j] * c * cf2[index++];
frequency += m[j] * s * cf2[index++];
}
cosLatI *= cosLat1; // Eq. 58
}
return frequency;
}
/**
* This method computes the Maximum Usable Frequency factor.
* @param modip modified DIP latitude, in degrees
* @param cm3 Fourier time series for M(3000)F2
* @param latitude latitude in radians
* @param scL sines/cosines array of longitude argument
* @return the Maximum Usable Frequency factor
*/
private double computeMF2(final double modip, final double[] cm3,
final double latitude, final double[] scL) {
// Legendre grades (Eq. 71)
final int[] r = new int[] {
7, 8, 6, 3, 2, 1, 1
};
double m3000 = cm3[0];
// ModipGrid coefficients Eq. 57
final double sinMODIP = FastMath.sin(FastMath.toRadians(modip));
final double[] m = new double[12];
m[0] = 1.0;
for (int i = 1; i < 12; i++) {
m[i] = sinMODIP * m[i - 1];
if (i < 7) {
m3000 += m[i] * cm3[i];
}
}
// latitude and longitude terms
int index = 7;
final double cosLat1 = FastMath.cos(latitude);
double cosLatI = cosLat1;
for (int i = 1; i < r.length; i++) {
final double c = cosLatI * scL[2 * i - 1];
final double s = cosLatI * scL[2 * i - 2];
for (int j = 0; j < r[i]; j++) {
m3000 += m[j] * c * cm3[index++];
m3000 += m[j] * s * cm3[index++];
}
cosLatI *= cosLat1; // Eq. 58
}
return m3000;
}
/**
* This method computes the F1 layer critical frequency.
* <p>
* This computation performs the algorithm exposed in Annex F
* of the reference document.
* </p>
* @param foE the E layer critical frequency, MHz
* @return the F1 layer critical frequency, MHz
*/
private double computefoF1(final double foE) {
final double temp = join(1.4 * foE, 0.0, 1000.0, foE - 2.0);
final double temp2 = join(0.0, temp, 1000.0, foE - temp);
final double value = join(temp2, 0.85 * temp2, 60.0, 0.85 * foF2 - temp2);
if (value < 1.0E-6) {
return 0.0;
} else {
return value;
}
}
/**
* This method allows the computation of the F2, F1 and E layer amplitudes.
* <p>
* The resulting element is an array having the following form:
* <ul>
* <li>double[0] = A1 → F2 layer amplitude
* <li>double[1] = A2 → F1 layer amplitude
* <li>double[2] = A3 → E layer amplitude
* </ul>
* </p>
* @param nmE E layer maximum density in 10^11 m⁻³
* @param nmF1 F1 layer maximum density in 10^11 m⁻³
* @param foF1 F1 layer critical frequency in MHz
* @return a three components array containing the layer amplitudes
*/
private double[] computeLayerAmplitudes(final double nmE, final double nmF1, final double foF1) {
// Initialize array
final double[] amplitude = new double[3];
// F2 layer amplitude (Eq. 90)
final double a1 = 4.0 * nmF2;
amplitude[0] = a1;
// F1 and E layer amplitudes (Eq. 91 to 98)
if (foF1 < 0.5) {
amplitude[1] = 0.0;
amplitude[2] = 4.0 * (nmE - epst(a1, hmF2, b2Bot, hmE));
} else {
double a2a = 0.0;
double a3a = 4.0 * nmE;
for (int i = 0; i < 5; i++) {
a2a = 4.0 * (nmF1 - epst(a1, hmF2, b2Bot, hmF1) - epst(a3a, hmE, beTop, hmF1));
a2a = join(a2a, 0.8 * nmF1, 1.0, a2a - 0.8 * nmF1);
a3a = 4.0 * (nmE - epst(a2a, hmF1, b1Bot, hmE) - epst(a1, hmF2, b2Bot, hmE));
}
amplitude[1] = a2a;
amplitude[2] = join(a3a, 0.05, 60.0, a3a - 0.005);
}
return amplitude;
}
/**
* A clipped exponential function.
* <p>
* This function, describe in section F.2.12.2 of the reference document, is
* recommanded for the computation of exponential values.
* </p>
* @param power power for exponential function
* @return clipped exponential value
*/
private double clipExp(final double power) {
if (power > 80.0) {
return 5.5406E34;
} else if (power < -80) {
return 1.8049E-35;
} else {
return FastMath.exp(power);
}
}
/**
* Allows smooth joining of functions f1 and f2
* (i.e. continuous first derivatives) at origin.
* <p>
* This function, describe in section F.2.12.1 of the reference document, is
* recommended for computational efficiency.
* </p>
* @param dF1 first function
* @param dF2 second function
* @param dA width of transition region
* @param dX x value
* @return the computed value
*/
double join(final double dF1, final double dF2, final double dA, final double dX) {
final double ee = clipExp(dA * dX);
return (dF1 * ee + dF2) / (ee + 1.0);
}
/**
* The Epstein function.
* <p>
* This function, describe in section 2.5.1 of the reference document, is used
* as a basis analytical function in NeQuick for the construction of the ionospheric layers.
* </p>
* @param x x parameter
* @param y y parameter
* @param z z parameter
* @param w w parameter
* @return value of the epstein function
*/
private double epst(final double x, final double y,
final double z, final double w) {
final double ex = clipExp((w - y) / z);
final double opex = 1.0 + ex;
return x * ex / (opex * opex);
}
}