/* Copyright 2002-2026 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.forces.radiation;
  
  import java.util.List;
  
  import org.hipparchus.CalculusFieldElement;
  import org.hipparchus.analysis.polynomials.PolynomialFunction;
  import org.hipparchus.analysis.polynomials.PolynomialsUtils;
  import org.hipparchus.geometry.euclidean.threed.FieldRotation;
  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  import org.hipparchus.geometry.euclidean.threed.Rotation;
  import org.hipparchus.geometry.euclidean.threed.RotationConvention;
  import org.hipparchus.geometry.euclidean.threed.Vector3D;
  import org.hipparchus.util.FastMath;
  import org.hipparchus.util.FieldSinCos;
  import org.hipparchus.util.MathUtils;
  import org.hipparchus.util.SinCos;
  import org.orekit.annotation.DefaultDataContext;
  import org.orekit.data.DataContext;
  import org.orekit.forces.ForceModel;
  import org.orekit.frames.Frame;
  import org.orekit.propagation.FieldSpacecraftState;
  import org.orekit.propagation.SpacecraftState;
  import org.orekit.time.AbsoluteDate;
  import org.orekit.time.FieldAbsoluteDate;
  import org.orekit.time.TimeScale;
  import org.orekit.utils.Constants;
  import org.orekit.utils.ExtendedPositionProvider;
  import org.orekit.utils.drivers.ParameterDriver;
  
  /** The Knocke Earth Albedo and IR emission force model.
   * <p>
   * This model is based on "EARTH RADIATION PRESSURE EFFECTS ON SATELLITES", 1988, by P. C. Knocke, J. C. Ries, and B. D. Tapley.
   * </p> <p>
   * This model represents the effects of radiation pressure coming from the Earth.
   * It considers Solar radiation which has been reflected by Earth (albedo) and Earth infrared emissions.
   * The planet is considered as a sphere and is divided into elementary areas.
   * Each elementary area is considered as a plane and emits radiation according to Lambert's law.
   * The flux the satellite receives is then equal to the sum of the elementary fluxes coming from Earth.
   * </p> <p>
   * The radiative model of the satellite, and its ability to diffuse, reflect  or absorb radiation is handled
   * by a {@link RadiationSensitive radiation sensitive model}.
   * </p> <p>
   * <b>Caution:</b> This model is only suitable for Earth. Using it with another central body is prone to error..
   * </p>
   *
   * @author Thomas Paulet
   * @since 10.3
   */
  public class KnockeRediffusedForceModel implements ForceModel {
  
      /** Earth rotation around Sun pulsation in rad/sec. */
      public static final double EARTH_AROUND_SUN_PULSATION = MathUtils.TWO_PI / Constants.JULIAN_YEAR;
  
      /** Coefficient for solar irradiance computation. */
      public static final double ES_COEFF = 4.5606E-6;
  
      /** First coefficient for albedo computation. */
      public static final double A0 = 0.34;
  
      /** Second coefficient for albedo computation. */
      public static final double C0 = 0.;
  
      /** Third coefficient for albedo computation. */
      public static final double C1 = 0.10;
  
      /** Fourth coefficient for albedo computation. */
      public static final double C2 = 0.;
  
      /** Fifth coefficient for albedo computation. */
      public static final double A2 = 0.29;
  
      /** First coefficient for Earth emissivity computation. */
      public static final double E0 = 0.68;
  
      /** Second coefficient for Earth emissivity computation. */
      public static final double K0 = 0.;
  
      /** Third coefficient for Earth emissivity computation. */
      public static final double K1 = -0.07;
  
      /** Fourth coefficient for Earth emissivity computation. */
      public static final double K2 = 0.;
  
     /** Fifth coefficient for Earth emissivity computation. */
     public static final double E2 = -0.18;
 
     /** Sun model. */
     private final ExtendedPositionProvider sun;
 
     /** Spacecraft. */
     private final RadiationSensitive spacecraft;
 
     /** Angular resolution for emissivity and albedo computation in rad. */
     private final double angularResolution;
 
     /** Earth equatorial radius in m. */
     private final double equatorialRadius;
 
     /** Reference date for periodic terms: December 22nd 1981.
      * Without more precision, the choice is to set it at midnight, UTC. */
     private final AbsoluteDate referenceEpoch;
 
     /** Default Constructor.
      * <p>This constructor uses the {@link DataContext#getDefault() default data context}</p>.
      * @param sun Sun model
      * @param spacecraft the object physical and geometrical information
      * @param equatorialRadius the Earth equatorial radius in m
      * @param angularResolution angular resolution in rad
      */
     @DefaultDataContext
     public KnockeRediffusedForceModel (final ExtendedPositionProvider sun,
                                        final RadiationSensitive spacecraft,
                                        final double equatorialRadius,
                                        final double angularResolution) {
 
         this(sun, spacecraft, equatorialRadius, angularResolution, DataContext.getDefault().getTimeScales().getUTC());
     }
 
     /** General constructor.
      * @param sun Sun model
      * @param spacecraft the object physical and geometrical information
      * @param equatorialRadius the Earth equatorial radius in m
      * @param angularResolution angular resolution in rad
      * @param utc the UTC time scale to define reference epoch
      */
     public KnockeRediffusedForceModel (final ExtendedPositionProvider sun,
                                        final RadiationSensitive spacecraft,
                                        final double equatorialRadius,
                                        final double angularResolution,
                                        final TimeScale utc) {
         this.sun               = sun;
         this.spacecraft        = spacecraft;
         this.equatorialRadius  = equatorialRadius;
         this.angularResolution = angularResolution;
         this.referenceEpoch    = new AbsoluteDate(1981, 12, 22, 0, 0, 0.0, utc);
     }
 
     /** {@inheritDoc} */
     @Override
     public Vector3D acceleration(final SpacecraftState s,
                                  final double[] parameters) {
 
         // Get date
         final AbsoluteDate date = s.getDate();
 
         // Get frame
         final Frame frame = s.getFrame();
 
         // Get satellite position
         final Vector3D satellitePosition = s.getPosition();
 
         // Get Sun position
         final Vector3D sunPosition = sun.getPosition(date, frame);
 
         // Project satellite on Earth as vector
         final Vector3D projectedToGround = satellitePosition.normalize().scalarMultiply(equatorialRadius);
 
         // Get elementary vector east for Earth browsing using rotations
         final double q = 1.0 / FastMath.hypot(satellitePosition.getX(), satellitePosition.getY());
         final Vector3D east = new Vector3D(-q * satellitePosition.getY(), q * satellitePosition.getX(), 0);
 
         // Initialize rediffused flux with elementary flux coming from the circular area around the projected satellite
         final double centerArea = MathUtils.TWO_PI * equatorialRadius * equatorialRadius *
                                  (1.0 - FastMath.cos(angularResolution));
         Vector3D rediffusedFlux = computeElementaryFlux(s, projectedToGround, sunPosition, centerArea);
 
         // Sectorize the part of Earth which is seen by the satellite into crown sectors with constant angular resolution
         for (double eastAxisOffset = 1.5 * angularResolution;
              eastAxisOffset < FastMath.acos(equatorialRadius / satellitePosition.getNorm());
              eastAxisOffset = eastAxisOffset + angularResolution) {
 
             // Build rotation transformations to get first crown elementary sector center
             final Rotation eastRotation = new Rotation(east, eastAxisOffset, RotationConvention.VECTOR_OPERATOR);
 
             // Get first elementary crown sector center
             final Vector3D firstCrownSectorCenter = eastRotation.applyTo(projectedToGround);
 
             // Compute current elementary crown sector area, it results of the integration of an elementary crown sector
             // over the angular resolution
             final double sectorArea = equatorialRadius * equatorialRadius *
                                       2.0 * angularResolution * FastMath.sin(0.5 * angularResolution) *
                                       FastMath.sin(eastAxisOffset);
 
             // Browse the entire crown
             for (double radialAxisOffset = 0.5 * angularResolution;
                  radialAxisOffset < MathUtils.TWO_PI;
                  radialAxisOffset = radialAxisOffset + angularResolution) {
 
                 // Build rotation transformations to get elementary area center
                 final Rotation radialRotation  = new Rotation(projectedToGround, radialAxisOffset, RotationConvention.VECTOR_OPERATOR);
 
                 // Get current elementary crown sector center
                 final Vector3D currentCenter = radialRotation.applyTo(firstCrownSectorCenter);
 
                 // Add current sector contribution to total rediffused flux
                 rediffusedFlux = rediffusedFlux.add(computeElementaryFlux(s, currentCenter, sunPosition, sectorArea));
             }
         }
 
         return spacecraft.radiationPressureAcceleration(s, rediffusedFlux, parameters);
     }
 
 
     /** {@inheritDoc} */
     @Override
     public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(final FieldSpacecraftState<T> s,
                                                                              final T[] parameters) {
         // Get date
         final FieldAbsoluteDate<T> date = s.getDate();
 
         // Get frame
         final Frame frame = s.getFrame();
 
         // Get zero
         final T zero = date.getField().getZero();
 
         // Get satellite position
         final FieldVector3D<T> satellitePosition = s.getPosition();
 
         // Get Sun position
         final FieldVector3D<T> sunPosition = sun.getPosition(date, frame);
 
         // Project satellite on Earth as vector
         final FieldVector3D<T> projectedToGround = satellitePosition.normalize().scalarMultiply(equatorialRadius);
 
         // Get elementary vector east for Earth browsing using rotations
         final T q = FastMath.hypot(satellitePosition.getX(), satellitePosition.getY()).reciprocal();
         final FieldVector3D<T> east = new FieldVector3D<>(q.negate().multiply(satellitePosition.getY()),
                                                           q.multiply(satellitePosition.getX()),
                                                           zero);
 
         // Initialize rediffused flux with elementary flux coming from the circular area around the projected satellite
         final T centerArea = zero.getPi().multiply(2.0).multiply(equatorialRadius).multiply(equatorialRadius).
                         multiply(1.0 - FastMath.cos(angularResolution));
         FieldVector3D<T> rediffusedFlux = computeElementaryFlux(s, projectedToGround, sunPosition, centerArea);
 
         // Sectorize the part of Earth which is seen by the satellite into crown sectors with constant angular resolution
         for (double eastAxisOffset = 1.5 * angularResolution;
              eastAxisOffset < FastMath.acos(equatorialRadius / satellitePosition.getNorm().getReal());
              eastAxisOffset = eastAxisOffset + angularResolution) {
 
             // Build rotation transformations to get first crown elementary sector center
             final FieldRotation<T> eastRotation = new FieldRotation<>(east, zero.newInstance(eastAxisOffset),
                                                                       RotationConvention.VECTOR_OPERATOR);
 
             // Get first elementary crown sector center
             final FieldVector3D<T> firstCrownSectorCenter = eastRotation.applyTo(projectedToGround);
 
             // Compute current elementary crown sector area, it results of the integration of an elementary crown sector
             // over the angular resolution
             final T sectorArea = zero.newInstance(equatorialRadius * equatorialRadius *
                                                   2.0 * angularResolution * FastMath.sin(0.5 * angularResolution) *
                                                   FastMath.sin(eastAxisOffset));
 
             // Browse the entire crown
             for (double radialAxisOffset = 0.5 * angularResolution;
                  radialAxisOffset < MathUtils.TWO_PI;
                  radialAxisOffset = radialAxisOffset + angularResolution) {
 
                 // Build rotation transformations to get elementary area center
                 final FieldRotation<T> radialRotation  = new FieldRotation<>(projectedToGround,
                                                                              zero.newInstance(radialAxisOffset),
                                                                              RotationConvention.VECTOR_OPERATOR);
 
                 // Get current elementary crown sector center
                 final FieldVector3D<T> currentCenter = radialRotation.applyTo(firstCrownSectorCenter);
 
                 // Add current sector contribution to total rediffused flux
                 rediffusedFlux = rediffusedFlux.add(computeElementaryFlux(s, currentCenter, sunPosition, sectorArea));
             }
         }
 
         return spacecraft.radiationPressureAcceleration(s, rediffusedFlux, parameters);
     }
 
 
     /** {@inheritDoc} */
     @Override
     public List<ParameterDriver> getParametersDrivers() {
         return spacecraft.getRadiationParametersDrivers();
     }
 
     /** Compute Earth albedo.
      * Albedo value represents the fraction of solar radiative flux that is reflected by Earth.
      * Its value is in [0;1].
      * @param date the date
      * @param phi the equatorial latitude in rad
      * @return the albedo in [0;1]
      */
     public double computeAlbedo(final AbsoluteDate date, final double phi) {
 
         // Get duration since coefficient reference epoch
         final double deltaT = date.durationFrom(referenceEpoch);
 
         // Compute 1rst Legendre polynomial coeficient
         final SinCos sc = FastMath.sinCos(EARTH_AROUND_SUN_PULSATION * deltaT);
         final double A1 = C0 +
                           C1 * sc.cos() +
                           C2 * sc.sin();
 
         // Get 1rst and 2nd order Legendre polynomials
         final PolynomialFunction firstLegendrePolynomial  = PolynomialsUtils.createLegendrePolynomial(1);
         final PolynomialFunction secondLegendrePolynomial = PolynomialsUtils.createLegendrePolynomial(2);
 
         // Get latitude sinus
         final double sinPhi = FastMath.sin(phi);
 
         // Compute albedo
         return A0 +
                A1 * firstLegendrePolynomial.value(sinPhi) +
                A2 * secondLegendrePolynomial.value(sinPhi);
 
     }
 
 
     /** Compute Earth albedo.
      * Albedo value represents the fraction of solar radiative flux that is reflected by Earth.
      * Its value is in [0;1].
      * @param date the date
      * @param phi the equatorial latitude in rad
      * @param <T> extends CalculusFieldElement
      * @return the albedo in [0;1]
      */
     public <T extends CalculusFieldElement<T>> T computeAlbedo(final FieldAbsoluteDate<T> date, final T phi) {
 
         // Get duration since coefficient reference epoch
         final T deltaT = date.durationFrom(referenceEpoch);
 
         // Compute 1rst Legendre polynomial coeficient
         final FieldSinCos<T> sc = FastMath.sinCos(deltaT.multiply(EARTH_AROUND_SUN_PULSATION));
         final T A1 = sc.cos().multiply(C1).add(
                      sc.sin().multiply(C2)).add(C0);
 
         // Get 1rst and 2nd order Legendre polynomials
         final PolynomialFunction firstLegendrePolynomial  = PolynomialsUtils.createLegendrePolynomial(1);
         final PolynomialFunction secondLegendrePolynomial = PolynomialsUtils.createLegendrePolynomial(2);
 
         // Get latitude sinus
         final T sinPhi = FastMath.sin(phi);
 
         // Compute albedo
         return firstLegendrePolynomial.value(sinPhi).multiply(A1).add(
                secondLegendrePolynomial.value(sinPhi).multiply(A2)).add(A0);
 
     }
 
     /** Compute Earth emisivity.
      * Emissivity is used to compute the infrared flux that is emitted by Earth.
      * Its value is in [0;1].
      * @param date the date
      * @param phi the equatorial latitude in rad
      * @return the emissivity in [0;1]
      */
     public double computeEmissivity(final AbsoluteDate date, final double phi) {
 
         // Get duration since coefficient reference epoch
         final double deltaT = date.durationFrom(referenceEpoch);
 
         // Compute 1rst Legendre polynomial coefficient
         final SinCos sc = FastMath.sinCos(EARTH_AROUND_SUN_PULSATION * deltaT);
         final double E1 = K0 +
                           K1 * sc.cos() +
                           K2 * sc.sin();
 
         // Get 1rst and 2nd order Legendre polynomials
         final PolynomialFunction firstLegendrePolynomial  = PolynomialsUtils.createLegendrePolynomial(1);
         final PolynomialFunction secondLegendrePolynomial = PolynomialsUtils.createLegendrePolynomial(2);
 
         // Get latitude sinus
         final double sinPhi = FastMath.sin(phi);
 
         // Compute albedo
         return E0 +
                E1 * firstLegendrePolynomial.value(sinPhi) +
                E2 * secondLegendrePolynomial.value(sinPhi);
 
     }
 
 
     /** Compute Earth emisivity.
      * Emissivity is used to compute the infrared flux that is emitted by Earth.
      * Its value is in [0;1].
      * @param date the date
      * @param phi the equatorial latitude in rad
      * @param <T> extends CalculusFieldElement
      * @return the emissivity in [0;1]
      */
     public <T extends CalculusFieldElement<T>> T computeEmissivity(final FieldAbsoluteDate<T> date, final T phi) {
 
         // Get duration since coefficient reference epoch
         final T deltaT = date.durationFrom(referenceEpoch);
 
         // Compute 1rst Legendre polynomial coeficient
         final FieldSinCos<T> sc = FastMath.sinCos(deltaT.multiply(EARTH_AROUND_SUN_PULSATION));
         final T E1 = sc.cos().multiply(K1).add(
                      sc.sin().multiply(K2)).add(K0);
 
         // Get 1rst and 2nd order Legendre polynomials
         final PolynomialFunction firstLegendrePolynomial  = PolynomialsUtils.createLegendrePolynomial(1);
         final PolynomialFunction secondLegendrePolynomial = PolynomialsUtils.createLegendrePolynomial(2);
 
         // Get latitude sinus
         final T sinPhi = FastMath.sin(phi);
 
         // Compute albedo
         return firstLegendrePolynomial.value(sinPhi).multiply(E1).add(
                secondLegendrePolynomial.value(sinPhi).multiply(E2)).add(E0);
 
     }
 
     /** Compute total solar flux impacting Earth.
      * @param sunPosition the Sun position in an Earth centered frame
      * @return the total solar flux impacting Earth in J/m^3
      */
     public double computeSolarFlux(final Vector3D sunPosition) {
 
         // Compute Earth - Sun distance in UA
         final double earthSunDistance = sunPosition.getNorm() / Constants.JPL_SSD_ASTRONOMICAL_UNIT;
 
         // Compute Solar flux
         return ES_COEFF * Constants.SPEED_OF_LIGHT / (earthSunDistance * earthSunDistance);
     }
 
 
     /** Compute total solar flux impacting Earth.
      * @param sunPosition the Sun position in an Earth centered frame
      * @param <T> extends CalculusFieldElement
      * @return the total solar flux impacting Earth in J/m^3
      */
     public <T extends CalculusFieldElement<T>> T computeSolarFlux(final FieldVector3D<T> sunPosition) {
 
         // Compute Earth - Sun distance in UA
         final T earthSunDistance = sunPosition.getNorm().divide(Constants.JPL_SSD_ASTRONOMICAL_UNIT);
 
         // Compute Solar flux
         return earthSunDistance.multiply(earthSunDistance).reciprocal().multiply(ES_COEFF * Constants.SPEED_OF_LIGHT);
     }
 
 
     /** Compute elementary rediffused flux on satellite.
      * @param state the current spacecraft state
      * @param elementCenter the position of the considered area center
      * @param sunPosition the position of the Sun in the spacecraft frame
      * @param elementArea the area of the current element
      * @return the rediffused flux from considered element on the spacecraft
      */
     public Vector3D computeElementaryFlux(final SpacecraftState state,
                                           final Vector3D elementCenter,
                                           final Vector3D sunPosition,
                                           final double elementArea) {
 
         // Get satellite position
         final Vector3D satellitePosition = state.getPosition();
 
         // Get current date
         final AbsoluteDate date = state.getDate();
 
         // Get solar flux impacting Earth
         final double solarFlux = computeSolarFlux(sunPosition);
 
         // Get satellite viewing angle as seen from current elementary area
         final double centerNorm = elementCenter.getNorm();
         final double cosAlpha   = Vector3D.dotProduct(elementCenter, satellitePosition) /
                                   (centerNorm * satellitePosition.getNorm());
 
         // Check that satellite sees the current area
         if (cosAlpha > 0) {
 
             // Get current elementary area center latitude
             final double currentLatitude = elementCenter.getDelta();
 
             // Compute Earth emissivity value
             final double e = computeEmissivity(date, currentLatitude);
 
             // Initialize albedo
             double a = 0.0;
 
             // Check if elementary area is in daylight
             final double cosSunAngle = Vector3D.dotProduct(elementCenter, sunPosition) /
                                        (centerNorm * sunPosition.getNorm());
 
             if (cosSunAngle > 0) {
                 // Elementary area is in daylight, compute albedo value
                 a = computeAlbedo(date, currentLatitude);
             }
 
             // Compute elementary area contribution to rediffused flux
             final double albedoAndIR = a * solarFlux * cosSunAngle + e * solarFlux * 0.25;
 
             // Compute elementary area - satellite vector and distance
             final Vector3D r = satellitePosition.subtract(elementCenter);
             final double rNorm = r.getNorm();
 
             // Compute attenuated projected elementary area vector
             final Vector3D projectedAreaVector = r.scalarMultiply(elementArea * cosAlpha /
                                                                  (FastMath.PI * rNorm * rNorm * rNorm));
 
             // Compute elementary radiation flux from current elementary area
             return projectedAreaVector.scalarMultiply(albedoAndIR / Constants.SPEED_OF_LIGHT);
 
         } else {
 
             // Spacecraft does not see the elementary area
             return new Vector3D(0.0, 0.0, 0.0);
         }
 
     }
 
 
     /** Compute elementary rediffused flux on satellite.
      * @param state the current spacecraft state
      * @param elementCenter the position of the considered area center
      * @param sunPosition the position of the Sun in the spacecraft frame
      * @param elementArea the area of the current element
      * @param <T> extends CalculusFieldElement
      * @return the rediffused flux from considered element on the spacecraft
      */
     public <T extends CalculusFieldElement<T>> FieldVector3D<T> computeElementaryFlux(final FieldSpacecraftState<T> state,
                                                                                       final FieldVector3D<T> elementCenter,
                                                                                       final FieldVector3D<T> sunPosition,
                                                                                       final T elementArea) {
 
         // Get satellite position
         final FieldVector3D<T> satellitePosition = state.getPosition();
 
         // Get current date
         final FieldAbsoluteDate<T> date = state.getDate();
 
         // Get zero
         final T zero = date.getField().getZero();
 
         // Get solar flux impacting Earth
         final T solarFlux = computeSolarFlux(sunPosition);
 
         // Get satellite viewing angle as seen from current elementary area
         final T centerNorm = elementCenter.getNorm();
         final T cosAlpha   = FieldVector3D.dotProduct(elementCenter, satellitePosition).
                              divide(centerNorm.multiply(satellitePosition.getNorm()));
 
         // Check that satellite sees the current area
         if (cosAlpha.getReal() > 0) {
 
             // Get current elementary area center latitude
             final T currentLatitude = elementCenter.getDelta();
 
             // Compute Earth emissivity value
             final T e = computeEmissivity(date, currentLatitude);
 
             // Initialize albedo
             T a = zero;
 
             // Check if elementary area is in daylight
             final T cosSunAngle = FieldVector3D.dotProduct(elementCenter, sunPosition).
                                   divide(centerNorm.multiply(sunPosition.getNorm()));
 
             if (cosSunAngle.getReal() > 0) {
                 // Elementary area is in daylight, compute albedo value
                 a = computeAlbedo(date, currentLatitude);
             }
 
             // Compute elementary area contribution to rediffused flux
             final T albedoAndIR = a.multiply(solarFlux).multiply(cosSunAngle).
                                   add(e.multiply(solarFlux).multiply(0.25));
 
             // Compute elementary area - satellite vector and distance
             final FieldVector3D<T> r = satellitePosition.subtract(elementCenter);
             final T rNorm = r.getNorm();
 
             // Compute attenuated projected elementary area vector
             final FieldVector3D<T> projectedAreaVector = r.scalarMultiply(elementArea.multiply(cosAlpha).
                                                                           divide(rNorm.square().multiply(rNorm).multiply(zero.getPi())));
 
             // Compute elementary radiation flux from current elementary area
             return projectedAreaVector.scalarMultiply(albedoAndIR.divide(Constants.SPEED_OF_LIGHT));
 
         } else {
 
             // Spacecraft does not see the elementary area
             return new FieldVector3D<>(zero, zero, zero);
         }
 
     }
 
 }