NeQuickModel.java

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

import java.io.BufferedReader;
import java.io.IOException;
import java.io.InputStream;
import java.io.InputStreamReader;
import java.nio.charset.StandardCharsets;
import java.util.Collections;
import java.util.List;
import java.util.regex.Pattern;

import org.hipparchus.CalculusFieldElement;
import org.hipparchus.Field;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.FieldSinCos;
import org.hipparchus.util.MathArrays;
import org.hipparchus.util.SinCos;
import org.orekit.annotation.DefaultDataContext;
import org.orekit.bodies.BodyShape;
import org.orekit.bodies.FieldGeodeticPoint;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.data.DataContext;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitMessages;
import org.orekit.frames.TopocentricFrame;
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;

/**
 * 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."
 *
 * @since 10.1
 */
public class NeQuickModel implements IonosphericModel {

    /** NeQuick resources base directory. */
    private static final String NEQUICK_BASE = "/assets/org/orekit/nequick/";

    /** Pattern for delimiting regular expressions. */
    private static final Pattern SEPARATOR = Pattern.compile("\\s+");

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

    /** Meters to kilometers converter. */
    private static final double M_TO_KM = 0.001;

    /** 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;

    /** The three ionospheric coefficients broadcast in the Galileo navigation message. */
    private final double[] alpha;

    /** MODIP grid. */
    private final double[][] stModip;

    /** Month used for loading CCIR coefficients. */
    private int month;

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

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

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

    /**
     * Build a new instance.
     *
     * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
     *
     * @param alpha effective ionisation level coefficients
     * @see #NeQuickModel(double[], TimeScale)
     */
    @DefaultDataContext
    public NeQuickModel(final double[] alpha) {
        this(alpha, DataContext.getDefault().getTimeScales().getUTC());
    }

    /**
     * Build a new instance.
     * @param alpha effective ionisation level coefficients
     * @param utc UTC time scale.
     * @since 10.1
     */
    public NeQuickModel(final double[] alpha,
                        final TimeScale utc) {
        // F2 layer values
        this.month      = 0;
        this.flattenF2  = null;
        this.flattenFm3 = null;
        // Read modip grid
        final MODIPLoader parser = new MODIPLoader();
        parser.loadMODIPGrid();
        this.stModip = parser.getMODIPGrid();
        // Ionisation level coefficients
        this.alpha = alpha.clone();
        this.utc = utc;
    }

    @Override
    public double pathDelay(final SpacecraftState state, final TopocentricFrame baseFrame,
                            final double frequency, final double[] parameters) {
        // Point
        final GeodeticPoint recPoint = baseFrame.getPoint();
        // Date
        final AbsoluteDate date = state.getDate();

        // Reference body shape
        final BodyShape ellipsoid = baseFrame.getParentShape();
        // Satellite geodetic coordinates
        final GeodeticPoint satPoint = ellipsoid.transform(state.getPosition(ellipsoid.getBodyFrame()), ellipsoid.getBodyFrame(), state.getDate());

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

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

    @Override
    public <T extends CalculusFieldElement<T>> T pathDelay(final FieldSpacecraftState<T> state, final TopocentricFrame baseFrame,
                                                       final double frequency, final T[] parameters) {
        // Date
        final FieldAbsoluteDate<T> date = state.getDate();
        // Point
        final FieldGeodeticPoint<T> recPoint = baseFrame.getPoint(date.getField());


        // Reference body shape
        final BodyShape ellipsoid = baseFrame.getParentShape();
        // Satellite geodetic coordinates
        final FieldGeodeticPoint<T> satPoint = ellipsoid.transform(state.getPosition(ellipsoid.getBodyFrame()), ellipsoid.getBodyFrame(), 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);
    }

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

    /**
     * This method allows the computation of the Stant Total Electron Content (STEC).
     * <p>
     * This method follows the Gauss algorithm exposed in section 2.5.8.2.8 of
     * the reference document.
     * </p>
     * @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) {

        // Ray-perigee parameters
        final Ray ray = new Ray(recP, satP);

        // Load the correct CCIR file
        final DateTimeComponents dateTime = date.getComponents(utc);
        loadsIfNeeded(dateTime.getDate());

        // Tolerance for the integration accuracy. Defined inside the reference document, section 2.5.8.1.
        final double h1 = recP.getAltitude();
        final double tolerance;
        if (h1 < 1000000.0) {
            tolerance = 0.001;
        } else {
            tolerance = 0.01;
        }

        // Integration
        int n = 8;
        final Segment seg1 = new Segment(n, ray);
        double gn1 = stecIntegration(seg1, dateTime);
        n *= 2;
        final Segment seg2 = new Segment(n, ray);
        double gn2 = stecIntegration(seg2, dateTime);

        int count = 1;
        while (FastMath.abs(gn2 - gn1) > tolerance * FastMath.abs(gn1) && count < 20) {
            gn1 = gn2;
            n *= 2;
            final Segment seg = new Segment(n, ray);
            gn2 = stecIntegration(seg, dateTime);
            count += 1;
        }

        // If count > 20 the integration did not converge
        if (count == 20) {
            throw new OrekitException(OrekitMessages.STEC_INTEGRATION_DID_NOT_CONVERGE);
        }

        // Eq. 202
        return (gn2 + ((gn2 - gn1) / 15.0)) * 1.0e-16;
    }

    /**
     * This method allows the computation of the Stant Total Electron Content (STEC).
     * <p>
     * This method follows the Gauss algorithm exposed in section 2.5.8.2.8 of
     * the reference document.
     * </p>
     * @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) {

        // Field
        final Field<T> field = date.getField();

        // Ray-perigee parameters
        final FieldRay<T> ray = new FieldRay<>(field, recP, satP);

        // Load the correct CCIR file
        final DateTimeComponents dateTime = date.getComponents(utc);
        loadsIfNeeded(dateTime.getDate());

        // Tolerance for the integration accuracy. Defined inside the reference document, section 2.5.8.1.
        final T h1 = recP.getAltitude();
        final double tolerance;
        if (h1.getReal() < 1000000.0) {
            tolerance = 0.001;
        } else {
            tolerance = 0.01;
        }

        // Integration
        int n = 8;
        final FieldSegment<T> seg1 = new FieldSegment<>(field, n, ray);
        T gn1 = stecIntegration(field, seg1, dateTime);
        n *= 2;
        final FieldSegment<T> seg2 = new FieldSegment<>(field, n, ray);
        T gn2 = stecIntegration(field, seg2, dateTime);

        int count = 1;
        while (FastMath.abs(gn2.subtract(gn1)).getReal() > FastMath.abs(gn1).multiply(tolerance).getReal() && count < 20) {
            gn1 = gn2;
            n *= 2;
            final FieldSegment<T> seg = new FieldSegment<>(field, n, ray);
            gn2 = stecIntegration(field, seg, dateTime);
            count += 1;
        }

        // If count > 20 the integration did not converge
        if (count == 20) {
            throw new OrekitException(OrekitMessages.STEC_INTEGRATION_DID_NOT_CONVERGE);
        }

        // Eq. 202
        return gn2.add(gn2.subtract(gn1).divide(15.0)).multiply(1.0e-16);
    }

    /**
     * This method perfoms the STEC integration.
     * @param seg coordinates along the integration path
     * @param dateTime current date and time componentns
     * @return result of the integration
     */
    private double stecIntegration(final Segment seg, final DateTimeComponents dateTime) {
        // Integration points
        final double[] heightS    = seg.getHeights();
        final double[] latitudeS  = seg.getLatitudes();
        final double[] longitudeS = seg.getLongitudes();

        // Compute electron density
        double density = 0.0;
        for (int i = 0; i < heightS.length; i++) {
            final NeQuickParameters parameters = new NeQuickParameters(dateTime, flattenF2, flattenFm3,
                                                                       latitudeS[i], longitudeS[i],
                                                                       alpha, stModip);
            density += electronDensity(heightS[i], parameters);
        }

        return 0.5 * seg.getInterval() * density;
    }

    /**
     * This method perfoms the STEC integration.
     * @param <T> type of the elements
     * @param field field of the elements
     * @param seg coordinates along the integration path
     * @param dateTime current date and time componentns
     * @return result of the integration
     */
    private <T extends CalculusFieldElement<T>> T stecIntegration(final Field<T> field,
                                                                  final FieldSegment<T> seg,
                                                                  final DateTimeComponents dateTime) {
        // Integration points
        final T[] heightS    = seg.getHeights();
        final T[] latitudeS  = seg.getLatitudes();
        final T[] longitudeS = seg.getLongitudes();

        // Compute electron density
        T density = field.getZero();
        for (int i = 0; i < heightS.length; i++) {
            final FieldNeQuickParameters<T> parameters = new FieldNeQuickParameters<>(dateTime, flattenF2, flattenFm3,
                                                                                      latitudeS[i], longitudeS[i],
                                                                                      alpha, stModip);
            density = density.add(electronDensity(field, heightS[i], parameters));
        }

        return seg.getInterval().multiply(density).multiply(0.5);
    }

    /**
     * Computes the electron density at a given height.
     * @param h height in m
     * @param parameters NeQuick model parameters
     * @return electron density [m^-3]
     */
    private double electronDensity(final double h, final NeQuickParameters parameters) {
        // Convert height in kilometers
        final double hInKm = h * M_TO_KM;
        // Electron density
        final double n;
        if (hInKm <= parameters.getHmF2()) {
            n = bottomElectronDensity(hInKm, parameters);
        } else {
            n = topElectronDensity(hInKm, parameters);
        }
        return n;
    }

    /**
     * Computes the electron density at a given height.
     * @param <T> type of the elements
     * @param field field of the elements
     * @param h height in m
     * @param parameters NeQuick model parameters
     * @return electron density [m^-3]
     */
    private <T extends CalculusFieldElement<T>> T electronDensity(final Field<T> field,
                                                              final T h,
                                                              final FieldNeQuickParameters<T> parameters) {
        // Convert height in kilometers
        final T hInKm = h.multiply(M_TO_KM);
        // Electron density
        final T n;
        if (hInKm.getReal() <= parameters.getHmF2().getReal()) {
            n = bottomElectronDensity(field, hInKm, parameters);
        } else {
            n = topElectronDensity(field, 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-3
     */
    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;
        if (h > parameters.getHmE()) {
            be = parameters.getBETop();
        } else {
            be = parameters.getBEBot();
        }
        final double bf1;
        if (h > parameters.getHmF1()) {
            bf1 = parameters.getB1Top();
        } else {
            bf1 = parameters.getB1Bot();
        }
        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)
        // If h < 100km we use h = 100km as recommended in the reference document
        final double   hTemp = FastMath.max(100.0, h);
        final double   exp   = clipExp(10.0 / (1.0 + FastMath.abs(hTemp - parameters.getHmF2())));
        final double[] arguments = new double[3];
        arguments[0] = (hTemp - parameters.getHmF2()) / bf2;
        arguments[1] = ((hTemp - parameters.getHmF1()) / bf1) * exp;
        arguments[2] = ((hTemp - parameters.getHmE()) / be) * exp;

        // 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 = MathArrays.linearCombination(s[0], 1.0, s[1], 1.0, s[2], 1.0);
        if (h >= 100) {
            return aNo * DENSITY_FACTOR;
        } else {
            // 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;
        }
    }

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

        // Zero and One
        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;
        if (h.getReal() > parameters.getHmE().getReal()) {
            be = parameters.getBETop();
        } else {
            be = parameters.getBEBot();
        }
        final T bf1;
        if (h.getReal() > parameters.getHmF1().getReal()) {
            bf1 = parameters.getB1Top();
        } else {
            bf1 = parameters.getB1Bot();
        }
        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)
        // If h < 100km we use h = 100km as recommended in the reference document
        final T   hTemp = FastMath.max(zero.newInstance(100.0), h);
        final T   exp   = clipExp(field, FastMath.abs(hTemp.subtract(parameters.getHmF2())).add(1.0).divide(10.0).reciprocal());
        final T[] arguments = MathArrays.buildArray(field, 3);
        arguments[0] = hTemp.subtract(parameters.getHmF2()).divide(bf2);
        arguments[1] = hTemp.subtract(parameters.getHmF1()).divide(bf1).multiply(exp);
        arguments[2] = hTemp.subtract(parameters.getHmE()).divide(be).multiply(exp);

        // 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(field, 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 = one.linearCombination(s[0], one, s[1], one, s[2], one);
        if (h.getReal() >= 100) {
            return aNo.multiply(DENSITY_FACTOR);
        } else {
            // Chapman parameters (Eq. 119 and 120)
            final T bc = s[0].multiply(ds[0]).add(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(field, bc.multiply(z).add(clipExp(field, z.negate())).negate().add(1.0))).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-3
     */
    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 z      = deltaH / (parameters.getH0() * (1.0 + (r * g * deltaH) / (r * parameters.getH0() + 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 field field of the elements
     * @param h height in km
     * @param parameters NeQuick model parameters
     * @return the electron density N in m-3
     */
    private <T extends CalculusFieldElement<T>> T topElectronDensity(final Field<T> field,
                                                                 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 z      = deltaH.divide(parameters.getH0().multiply(deltaH.multiply(r).multiply(g).divide(parameters.getH0().multiply(r).add(deltaH.multiply(g))).add(1.0)));

        // Exponential (Eq. 126)
        final T ee = clipExp(field, 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(field.getOne());
            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) {

        // Current month
        final int currentMonth = date.getMonth();

        // Check if date have changed or if f2 and fm3 arrays are null
        if (currentMonth != month || flattenF2 == null || flattenFm3 == null) {
            this.month = currentMonth;

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

            // Update arrays
            this.flattenF2  = flatten(loader.getF2());
            this.flattenFm3 = 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
     * 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);
        }
    }

    /**
     * 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 <T> type of the elements
     * @param field field of the elements
     * @param power power for exponential function
     * @return clipped exponential value
     */
    private <T extends CalculusFieldElement<T>> T clipExp(final Field<T> field, final T power) {
        final T zero = field.getZero();
        if (power.getReal() > 80.0) {
            return zero.newInstance(5.5406E34);
        } else if (power.getReal() < -80) {
            return zero.newInstance(1.8049E-35);
        } else {
            return FastMath.exp(power);
        }
    }

    /** Get a data stream.
     * @param name file name of the resource stream
     * @return stream
     */
    private static InputStream getStream(final String name) {
        return NeQuickModel.class.getResourceAsStream(name);
    }

    /**
     * Parser for Modified Dip Latitude (MODIP) grid file.
     * <p>
     * The MODIP grid allows to estimate MODIP μ [deg] at a given point (φ,λ)
     * by interpolation of the relevant values contained in the support file.
     * </p> <p>
     * The file contains the values of MODIP (expressed in degrees) on a geocentric grid
     * from 90°S to 90°N with a 5-degree step in latitude and from 180°W to 180°E with a
     * 10-degree in longitude.
     * </p>
     */
    private static class MODIPLoader {

        /** Supported name for MODIP grid. */
        private static final String SUPPORTED_NAME = NEQUICK_BASE + "modip.txt";

        /** MODIP grid. */
        private double[][] grid;

        /**
         * Build a new instance.
         */
        MODIPLoader() {
            this.grid = null;
        }

        /** Returns the MODIP grid array.
         * @return the MODIP grid array
         */
        public double[][] getMODIPGrid() {
            return grid.clone();
        }

        /**
         * Load the data using supported names.
         */
        public void loadMODIPGrid() {
            try (InputStream in = getStream(SUPPORTED_NAME)) {
                loadData(in, SUPPORTED_NAME);
            } catch (IOException e) {
                throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
            }

            // Throw an exception if MODIP grid was not loaded properly
            if (grid == null) {
                throw new OrekitException(OrekitMessages.MODIP_GRID_NOT_LOADED, SUPPORTED_NAME);
            }
        }

        /**
         * Load data from a stream.
         * @param input input stream
         * @param name name of the file
         * @throws IOException if data can't be read
         */
        public void loadData(final InputStream input, final String name)
            throws IOException {

            // Grid size
            final int size = 39;

            // Initialize array
            final double[][] array = new double[size][size];

            // Open stream and parse data
            int   lineNumber = 0;
            String line      = null;
            try (InputStreamReader isr = new InputStreamReader(input, StandardCharsets.UTF_8);
                 BufferedReader    br = new BufferedReader(isr)) {

                for (line = br.readLine(); line != null; line = br.readLine()) {
                    ++lineNumber;
                    line = line.trim();

                    // Read grid data
                    if (!line.isEmpty()) {
                        final String[] modip_line = SEPARATOR.split(line);
                        for (int column = 0; column < modip_line.length; column++) {
                            array[lineNumber - 1][column] = Double.parseDouble(modip_line[column]);
                        }
                    }

                }

            } catch (NumberFormatException nfe) {
                throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
                                          lineNumber, name, line);
            }

            // Clone parsed grid
            grid = array.clone();

        }
    }

    /**
     * Parser for CCIR files.
     * <p>
     * Numerical grid maps which describe the regular variation of the ionosphere.
     * They are used to derive other variables such as critical frequencies and transmission factors.
     * </p> <p>
     * The coefficients correspond to low and high solar activity conditions.
     * </p> <p>
     * The CCIR file naming convention is ccirXX.asc where each XX means month + 10.
     * </p> <p>
     * Coefficients are store into tow arrays, F2 and Fm3. F2 coefficients are used for the computation
     * of the F2 layer critical frequency. Fm3 for the computation of the F2 layer maximum usable frequency factor.
     * The size of these two arrays is fixed and discussed into the section 2.5.3.2
     * of the reference document.
     * </p>
     */
    private static class CCIRLoader {

        /** Default supported files name pattern. */
        public static final String DEFAULT_SUPPORTED_NAME = "ccir**.asc";

        /** Total number of F2 coefficients contained in the file. */
        private static final int NUMBER_F2_COEFFICIENTS = 1976;

        /** Rows number for F2 and Fm3 arrays. */
        private static final int ROWS = 2;

        /** Columns number for F2 array. */
        private static final int TOTAL_COLUMNS_F2 = 76;

        /** Columns number for Fm3 array. */
        private static final int TOTAL_COLUMNS_FM3 = 49;

        /** Depth of F2 array. */
        private static final int DEPTH_F2 = 13;

        /** Depth of Fm3 array. */
        private static final int DEPTH_FM3 = 9;

        /** Regular expression for supported file name. */
        private String supportedName;

        /** F2 coefficients used for the computation of the F2 layer critical frequency. */
        private double[][][] f2Loader;

        /** Fm3 coefficients used for the computation of the F2 layer maximum usable frequency factor. */
        private double[][][] fm3Loader;

        /**
         * Build a new instance.
         */
        CCIRLoader() {
            this.supportedName = DEFAULT_SUPPORTED_NAME;
            this.f2Loader  = null;
            this.fm3Loader = null;
        }

        /**
         * Get the F2 coefficients used for the computation of the F2 layer critical frequency.
         * @return the F2 coefficients
         */
        public double[][][] getF2() {
            return f2Loader.clone();
        }

        /**
         * Get the Fm3 coefficients used for the computation of the F2 layer maximum usable frequency factor.
         * @return the F2 coefficients
         */
        public double[][][] getFm3() {
            return fm3Loader.clone();
        }

        /** Load the data for a given month.
         * @param dateComponents month given but its DateComponents
         */
        public void loadCCIRCoefficients(final DateComponents dateComponents) {

            // The files are named ccirXX.asc where XX substitute the month of the year + 10
            final int currentMonth = dateComponents.getMonth();
            this.supportedName = NEQUICK_BASE + String.format("ccir%02d.asc", currentMonth + 10);
            try (InputStream in = getStream(supportedName)) {
                loadData(in, supportedName);
            } catch (IOException e) {
                throw new OrekitException(OrekitMessages.INTERNAL_ERROR, e);
            }
            // Throw an exception if F2 or Fm3 were not loaded properly
            if (f2Loader == null || fm3Loader == null) {
                throw new OrekitException(OrekitMessages.NEQUICK_F2_FM3_NOT_LOADED, supportedName);
            }

        }

        /**
         * Load data from a stream.
         * @param input input stream
         * @param name name of the file
         * @throws IOException if data can't be read
         */
        public void loadData(final InputStream input, final String name)
            throws IOException {

            // Initialize arrays
            final double[][][] f2Temp  = new double[ROWS][TOTAL_COLUMNS_F2][DEPTH_F2];
            final double[][][] fm3Temp = new double[ROWS][TOTAL_COLUMNS_FM3][DEPTH_FM3];

            // Placeholders for parsed data
            int    lineNumber       = 0;
            int    index            = 0;
            int    currentRowF2     = 0;
            int    currentColumnF2  = 0;
            int    currentDepthF2   = 0;
            int    currentRowFm3    = 0;
            int    currentColumnFm3 = 0;
            int    currentDepthFm3  = 0;
            String line             = null;

            try (InputStreamReader isr = new InputStreamReader(input, StandardCharsets.UTF_8);
                 BufferedReader    br = new BufferedReader(isr)) {

                for (line = br.readLine(); line != null; line = br.readLine()) {
                    ++lineNumber;
                    line = line.trim();

                    // Read grid data
                    if (!line.isEmpty()) {
                        final String[] ccir_line = SEPARATOR.split(line);
                        for (final String field : ccir_line) {

                            if (index < NUMBER_F2_COEFFICIENTS) {
                                // Parse F2 coefficients
                                if (currentDepthF2 >= DEPTH_F2 && currentColumnF2 < (TOTAL_COLUMNS_F2 - 1)) {
                                    currentDepthF2 = 0;
                                    currentColumnF2++;
                                } else if (currentDepthF2 >= DEPTH_F2 && currentColumnF2 >= (TOTAL_COLUMNS_F2 - 1)) {
                                    currentDepthF2 = 0;
                                    currentColumnF2 = 0;
                                    currentRowF2++;
                                }
                                f2Temp[currentRowF2][currentColumnF2][currentDepthF2++] = Double.parseDouble(field);
                                index++;
                            } else {
                                // Parse Fm3 coefficients
                                if (currentDepthFm3 >= DEPTH_FM3 && currentColumnFm3 < (TOTAL_COLUMNS_FM3 - 1)) {
                                    currentDepthFm3 = 0;
                                    currentColumnFm3++;
                                } else if (currentDepthFm3 >= DEPTH_FM3 && currentColumnFm3 >= (TOTAL_COLUMNS_FM3 - 1)) {
                                    currentDepthFm3 = 0;
                                    currentColumnFm3 = 0;
                                    currentRowFm3++;
                                }
                                fm3Temp[currentRowFm3][currentColumnFm3][currentDepthFm3++] = Double.parseDouble(field);
                                index++;
                            }

                        }
                    }

                }

            } catch (NumberFormatException nfe) {
                throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
                                          lineNumber, name, line);
            }

            f2Loader  = f2Temp.clone();
            fm3Loader = fm3Temp.clone();

        }

    }

    /**
     * Container for ray-perigee parameters.
     * By convention, point 1 is at lower height.
     */
    private static class Ray {

        /** Threshold for ray-perigee parameters computation. */
        private static final double THRESHOLD = 1.0e-10;

        /** Distance of the first point from the ray perigee [m]. */
        private final double s1;

        /** Distance of the second point from the ray perigee [m]. */
        private final double s2;

        /** Ray-perigee radius [m]. */
        private final double rp;

        /** Ray-perigee latitude [rad]. */
        private final double latP;

        /** Ray-perigee longitude [rad]. */
        private final double lonP;

        /** Sine and cosine of ray-perigee latitude. */
        private final SinCos scLatP;

        /** Sine of azimuth of satellite as seen from ray-perigee. */
        private final double sinAzP;

        /** Cosine of azimuth of satellite as seen from ray-perigee. */
        private final double cosAzP;

        /**
         * Constructor.
         * @param recP receiver position
         * @param satP satellite position
         */
        Ray(final GeodeticPoint recP, final GeodeticPoint satP) {

            // Integration limits in meters (Eq. 140 and 141)
            final double r1 = RE + recP.getAltitude();
            final double r2 = RE + satP.getAltitude();

            // Useful parameters
            final double lat1     = recP.getLatitude();
            final double lat2     = satP.getLatitude();
            final double lon1     = recP.getLongitude();
            final double lon2     = satP.getLongitude();
            final SinCos scLatSat = FastMath.sinCos(lat2);
            final SinCos scLatRec = FastMath.sinCos(lat1);
            final SinCos scLon21  = FastMath.sinCos(lon2 - lon1);

            // Zenith angle computation (Eq. 153 to 155)
            // with added protection against numerical noise near zenith observation
            final double cosD = FastMath.min(1.0,
                                             scLatRec.sin() * scLatSat.sin() +
                                             scLatRec.cos() * scLatSat.cos() * scLon21.cos());
            final double sinD = FastMath.sqrt(1.0 - cosD * cosD);
            final double z    = FastMath.atan2(sinD, cosD - (r1 / r2));
            final SinCos scZ  = FastMath.sinCos(z);

            // Ray-perigee computation in meters (Eq. 156)
            this.rp = r1 * scZ.sin();

            // Ray-perigee latitude and longitude
            if (FastMath.abs(FastMath.abs(lat1) - 0.5 * FastMath.PI) < THRESHOLD) {
                // receiver is almost at North or South pole

                // Ray-perigee latitude (Eq. 157)
                this.latP = FastMath.copySign(z, lat1);

                // Ray-perigee longitude (Eq. 164)
                if (z < 0) {
                    this.lonP = lon2;
                } else {
                    this.lonP = lon2 + FastMath.PI;
                }

            } else if (FastMath.abs(scZ.sin()) < THRESHOLD) {
                // satellite is almost on receiver zenith

                this.latP = recP.getLatitude();
                this.lonP = recP.getLongitude();

            } else {

                // Ray-perigee latitude (Eq. 158 to 163)
                final double sinAz    = scLon21.sin() * scLatSat.cos() / sinD;
                final double cosAz    = (scLatSat.sin() - cosD * scLatRec.sin()) / (sinD * scLatRec.cos());
                final double sinLatP  = scLatRec.sin() * scZ.sin() - scLatRec.cos() * scZ.cos() * cosAz;
                final double cosLatP  = FastMath.sqrt(1.0 - sinLatP * sinLatP);
                this.latP = FastMath.atan2(sinLatP, cosLatP);

                // Ray-perigee longitude (Eq. 165 to 167)
                final double sinLonP = -sinAz * scZ.cos() / cosLatP;
                final double cosLonP = (scZ.sin() - scLatRec.sin() * sinLatP) / (scLatRec.cos() * cosLatP);
                this.lonP = FastMath.atan2(sinLonP, cosLonP) + lon1;

            }

            // Sine and cosine of ray-perigee latitude
            this.scLatP = FastMath.sinCos(latP);

            if (FastMath.abs(FastMath.abs(latP) - 0.5 * FastMath.PI) < THRESHOLD ||
                FastMath.abs(scZ.sin()) < THRESHOLD) {
                // Eq. 172 and 173
                this.sinAzP = 0.0;
                this.cosAzP = -FastMath.copySign(1, latP);
            } else {
                final SinCos scLon = FastMath.sinCos(lon2 - lonP);
                // Sine and cosine of azimuth of satellite as seen from ray-perigee
                final SinCos scPsi = FastMath.sinCos(greatCircleAngle(scLatSat, scLon));
                // Eq. 174 and 175
                this.sinAzP =  scLatSat.cos() * scLon.sin()                 /  scPsi.sin();
                this.cosAzP = (scLatSat.sin() - scLatP.sin() * scPsi.cos()) / (scLatP.cos() * scPsi.sin());
            }

            // Integration en points s1 and s2 in meters (Eq. 176 and 177)
            this.s1 = FastMath.sqrt(r1 * r1 - rp * rp);
            this.s2 = FastMath.sqrt(r2 * r2 - rp * rp);
        }

        /**
         * Get the distance of the first point from the ray perigee.
         * @return s1 in meters
         */
        public double getS1() {
            return s1;
        }

        /**
         * Get the distance of the second point from the ray perigee.
         * @return s2 in meters
         */
        public double getS2() {
            return s2;
        }

        /**
         * Get the ray-perigee radius.
         * @return the ray-perigee radius in meters
         */
        public double getRadius() {
            return rp;
        }

        /**
         * Get the ray-perigee latitude.
         * @return the ray-perigee latitude in radians
         */
        public double getLatitude() {
            return latP;
        }

        /**
         * Get the ray-perigee longitude.
         * @return the ray-perigee longitude in radians
         */
        public double getLongitude() {
            return lonP;
        }

        /**
         * Get the sine of azimuth of satellite as seen from ray-perigee.
         * @return the sine of azimuth
         */
        public double getSineAz() {
            return sinAzP;
        }

        /**
         * Get the cosine of azimuth of satellite as seen from ray-perigee.
         * @return the cosine of azimuth
         */
        public double getCosineAz() {
            return cosAzP;
        }

        /**
         * Compute the great circle angle from ray-perigee to satellite.
         * <p>
         * This method used the equations 168 to 171 pf the reference document.
         * </p>
         * @param scLat sine and cosine of satellite latitude
         * @param scLon sine and cosine of satellite longitude minus receiver longitude
         * @return the great circle angle in radians
         */
        private double greatCircleAngle(final SinCos scLat, final SinCos scLon) {
            if (FastMath.abs(FastMath.abs(latP) - 0.5 * FastMath.PI) < THRESHOLD) {
                return FastMath.abs(FastMath.asin(scLat.sin()) - latP);
            } else {
                final double cosPhi = scLatP.sin() * scLat.sin() +
                                      scLatP.cos() * scLat.cos() * scLon.cos();
                final double sinPhi = FastMath.sqrt(1.0 - cosPhi * cosPhi);
                return FastMath.atan2(sinPhi, cosPhi);
            }
        }
    }

    /**
     * Container for ray-perigee parameters.
     * By convention, point 1 is at lower height.
     */
    private static class FieldRay <T extends CalculusFieldElement<T>> {

        /** Threshold for ray-perigee parameters computation. */
        private static final double THRESHOLD = 1.0e-10;

        /** Distance of the first point from the ray perigee [m]. */
        private final T s1;

        /** Distance of the second point from the ray perigee [m]. */
        private final T s2;

        /** Ray-perigee radius [m]. */
        private final T rp;

        /** Ray-perigee latitude [rad]. */
        private final T latP;

        /** Ray-perigee longitude [rad]. */
        private final T lonP;

        /** Sine and cosine of ray-perigee latitude. */
        private final FieldSinCos<T> scLatP;

        /** Sine of azimuth of satellite as seen from ray-perigee. */
        private final T sinAzP;

        /** Cosine of azimuth of satellite as seen from ray-perigee. */
        private final T cosAzP;

        /**
         * Constructor.
         * @param field field of the elements
         * @param recP receiver position
         * @param satP satellite position
         */
        FieldRay(final Field<T> field, final FieldGeodeticPoint<T> recP, final FieldGeodeticPoint<T> satP) {

            // Integration limits in meters (Eq. 140 and 141)
            final T r1 = recP.getAltitude().add(RE);
            final T r2 = satP.getAltitude().add(RE);

            // Useful parameters
            final T pi   = r1.getPi();
            final T lat1 = recP.getLatitude();
            final T lat2 = satP.getLatitude();
            final T lon1 = recP.getLongitude();
            final T lon2 = satP.getLongitude();
            final FieldSinCos<T> scLatSat = FastMath.sinCos(lat2);
            final FieldSinCos<T> scLatRec = FastMath.sinCos(lat1);
            final FieldSinCos<T> scLon21  = FastMath.sinCos(lon2.subtract(lon1));

            // Zenith angle computation (Eq. 153 to 155)
            final T cosD = scLatRec.sin().multiply(scLatSat.sin()).
                            add(scLatRec.cos().multiply(scLatSat.cos()).multiply(scLon21.cos()));
            final T sinD = FastMath.sqrt(cosD.multiply(cosD).negate().add(1.0));
            final T z = FastMath.atan2(sinD, cosD.subtract(r1.divide(r2)));
            final FieldSinCos<T> scZ  = FastMath.sinCos(z);

            // Ray-perigee computation in meters (Eq. 156)
            this.rp = r1.multiply(scZ.sin());

            // Ray-perigee latitude and longitude
            if (FastMath.abs(FastMath.abs(lat1).subtract(pi.multiply(0.5)).getReal()) < THRESHOLD) {

                // Ray-perigee latitude (Eq. 157)
                this.latP = FastMath.copySign(z, lat1);

                // Ray-perigee longitude (Eq. 164)
                if (z.getReal() < 0) {
                    this.lonP = lon2;
                } else {
                    this.lonP = lon2.add(pi);
                }

            } else if (FastMath.abs(scZ.sin().getReal()) < THRESHOLD) {
                // satellite is almost on receiver zenith

                this.latP = recP.getLatitude();
                this.lonP = recP.getLongitude();

            } else {

                // Ray-perigee latitude (Eq. 158 to 163)
                final T sinAz    = FastMath.sin(lon2.subtract(lon1)).multiply(scLatSat.cos()).divide(sinD);
                final T cosAz    = scLatSat.sin().subtract(cosD.multiply(scLatRec.sin())).divide(sinD.multiply(scLatRec.cos()));
                final T sinLatP  = scLatRec.sin().multiply(scZ.sin()).subtract(scLatRec.cos().multiply(scZ.cos()).multiply(cosAz));
                final T cosLatP  = FastMath.sqrt(sinLatP.multiply(sinLatP).negate().add(1.0));
                this.latP = FastMath.atan2(sinLatP, cosLatP);

                // Ray-perigee longitude (Eq. 165 to 167)
                final T sinLonP = sinAz.negate().multiply(scZ.cos()).divide(cosLatP);
                final T cosLonP = scZ.sin().subtract(scLatRec.sin().multiply(sinLatP)).divide(scLatRec.cos().multiply(cosLatP));
                this.lonP = FastMath.atan2(sinLonP, cosLonP).add(lon1);

            }

            // Sine and cosine of ray-perigee latitude
            this.scLatP = FastMath.sinCos(latP);

            if (FastMath.abs(FastMath.abs(latP).subtract(pi.multiply(0.5)).getReal()) < THRESHOLD ||
                FastMath.abs(scZ.sin().getReal()) < THRESHOLD) {
                // Eq. 172 and 173
                this.sinAzP = field.getZero();
                this.cosAzP = FastMath.copySign(field.getOne(), latP).negate();
            } else {
                final FieldSinCos<T> scLon = FastMath.sinCos(lon2.subtract(lonP));
                // Sine and cosine of azimuth of satellite as seen from ray-perigee
                final FieldSinCos<T> scPsi = FastMath.sinCos(greatCircleAngle(scLatSat, scLon));
                // Eq. 174 and 175
                this.sinAzP = scLatSat.cos().multiply(scLon.sin()).divide(scPsi.sin());
                this.cosAzP = scLatSat.sin().subtract(scLatP.sin().multiply(scPsi.cos())).divide(scLatP.cos().multiply(scPsi.sin()));
            }

            // Integration en points s1 and s2 in meters (Eq. 176 and 177)
            this.s1 = FastMath.sqrt(r1.multiply(r1).subtract(rp.multiply(rp)));
            this.s2 = FastMath.sqrt(r2.multiply(r2).subtract(rp.multiply(rp)));
        }

        /**
         * Get the distance of the first point from the ray perigee.
         * @return s1 in meters
         */
        public T getS1() {
            return s1;
        }

        /**
         * Get the distance of the second point from the ray perigee.
         * @return s2 in meters
         */
        public T getS2() {
            return s2;
        }

        /**
         * Get the ray-perigee radius.
         * @return the ray-perigee radius in meters
         */
        public T getRadius() {
            return rp;
        }

        /**
         * Get the ray-perigee latitude.
         * @return the ray-perigee latitude in radians
         */
        public T getLatitude() {
            return latP;
        }

        /**
         * Get the ray-perigee longitude.
         * @return the ray-perigee longitude in radians
         */
        public T getLongitude() {
            return lonP;
        }

        /**
         * Get the sine of azimuth of satellite as seen from ray-perigee.
         * @return the sine of azimuth
         */
        public T getSineAz() {
            return sinAzP;
        }

        /**
         * Get the cosine of azimuth of satellite as seen from ray-perigee.
         * @return the cosine of azimuth
         */
        public T getCosineAz() {
            return cosAzP;
        }

        /**
         * Compute the great circle angle from ray-perigee to satellite.
         * <p>
         * This method used the equations 168 to 171 pf the reference document.
         * </p>
         * @param scLat sine and cosine of satellite latitude
         * @param scLon sine and cosine of satellite longitude minus receiver longitude
         * @return the great circle angle in radians
         */
        private T greatCircleAngle(final FieldSinCos<T> scLat, final FieldSinCos<T> scLon) {
            if (FastMath.abs(FastMath.abs(latP).getReal() - 0.5 * FastMath.PI) < THRESHOLD) {
                return FastMath.abs(FastMath.asin(scLat.sin()).subtract(latP));
            } else {
                final T cosPhi = scLatP.sin().multiply(scLat.sin()).add(
                                 scLatP.cos().multiply(scLat.cos()).multiply(scLon.cos()));
                final T sinPhi = FastMath.sqrt(cosPhi.multiply(cosPhi).negate().add(1.0));
                return FastMath.atan2(sinPhi, cosPhi);
            }
        }
    }

    /** Performs the computation of the coordinates along the integration path. */
    private static class Segment {

        /** Threshold for zenith segment. */
        private static final double THRESHOLD = 1.0e-3;

        /** Latitudes [rad]. */
        private final double[] latitudes;

        /** Longitudes [rad]. */
        private final double[] longitudes;

        /** Heights [m]. */
        private final double[] heights;

        /** Integration step [m]. */
        private final double deltaN;

        /**
         * Constructor.
         * @param n number of points used for the integration
         * @param ray ray-perigee parameters
         */
        Segment(final int n, final Ray ray) {
            // Integration en points
            final double s1 = ray.getS1();
            final double s2 = ray.getS2();

            // Integration step (Eq. 195)
            this.deltaN = (s2 - s1) / n;

            // Segments
            final double[] s = getSegments(n, s1);

            // Useful parameters
            final double rp = ray.getRadius();
            final SinCos scLatP = FastMath.sinCos(ray.getLatitude());

            // Geodetic coordinates
            final int size = s.length;
            heights    = new double[size];
            latitudes  = new double[size];
            longitudes = new double[size];
            for (int i = 0; i < size; i++) {
                // Heights (Eq. 178)
                heights[i] = FastMath.sqrt(s[i] * s[i] + rp * rp) - RE;

                if (rp < THRESHOLD) {
                    // zenith segment
                    latitudes[i]  = ray.getLatitude();
                    longitudes[i] = ray.getLongitude();
                } else {
                    // Great circle parameters (Eq. 179 to 181)
                    final double tanDs = s[i] / rp;
                    final double cosDs = 1.0 / FastMath.sqrt(1.0 + tanDs * tanDs);
                    final double sinDs = tanDs * cosDs;

                    // Latitude (Eq. 182 to 183)
                    final double sinLatS = scLatP.sin() * cosDs + scLatP.cos() * sinDs * ray.getCosineAz();
                    final double cosLatS = FastMath.sqrt(1.0 - sinLatS * sinLatS);
                    latitudes[i] = FastMath.atan2(sinLatS, cosLatS);

                    // Longitude (Eq. 184 to 187)
                    final double sinLonS = sinDs * ray.getSineAz() * scLatP.cos();
                    final double cosLonS = cosDs - scLatP.sin() * sinLatS;
                    longitudes[i] = FastMath.atan2(sinLonS, cosLonS) + ray.getLongitude();
                }
            }
        }

        /**
         * Computes the distance of a point from the ray-perigee.
         * @param n number of points used for the integration
         * @param s1 lower boundary
         * @return the distance of a point from the ray-perigee in km
         */
        private double[] getSegments(final int n, final double s1) {
            // Eq. 196
            final double g = 0.5773502691896 * deltaN;
            // Eq. 197
            final double y = s1 + (deltaN - g) * 0.5;
            final double[] segments = new double[2 * n];
            int index = 0;
            for (int i = 0; i < n; i++) {
                // Eq. 198
                segments[index] = y + i * deltaN;
                index++;
                segments[index] = y + i * deltaN + g;
                index++;
            }
            return segments;
        }

        /**
         * Get the latitudes of the coordinates along the integration path.
         * @return the latitudes in radians
         */
        public double[] getLatitudes() {
            return latitudes;
        }

        /**
         * Get the longitudes of the coordinates along the integration path.
         * @return the longitudes in radians
         */
        public double[] getLongitudes() {
            return longitudes;
        }

        /**
         * Get the heights of the coordinates along the integration path.
         * @return the heights in m
         */
        public double[] getHeights() {
            return heights;
        }

        /**
         * Get the integration step.
         * @return the integration step in meters
         */
        public double getInterval() {
            return deltaN;
        }
    }

    /** Performs the computation of the coordinates along the integration path. */
    private static class FieldSegment <T extends CalculusFieldElement<T>> {

        /** Threshold for zenith segment. */
        private static final double THRESHOLD = 1.0e-3;

        /** Latitudes [rad]. */
        private final T[] latitudes;

        /** Longitudes [rad]. */
        private final T[] longitudes;

        /** Heights [m]. */
        private final T[] heights;

        /** Integration step [m]. */
        private final T deltaN;

        /**
         * Constructor.
         * @param field field of the elements
         * @param n number of points used for the integration
         * @param ray ray-perigee parameters
         */
        FieldSegment(final Field<T> field, final int n, final FieldRay<T> ray) {
            // Integration en points
            final T s1 = ray.getS1();
            final T s2 = ray.getS2();

            // Integration step (Eq. 195)
            this.deltaN = s2.subtract(s1).divide(n);

            // Segments
            final T[] s = getSegments(field, n, s1);

            // Useful parameters
            final T rp = ray.getRadius();
            final FieldSinCos<T> scLatP = FastMath.sinCos(ray.getLatitude());

            // Geodetic coordinates
            final int size = s.length;
            heights    = MathArrays.buildArray(field, size);
            latitudes  = MathArrays.buildArray(field, size);
            longitudes = MathArrays.buildArray(field, size);
            for (int i = 0; i < size; i++) {
                // Heights (Eq. 178)
                heights[i] = FastMath.sqrt(s[i].multiply(s[i]).add(rp.multiply(rp))).subtract(RE);

                if (rp.getReal() < THRESHOLD) {
                    // zenith segment
                    latitudes[i]  = ray.getLatitude();
                    longitudes[i] = ray.getLongitude();
                } else {
                    // Great circle parameters (Eq. 179 to 181)
                    final T tanDs = s[i].divide(rp);
                    final T cosDs = FastMath.sqrt(tanDs.multiply(tanDs).add(1.0)).reciprocal();
                    final T sinDs = tanDs.multiply(cosDs);

                    // Latitude (Eq. 182 to 183)
                    final T sinLatS = scLatP.sin().multiply(cosDs).add(scLatP.cos().multiply(sinDs).multiply(ray.getCosineAz()));
                    final T cosLatS = FastMath.sqrt(sinLatS.multiply(sinLatS).negate().add(1.0));
                    latitudes[i] = FastMath.atan2(sinLatS, cosLatS);

                    // Longitude (Eq. 184 to 187)
                    final T sinLonS = sinDs.multiply(ray.getSineAz()).multiply(scLatP.cos());
                    final T cosLonS = cosDs.subtract(scLatP.sin().multiply(sinLatS));
                    longitudes[i] = FastMath.atan2(sinLonS, cosLonS).add(ray.getLongitude());
                }
            }
        }

        /**
         * Computes the distance of a point from the ray-perigee.
         * @param field field of the elements
         * @param n number of points used for the integration
         * @param s1 lower boundary
         * @return the distance of a point from the ray-perigee in km
         */
        private T[] getSegments(final Field<T> field, final int n, final T s1) {
            // Eq. 196
            final T g = deltaN.multiply(0.5773502691896);
            // Eq. 197
            final T y = s1.add(deltaN.subtract(g).multiply(0.5));
            final T[] segments = MathArrays.buildArray(field, 2 * n);
            int index = 0;
            for (int i = 0; i < n; i++) {
                // Eq. 198
                segments[index] = y.add(deltaN.multiply(i));
                index++;
                segments[index] = y.add(deltaN.multiply(i)).add(g);
                index++;
            }
            return segments;
        }

        /**
         * Get the latitudes of the coordinates along the integration path.
         * @return the latitudes in radians
         */
        public T[] getLatitudes() {
            return latitudes;
        }

        /**
         * Get the longitudes of the coordinates along the integration path.
         * @return the longitudes in radians
         */
        public T[] getLongitudes() {
            return longitudes;
        }

        /**
         * Get the heights of the coordinates along the integration path.
         * @return the heights in m
         */
        public T[] getHeights() {
            return heights;
        }

        /**
         * Get the integration step.
         * @return the integration step in meters
         */
        public T getInterval() {
            return deltaN;
        }
    }

}