/* Copyright 2002-2025 Thales Alenia Space
    * Licensed to CS Communication & Systèmes (CS) under one or more
    * contributor license agreements.  See the NOTICE file distributed with
    * this work for additional information regarding copyright ownership.
    * CS licenses this file to You under the Apache License, Version 2.0
    * (the "License"); you may not use this file except in compliance with
    * the License.  You may obtain a copy of the License at
    *
    *   http://www.apache.org/licenses/LICENSE-2.0
   *
   * Unless required by applicable law or agreed to in writing, software
   * distributed under the License is distributed on an "AS IS" BASIS,
   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   * See the License for the specific language governing permissions and
   * limitations under the License.
   */
  package org.orekit.models.earth.ionosphere.nequick;
  
  import org.hipparchus.CalculusFieldElement;
  import org.hipparchus.util.FastMath;
  import org.hipparchus.util.MathArrays;
  import org.orekit.annotation.DefaultDataContext;
  import org.orekit.bodies.FieldGeodeticPoint;
  import org.orekit.bodies.GeodeticPoint;
  import org.orekit.data.DataContext;
  import org.orekit.data.DataSource;
  import org.orekit.errors.OrekitException;
  import org.orekit.errors.OrekitMessages;
  import org.orekit.time.DateTimeComponents;
  import org.orekit.time.TimeScale;
  
  /** Galileo-specific version of NeQuick engine.
   * @since 13.0
   * @author Luc Maisonobe
   */
  public class NeQuickGalileo extends NeQuickModel {
  
      /** Starting number of points for integration. */
      private static final int N_START = 8;
  
      /** Broadcast ionization engine coefficients. */
      private final double[] alpha;
  
      /**
       * Build a new instance.
       *
       * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
       *
       * @param alpha effective ionisation level coefficients
       * @see #NeQuickGalileo(double[], TimeScale)
       */
      @DefaultDataContext
      public NeQuickGalileo(final double[] alpha) {
          this(alpha, DataContext.getDefault().getTimeScales().getUTC());
      }
  
      /**
       * Build a new instance of the Galileo version of the NeQuick-2 model.
       * <p>
       * The Galileo version uses a loose modip grid and 3 broadcast parameters to compute
       * effective ionization level.
       * </p>
       * @param alpha broadcast effective ionisation level coefficients
       * @param utc UTC time scale.
       * @since 10.1
       */
      public NeQuickGalileo(final double[] alpha, final TimeScale utc) {
          super(utc);
          this.alpha = alpha.clone();
      }
  
      /** Get effective ionisation level coefficients.
       * @return effective ionisation level coefficients
       */
      public double[] getAlpha() {
          return alpha.clone();
      }
  
      /**
       * Compute effective ionization level.
       *
       * @param modip modified dip latitude at receiver location
       * @return effective ionization level (Az in Nequick Galileo, R12 in original Nequick ITU)
       */
      private double effectiveIonizationLevel(final double modip) {
          // Particular condition (Eq. 17)
          if (alpha[0] == 0.0 && alpha[1] == 0.0 && alpha[2] == 0.0) {
              return 63.7;
          } else {
              // Az = a0 + modip * a1 + modip² * a2 (Eq. 18)
              return FastMath.min(FastMath.max(alpha[0] + modip * (alpha[1] + modip * alpha[2]), 0.0), 400.0);
          }
      }
  
      /**
       * Compute effective ionization level.
       *
       * @param <T>   type of the field elements
       * @param modip modified dip latitude at receiver location
      * @return effective ionization level (Az in Nequick Galileo, R12 in original Nequick ITU)
      */
     private <T extends CalculusFieldElement<T>> T effectiveIonizationLevel(final T modip) {
         // Particular condition (Eq. 17)
         if (alpha[0] == 0.0 && alpha[1] == 0.0 && alpha[2] == 0.0) {
             return modip.newInstance(63.7);
         } else {
             // Az = a0 + modip * a1 + modip² * a2 (Eq. 18)
             return FastMath.min(FastMath.max(modip.multiply(alpha[2]).add(alpha[1]).multiply(modip).add(alpha[0]),
                                              0.0),
                                 400.0);
         }
     }
 
     /** {@inheritDoc} */
     @Override
     double stec(final DateTimeComponents dateTime, final Ray ray) {
 
         // Tolerance for the integration accuracy. Defined inside the reference document, section 2.5.8.1.
         final double h1 = ray.getRecH();
         final double tolerance;
         if (h1 < 1000000.0) {
             tolerance = 0.001;
         } else {
             tolerance = 0.01;
         }
 
         // Integration
         int n = N_START;
         final Segment seg1 = new Segment(n, ray, ray.getS1(), ray.getS2());
         double gn1 = stecIntegration(dateTime, seg1);
         n *= 2;
         final Segment seg2 = new Segment(n, ray, ray.getS1(), ray.getS2());
         double gn2 = stecIntegration(dateTime, seg2);
 
         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, ray.getS1(), ray.getS2());
             gn2 = stecIntegration(dateTime, seg);
             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;
 
     }
 
     /** {@inheritDoc} */
     @Override
     <T extends CalculusFieldElement<T>> T stec(final DateTimeComponents dateTime, final FieldRay<T> ray) {
 
         // Tolerance for the integration accuracy. Defined inside the reference document, section 2.5.8.1.
         final T h1 = ray.getRecH();
         final double tolerance;
         if (h1.getReal() < 1000000.0) {
             tolerance = 0.001;
         } else {
             tolerance = 0.01;
         }
 
         // Integration
         int n = N_START;
         final FieldSegment<T> seg1 = new FieldSegment<>(n, ray, ray.getS1(), ray.getS2());
         T gn1 = stecIntegration(dateTime, seg1);
         n *= 2;
         final FieldSegment<T> seg2 = new FieldSegment<>(n, ray, ray.getS1(), ray.getS2());
         T gn2 = stecIntegration(dateTime, seg2);
 
         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<>(n, ray, ray.getS1(), ray.getS2());
             gn2 = stecIntegration(dateTime, seg);
             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 performs the STEC integration.
      *
      * @param dateTime current date and time components
      * @param seg      coordinates along the integration path
      * @return result of the integration
      */
     private double stecIntegration(final DateTimeComponents dateTime, final Segment seg) {
 
         // Compute electron density
         double density = 0.0;
         for (int i = 0; i < seg.getNbPoints(); i++) {
             final GeodeticPoint gp = seg.getPoint(i);
             final double modip = GalileoHolder.INSTANCE.computeMODIP(gp.getLatitude(), gp.getLongitude());
             final double az = effectiveIonizationLevel(modip);
             density += electronDensity(dateTime, az, gp.getLatitude(), gp.getLongitude(), gp.getAltitude());
         }
 
         return 0.5 * seg.getInterval() * density;
     }
 
     /**
      * This method performs the STEC integration.
      *
      * @param <T>      type of the elements
      * @param dateTime current date and time components
      * @param seg      coordinates along the integration path
      * @return result of the integration
      */
     private <T extends CalculusFieldElement<T>> T stecIntegration(final DateTimeComponents dateTime,
                                                               final FieldSegment<T> seg) {
         // Compute electron density
         T density = seg.getInterval().getField().getZero();
         for (int i = 0; i < seg.getNbPoints(); i++) {
             final FieldGeodeticPoint<T> gp = seg.getPoint(i);
             final T modip = GalileoHolder.INSTANCE.computeMODIP(gp.getLatitude(), gp.getLongitude());
             final T az = effectiveIonizationLevel(modip);
             density = density.add(electronDensity(dateTime, az,
                                                   gp.getLatitude(), gp.getLongitude(), gp.getAltitude()));
         }
 
         return seg.getInterval().multiply(density).multiply(0.5);
     }
 
     /** {@inheritDoc} */
     @Override
     protected double computeMODIP(final double latitude, final double longitude) {
         return GalileoHolder.INSTANCE.computeMODIP(latitude, longitude);
     }
 
     /** {@inheritDoc} */
     @Override
     protected <T extends CalculusFieldElement<T>> T computeMODIP(final T latitude, final T longitude) {
         return GalileoHolder.INSTANCE.computeMODIP(latitude, longitude);
     }
 
     /** {@inheritDoc} */
     @Override
     boolean applyChapmanParameters(final double hInKm) {
         return hInKm < 100.0;
     }
 
     /** {@inheritDoc} */
     @Override
     double[] computeExponentialArguments(final double h, final NeQuickParameters parameters) {
 
         // If h < 100km we use h = 100km as recommended in the reference document
         // for equations 111 to 113
         final double clippedH = FastMath.max(h, 100.0);
 
         // Eq. 111 to 113
         final double   exp       = clipExp(10.0 / (1.0 + FastMath.abs(clippedH - parameters.getHmF2())));
         final double[] arguments = new double[3];
         arguments[0] =  (clippedH - parameters.getHmF2()) / parameters.getB2Bot();
         arguments[1] = ((clippedH - parameters.getHmF1()) / parameters.getBF1(h)) * exp;
         arguments[2] = ((clippedH - parameters.getHmE())  / parameters.getBE(h))  * exp;
         return arguments;
 
     }
 
     /** {@inheritDoc} */
     @Override
     <T extends CalculusFieldElement<T>> T[] computeExponentialArguments(final T h,
                                                                         final FieldNeQuickParameters<T> parameters) {
         // If h < 100km we use h = 100km as recommended in the reference document
         // for equations 111 to 113
         final T clippedH = FastMath.max(h, 100);
 
         // Eq. 111 to 113
         final T   exp   = clipExp(FastMath.abs(clippedH.subtract(parameters.getHmF2())).add(1.0).reciprocal().multiply(10.0));
         final T[] arguments = MathArrays.buildArray(h.getField(), 3);
         arguments[0] = clippedH.subtract(parameters.getHmF2()).divide(parameters.getB2Bot());
         arguments[1] = clippedH.subtract(parameters.getHmF1()).divide(parameters.getBF1(h)).multiply(exp);
         arguments[2] = clippedH.subtract(parameters.getHmE()).divide(parameters.getBE(h)).multiply(exp);
         return arguments;
 
     }
 
     /** {@inheritDoc} */
     @Override
     double computeH0(final NeQuickParameters parameters) {
 
         final int month = parameters.getDateTime().getDate().getMonth();
 
         // Auxiliary parameter ka (Eq. 99 and 100)
         final double ka;
         if (month > 3 && month < 10) {
             // month = 4,5,6,7,8,9
             ka = 6.705 - 0.014 * parameters.getAzr() - 0.008 * parameters.getHmF2();
         } else {
             // month = 1,2,3,10,11,12
             final double ratio = parameters.getHmF2() / parameters.getB2Bot();
             ka = -7.77 + 0.097 * ratio * ratio + 0.153 * parameters.getNmF2();
         }
 
         // Auxiliary parameter kb (Eq. 101 and 102)
         double kb = parameters.join(ka, 2.0, 1.0, ka - 2.0);
         kb = parameters.join(8.0, kb, 1.0, kb - 8.0);
 
         // Auxiliary parameter Ha (Eq. 103)
         final double hA = kb * parameters.getB2Bot();
 
         // Auxiliary parameters x and v (Eq. 104 and 105)
         final double x = 0.01 * (hA - 150.0);
         final double v = (0.041163 * x - 0.183981) * x + 1.424472;
 
         // Topside thickness parameter (Eq. 106)
         return hA / v;
 
     }
 
     /** {@inheritDoc} */
     @Override
     <T extends CalculusFieldElement<T>> T computeH0(final FieldNeQuickParameters<T> parameters) {
 
         final int month = parameters.getDateTime().getDate().getMonth();
 
         // One
         final T one = parameters.getAzr().getField().getOne();
 
         // Auxiliary parameter ka (Eq. 99 and 100)
         final T ka;
         if (month > 3 && month < 10) {
             // month = 4,5,6,7,8,9
             ka = parameters.getAzr().multiply(0.014).add(parameters.getHmF2().multiply(0.008)).negate().add(6.705);
         } else {
             // month = 1,2,3,10,11,12
             final T ratio = parameters.getHmF2().divide(parameters.getB2Bot());
             ka = ratio.multiply(ratio).multiply(0.097).add(parameters.getNmF2().multiply(0.153)).add(-7.77);
         }
 
         // Auxiliary parameter kb (Eq. 101 and 102)
         T kb = parameters.join(ka, one.newInstance(2.0), one, ka.subtract(2.0));
         kb = parameters.join(one.newInstance(8.0), kb, one, kb.subtract(8.0));
 
         // Auxiliary parameter Ha (Eq. 103)
         final T hA = kb.multiply(parameters.getB2Bot());
 
         // Auxiliary parameters x and v (Eq. 104 and 105)
         final T x = hA.subtract(150.0).multiply(0.01);
         final T v = x.multiply(0.041163).subtract(0.183981).multiply(x).add(1.424472);
 
         // Topside thickness parameter (Eq. 106)
         return hA.divide(v);
 
     }
 
     /** Holder for the Galileo-specific modip singleton.
      * <p>
      * We use the initialization on demand holder idiom to store the singleton,
      * as it is both thread-safe, efficient (no synchronization) and works with
      * all versions of java.
      * </p>
      */
     private static class GalileoHolder {
 
         /** Resource for modip grid. */
         private static final String MODIP_GRID = "/assets/org/orekit/nequick/modipNeQG_wrapped.asc";
 
         /** Unique instance. */
         private static final ModipGrid INSTANCE =
             new ModipGrid(36, 36,
                           new DataSource(MODIP_GRID, () -> GalileoHolder.class.getResourceAsStream(MODIP_GRID)),
                           true);
 
         /** Private constructor.
          * <p>This class is a utility class, it should neither have a public
          * nor a default constructor. This private constructor prevents
          * the compiler from generating one automatically.</p>
          */
         private GalileoHolder() {
         }
 
     }
 
 }