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