NeQuickModel.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 java.util.Collections;
import java.util.List;

import org.hipparchus.CalculusFieldElement;
import org.hipparchus.Field;
import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.MathArrays;
import org.orekit.bodies.BodyShape;
import org.orekit.bodies.FieldGeodeticPoint;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.frames.FieldStaticTransform;
import org.orekit.frames.Frame;
import org.orekit.frames.StaticTransform;
import org.orekit.frames.TopocentricFrame;
import org.orekit.models.earth.ionosphere.IonosphericDelayModel;
import org.orekit.models.earth.ionosphere.IonosphericModel;
import org.orekit.propagation.FieldSpacecraftState;
import org.orekit.propagation.SpacecraftState;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.DateComponents;
import org.orekit.time.DateTimeComponents;
import org.orekit.time.FieldAbsoluteDate;
import org.orekit.time.TimeScale;
import org.orekit.utils.ParameterDriver;
import org.orekit.utils.units.Unit;

/**
 * NeQuick ionospheric delay 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 abstract class NeQuickModel implements IonosphericModel, IonosphericDelayModel {

    /** Mean Earth radius in m (Ref Table 2.5.2). */
    public static final double RE = 6371200.0;

    /** Factor for the electron density computation. */
    private static final double DENSITY_FACTOR = 1.0e11;

    /** Factor for the path delay computation. */
    private static final double DELAY_FACTOR = 40.3e16;

    /** F2 coefficients used by the F2 layer (flatten array for cache efficiency). */
    private final double[][] flattenF2;

    /** Fm3 coefficients used by the M(3000)F2 layer(flatten array for cache efficiency). */
    private final double[][] flattenFm3;

    /** UTC time scale. */
    private final TimeScale utc;

    /** Simple constructor.
     * @param utc UTC time scale
     * @since 13.0
     */
    protected NeQuickModel(final TimeScale utc) {

        this.utc = utc;

        // F2 layer values
        this.flattenF2  = new double[12][];
        this.flattenFm3 = new double[12][];

    }

    /** Get UTC time scale.
     * @return UTC time scale
     * @since 13.0
     */
    public TimeScale getUtc() {
        return utc;
    }

    /** {@inheritDoc} */
    @Override
    public double pathDelay(final SpacecraftState state, final TopocentricFrame baseFrame,
                            final double frequency, final double[] parameters) {
        return pathDelay(state, baseFrame, state.getDate(), frequency, parameters);
    }

    /** {@inheritDoc} */
    @Override
    public double pathDelay(final SpacecraftState state,
                            final TopocentricFrame baseFrame, final AbsoluteDate receptionDate,
                            final double frequency, final double[] parameters) {

        // Reference body shape
        final BodyShape ellipsoid = baseFrame.getParentShape();

        // we use transform from body frame to inert frame and invert it
        // because base frame could be peered with inertial frame (hence improving performances with caching)
        // but the reverse is almost never used
        final Frame           bodyFrame  = ellipsoid.getBodyFrame();
        final StaticTransform body2Inert = bodyFrame.getStaticTransformTo(state.getFrame(), receptionDate);
        final Vector3D        position   = body2Inert.getInverse().transformPosition(state.getPosition());

        // Point
        final GeodeticPoint recPoint = baseFrame.getPoint();
        // Date
        final AbsoluteDate date = state.getDate();

        // Satellite geodetic coordinates
        final GeodeticPoint satPoint = ellipsoid.transform(position, bodyFrame, state.getDate());

        // Total Electron Content
        final double tec = stec(date, recPoint, satPoint);

        // Ionospheric delay
        final double factor = DELAY_FACTOR / (frequency * frequency);
        return factor * tec;
    }

    /** {@inheritDoc} */
    @Override
    public <T extends CalculusFieldElement<T>> T pathDelay(final FieldSpacecraftState<T> state,
                                                           final TopocentricFrame baseFrame,
                                                           final double frequency,
                                                           final T[] parameters) {
        return pathDelay(state, baseFrame, state.getDate(), frequency, parameters);
    }

    /** {@inheritDoc} */
    @Override
    public <T extends CalculusFieldElement<T>> T pathDelay(final FieldSpacecraftState<T> state,
                                                           final TopocentricFrame baseFrame,
                                                           final FieldAbsoluteDate<T> receptionDate,
                                                           final double frequency,
                                                           final T[] parameters) {
        // Reference body shape
        final BodyShape ellipsoid = baseFrame.getParentShape();

        // we use transform from body frame to inert frame and invert it
        // because base frame could be peered with inertial frame (hence improving performances with caching)
        // but the reverse is almost never used
        final Frame                   bodyFrame  = ellipsoid.getBodyFrame();
        final FieldStaticTransform<T> body2Inert = bodyFrame.getStaticTransformTo(state.getFrame(), receptionDate);
        final FieldVector3D<T>        position   = body2Inert.getInverse().transformPosition(state.getPosition());

        // Date
        final FieldAbsoluteDate<T> date = state.getDate();
        // Point
        final FieldGeodeticPoint<T> recPoint = baseFrame.getPoint(date.getField());

        // Satellite geodetic coordinates
        final FieldGeodeticPoint<T> satPoint = ellipsoid.transform(position, bodyFrame, state.getDate());

        // Total Electron Content
        final T tec = stec(date, recPoint, satPoint);

        // Ionospheric delay
        final double factor = DELAY_FACTOR / (frequency * frequency);
        return tec.multiply(factor);
    }

    /** {@inheritDoc} */
    @Override
    public List<ParameterDriver> getParametersDrivers() {
        return Collections.emptyList();
    }

    /**
     * This method allows the computation of the Slant Total Electron Content (STEC).
     * @param date current date
     * @param recP receiver position
     * @param satP satellite position
     * @return the STEC in TECUnits
     */
    public double stec(final AbsoluteDate date, final GeodeticPoint recP, final GeodeticPoint satP) {
        return stec(date.getComponents(utc), new Ray(recP, satP));
    }

    /**
     * This method allows the computation of the Slant Total Electron Content (STEC).
     * @param <T> type of the elements
     * @param date current date
     * @param recP receiver position
     * @param satP satellite position
     * @return the STEC in TECUnits
     */
    public <T extends CalculusFieldElement<T>> T stec(final FieldAbsoluteDate<T> date,
                                                      final FieldGeodeticPoint<T> recP,
                                                      final FieldGeodeticPoint<T> satP) {
        return stec(date.getComponents(utc), new FieldRay<>(recP, satP));
    }

    /** Compute modip for a location.
     * @param latitude latitude
     * @param longitude longitude
     * @return modip at specified location
     * @since 13.0
     */
    protected abstract double computeMODIP(double latitude, double longitude);

    /** Compute modip for a location.
     * @param <T> type of the field elements
     * @param latitude latitude
     * @param longitude longitude
     * @return modip at specified location
     * @since 13.0
     */
    protected abstract <T extends CalculusFieldElement<T>> T computeMODIP(T latitude, T longitude);

    /**
     * Compute Fourier time series.
     * @param dateTime current date time components
     * @param az effective ionisation level
     * @return Fourier time series
     * @since 13.0.1
     */
    public FourierTimeSeries computeFourierTimeSeries(final DateTimeComponents dateTime, final double az) {

         // Load the correct CCIR file
        loadsIfNeeded(dateTime.getDate());

        return new FourierTimeSeries(dateTime, az,
                                     flattenF2[dateTime.getDate().getMonth() - 1],
                                     flattenFm3[dateTime.getDate().getMonth() - 1]);

    }

    /**
     * Computes the electron density at a given height.
     * @param dateTime date
     * @param az effective ionization level
     * @param latitude latitude along the integration path
     * @param longitude longitude along the integration path
     * @param h height along the integration path in m
     * @return electron density [m⁻³]
     * @since 13.0
     */
    public double electronDensity(final DateTimeComponents dateTime, final double az,
                                  final double latitude, final double longitude, final double h) {
        return electronDensity(computeFourierTimeSeries(dateTime, az), latitude, longitude, h);
    }

    /**
     * Computes the electron density at a given height.
     * @param fourierTimeSeries Fourier time series for foF2 and M(3000)F2 layer (flatten array)
     * @param latitude latitude along the integration path
     * @param longitude longitude along the integration path
     * @param h height along the integration path in m
     * @return electron density [m⁻³]
     * @since 13.0.1
     */
    public double electronDensity(final FourierTimeSeries fourierTimeSeries,
                                  final double latitude, final double longitude, final double h) {

        final double modip = computeMODIP(latitude, longitude);
        final NeQuickParameters parameters = new NeQuickParameters(fourierTimeSeries, latitude, longitude, modip);

        // Convert height in kilometers
        final double hInKm = Unit.KILOMETRE.fromSI(h);

        // Electron density
        final double n;
        if (hInKm <= parameters.getHmF2()) {
            n = bottomElectronDensity(hInKm, parameters);
        } else {
            n = topElectronDensity(hInKm, parameters);
        }

        return n;

    }

    /**
     * Compute Fourier time series.
     * @param <T> type of the elements
     * @param dateTime current date time components
     * @param az effective ionisation level
     * @return Fourier time series
     * @since 13.0.1
     */
    public <T extends CalculusFieldElement<T>> FieldFourierTimeSeries<T> computeFourierTimeSeries(final DateTimeComponents dateTime,
                                                                                                  final T az) {

         // Load the correct CCIR file
        loadsIfNeeded(dateTime.getDate());

        return new FieldFourierTimeSeries<>(dateTime, az,
                                            flattenF2[dateTime.getDate().getMonth() - 1],
                                            flattenFm3[dateTime.getDate().getMonth() - 1]);

    }

    /**
     * Computes the electron density at a given height.
     * @param <T> type of the elements
     * @param dateTime date
     * @param az effective ionization level
     * @param latitude latitude along the integration path
     * @param longitude longitude along the integration path
     * @param h height along the integration path in m
     * @return electron density [m⁻³]
     * @since 13.0
     * CalculusFieldElement, CalculusFieldElement, CalculusFieldElement)}
     */
    public <T extends CalculusFieldElement<T>> T electronDensity(final DateTimeComponents dateTime, final T az,
                                                                 final T latitude, final T longitude, final T h) {
        return electronDensity(computeFourierTimeSeries(dateTime, az), latitude, longitude, h);
    }

    /**
     * Computes the electron density at a given height.
     * @param <T> type of the elements
     * @param fourierTimeSeries Fourier time series for foF2 and M(3000)F2 layer (flatten array)
     * @param latitude latitude along the integration path
     * @param longitude longitude along the integration path
     * @param h height along the integration path in m
     * @return electron density [m⁻³]
     * @since 13.0.1
     */
    public <T extends CalculusFieldElement<T>> T electronDensity(final FieldFourierTimeSeries<T> fourierTimeSeries,
                                                                 final T latitude, final T longitude, final T h) {

        final T modip = computeMODIP(latitude, longitude);
        final FieldNeQuickParameters<T> parameters = new FieldNeQuickParameters<>(fourierTimeSeries, latitude, longitude, modip);

        // Convert height in kilometers
        final T hInKm = Unit.KILOMETRE.fromSI(h);

        // Electron density
        final T n;
        if (hInKm.getReal() <= parameters.getHmF2().getReal()) {
            n = bottomElectronDensity(hInKm, parameters);
        } else {
            n = topElectronDensity(hInKm, parameters);
        }

        return n;

    }

    /**
     * Computes the electron density of the bottomside.
     * @param h height in km
     * @param parameters NeQuick model parameters
     * @return the electron density N in m⁻³
     */
    private double bottomElectronDensity(final double h, final NeQuickParameters parameters) {

        // Select the relevant B parameter for the current height (Eq. 109 and 110)
        final double be  = parameters.getBE(h);
        final double bf1 = parameters.getBF1(h);
        final double bf2 = parameters.getB2Bot();

        // Useful array of constants
        final double[] ct = new double[] {
            1.0 / bf2, 1.0 / bf1, 1.0 / be
        };

        // Compute the exponential argument for each layer (Eq. 111 to 113)
        final double[] arguments = computeExponentialArguments(h, parameters);

        // S coefficients
        final double[] s = new double[3];
        // Array of corrective terms
        final double[] ds = new double[3];

        // Layer amplitudes
        final double[] amplitudes = parameters.getLayerAmplitudes();

        // Fill arrays (Eq. 114 to 118)
        for (int i = 0; i < 3; i++) {
            if (FastMath.abs(arguments[i]) > 25.0) {
                s[i]  = 0.0;
                ds[i] = 0.0;
            } else {
                final double expA   = clipExp(arguments[i]);
                final double opExpA = 1.0 + expA;
                s[i]  = amplitudes[i] * (expA / (opExpA * opExpA));
                ds[i] = ct[i] * ((1.0 - expA) / (1.0 + expA));
            }
        }

        // Electron density
        final double aNo = s[0] + s[1] + s[2];
        if (applyChapmanParameters(h)) {
            // Chapman parameters (Eq. 119 and 120)
            final double bc = 1.0 - 10.0 * (MathArrays.linearCombination(s[0], ds[0], s[1], ds[1], s[2], ds[2]) / aNo);
            final double z  = 0.1 * (h - 100.0);
            // Electron density (Eq. 121)
            return aNo * clipExp(1.0 - bc * z - clipExp(-z)) * DENSITY_FACTOR;
        } else {
            return aNo * DENSITY_FACTOR;
        }

    }

    /**
     * Computes the electron density of the bottomside.
     * @param <T> type of the elements
     * @param h height in km
     * @param parameters NeQuick model parameters
     * @return the electron density N in m⁻³
     */
    private <T extends CalculusFieldElement<T>> T bottomElectronDensity(final T h,
                                                                        final FieldNeQuickParameters<T> parameters) {

        // Zero and One
        final Field<T> field = h.getField();
        final T zero = field.getZero();
        final T one  = field.getOne();

        // Select the relevant B parameter for the current height (Eq. 109 and 110)
        final T be  = parameters.getBE(h);
        final T bf1 = parameters.getBF1(h);
        final T bf2 = parameters.getB2Bot();

        // Useful array of constants
        final T[] ct = MathArrays.buildArray(field, 3);
        ct[0] = bf2.reciprocal();
        ct[1] = bf1.reciprocal();
        ct[2] = be.reciprocal();

        // Compute the exponential argument for each layer (Eq. 111 to 113)
        final T[] arguments = computeExponentialArguments(h, parameters);

        // S coefficients
        final T[] s = MathArrays.buildArray(field, 3);
        // Array of corrective terms
        final T[] ds = MathArrays.buildArray(field, 3);

        // Layer amplitudes
        final T[] amplitudes = parameters.getLayerAmplitudes();

        // Fill arrays (Eq. 114 to 118)
        for (int i = 0; i < 3; i++) {
            if (FastMath.abs(arguments[i]).getReal() > 25.0) {
                s[i]  = zero;
                ds[i] = zero;
            } else {
                final T expA   = clipExp(arguments[i]);
                final T opExpA = expA.add(1.0);
                s[i]  = amplitudes[i].multiply(expA.divide(opExpA.multiply(opExpA)));
                ds[i] = ct[i].multiply(expA.negate().add(1.0).divide(expA.add(1.0)));
            }
        }

        // Electron density
        final T aNo = s[0].add(s[1]).add(s[2]);
        if (applyChapmanParameters(h.getReal())) {
            // Chapman parameters (Eq. 119 and 120)
            final T bc = one.linearCombination(s[0], ds[0], s[1], ds[1], s[2], ds[2]).divide(aNo).multiply(10.0).negate().add(1.0);
            final T z  = h.subtract(100.0).multiply(0.1);
            // Electron density (Eq. 121)
            return aNo.multiply(clipExp(bc.multiply(z).add(clipExp(z.negate())).negate().add(1.0))).multiply(DENSITY_FACTOR);
        } else {
            return aNo.multiply(DENSITY_FACTOR);
        }

    }

    /**
     * Computes the electron density of the topside.
     * @param h height in km
     * @param parameters NeQuick model parameters
     * @return the electron density N in m⁻³
     */
    private double topElectronDensity(final double h, final NeQuickParameters parameters) {

        // Constant parameters (Eq. 122 and 123)
        final double g = 0.125;
        final double r = 100.0;

        // Arguments deltaH and z (Eq. 124 and 125)
        final double deltaH = h - parameters.getHmF2();
        final double h0     = computeH0(parameters);
        final double z      = deltaH / (h0 * (1.0 + (r * g * deltaH) / (r * h0 + g * deltaH)));

        // Exponential (Eq. 126)
        final double ee = clipExp(z);

        // Electron density (Eq. 127)
        if (ee > 1.0e11) {
            return (4.0 * parameters.getNmF2() / ee) * DENSITY_FACTOR;
        } else {
            final double opExpZ = 1.0 + ee;
            return ((4.0 * parameters.getNmF2() * ee) / (opExpZ * opExpZ)) * DENSITY_FACTOR;
        }
    }

    /**
     * Computes the electron density of the topside.
     * @param <T> type of the elements
     * @param h height in km
     * @param parameters NeQuick model parameters
     * @return the electron density N in m⁻³
     */
    private <T extends CalculusFieldElement<T>> T topElectronDensity(final T h,
                                                                     final FieldNeQuickParameters<T> parameters) {

        // Constant parameters (Eq. 122 and 123)
        final double g = 0.125;
        final double r = 100.0;

        // Arguments deltaH and z (Eq. 124 and 125)
        final T deltaH = h.subtract(parameters.getHmF2());
        final T h0     = computeH0(parameters);
        final T z      = deltaH.divide(h0.multiply(deltaH.multiply(r).multiply(g).divide(h0.multiply(r).add(deltaH.multiply(g))).add(1.0)));

        // Exponential (Eq. 126)
        final T ee = clipExp(z);

        // Electron density (Eq. 127)
        if (ee.getReal() > 1.0e11) {
            return parameters.getNmF2().multiply(4.0).divide(ee).multiply(DENSITY_FACTOR);
        } else {
            final T opExpZ = ee.add(1.0);
            return parameters.getNmF2().multiply(4.0).multiply(ee).divide(opExpZ.multiply(opExpZ)).multiply(DENSITY_FACTOR);
        }
    }

    /**
     * Lazy loading of CCIR data.
     * @param date current date components
     */
    private void loadsIfNeeded(final DateComponents date) {

        // Month index
        final int monthIndex = date.getMonth() - 1;

        // Check if CCIR has already been loaded for this month
        if (flattenF2[monthIndex] == null) {

            // Read file
            final CCIRLoader loader = new CCIRLoader();
            loader.loadCCIRCoefficients(date);

            // Update arrays
            this.flattenF2[monthIndex]  = flatten(loader.getF2());
            this.flattenFm3[monthIndex] = flatten(loader.getFm3());
        }
    }

    /** Flatten a 3-dimensions array.
     * <p>
     * This method convert 3-dimensions arrays into 1-dimension arrays
     * optimized to avoid cache misses when looping over all elements.
     * </p>
     * @param original original array a[i][j][k]
     * @return flatten array, for embedded loops on j (outer), k (intermediate), i (inner)
     */
    private double[] flatten(final double[][][] original) {
        final double[] flatten = new double[original.length * original[0].length * original[0][0].length];
        int index = 0;
        for (int j = 0; j < original[0].length; j++) {
            for (int k = 0; k < original[0][0].length; k++) {
                for (final double[][] doubles : original) {
                    flatten[index++] = doubles[j][k];
                }
            }
        }
        return flatten;
    }

    /**
     * A clipped exponential function.
     * <p>
     * This function, describe in section F.2.12.2 of the reference document, is
     * recommended for the computation of exponential values.
     * </p>
     * @param power power for exponential function
     * @return clipped exponential value
     */
    protected 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);
        }
    }

    /**
     * A clipped exponential function.
     * <p>
     * This function, describe in section F.2.12.2 of the reference document, is
     * recommended for the computation of exponential values.
     * </p>
     * @param <T> type of the elements
     * @param power power for exponential function
     * @return clipped exponential value
     */
    protected <T extends CalculusFieldElement<T>> T clipExp(final T power) {
        if (power.getReal() > 80.0) {
            return power.newInstance(5.5406E34);
        } else if (power.getReal() < -80) {
            return power.newInstance(1.8049E-35);
        } else {
            return FastMath.exp(power);
        }
    }

    /**
     * This method allows the computation of the Slant Total Electron Content (STEC).
     *
     * @param dateTime current date
     * @param ray      ray-perigee parameters
     * @return the STEC in TECUnits
     */
    abstract double stec(DateTimeComponents dateTime, Ray ray);

    /**
     * This method allows the computation of the Slant Total Electron Content (STEC).
     *
     * @param <T>      type of the field elements
     * @param dateTime current date
     * @param ray      ray-perigee parameters
     * @return the STEC in TECUnits
     */
    abstract <T extends CalculusFieldElement<T>> T stec(DateTimeComponents dateTime, FieldRay<T> ray);

    /**
     * Check if Chapman parameters should be applied.
     *
     * @param hInKm height in kilometers
     * @return true if Chapman parameters should be applied
     * @since 13.0
     */
    abstract boolean applyChapmanParameters(double hInKm);

    /**
     * Compute exponential arguments.
     * @param h height in km
     * @param parameters NeQuick model parameters
     * @return exponential arguments
     * @since 13.0
     */
    abstract double[] computeExponentialArguments(double h, NeQuickParameters parameters);

    /**
     * Compute exponential arguments.
     * @param <T>   type of the field elements
     * @param h height in km
     * @param parameters NeQuick model parameters
     * @return exponential arguments
     * @since 13.0
     */
    abstract <T extends CalculusFieldElement<T>> T[] computeExponentialArguments(T h,
                                                                                 FieldNeQuickParameters<T> parameters);

    /**
     * Compute topside thickness parameter.
     * @param parameters NeQuick model parameters
     * @return topside thickness parameter
     * @since 13.0
     */
    abstract double computeH0(NeQuickParameters parameters);

    /**
     * Compute topside thickness parameter.
     * @param <T>   type of the field elements
     * @param parameters NeQuick model parameters
     * @return topside thickness parameter
     * @since 13.0
     */
    abstract <T extends CalculusFieldElement<T>> T computeH0(FieldNeQuickParameters<T> parameters);

}