/* Copyright 2022-2025 Romain Serra
    * 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 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.orekit.frames.Frame;
  import org.orekit.propagation.FieldSpacecraftState;
  import org.orekit.propagation.SpacecraftState;
  import org.orekit.propagation.events.intervals.AdaptableInterval;
  import org.orekit.propagation.events.EventDetector;
  import org.orekit.propagation.events.EventDetectionSettings;
  import org.orekit.propagation.events.intervals.FieldAdaptableInterval;
  import org.orekit.propagation.events.FieldEventDetector;
  import org.orekit.propagation.events.FieldEventDetectionSettings;
  import org.orekit.propagation.events.handlers.EventHandler;
  import org.orekit.propagation.events.handlers.FieldEventHandler;
  import org.orekit.propagation.events.handlers.FieldResetDerivativesOnEvent;
  import org.orekit.propagation.events.handlers.ResetDerivativesOnEvent;
  import org.orekit.time.AbsoluteDate;
  import org.orekit.time.FieldAbsoluteDate;
  import org.orekit.utils.ExtendedPositionProvider;
  
  import java.util.ArrayList;
  import java.util.List;
  
  /**
   * Class defining a flux model from a single occulted body, casting a shadow on a spherical occulting body.
   * It cannot model oblate bodies or multiple occulting objects (for this, see {@link SolarRadiationPressure}).
   *
   * @author Romain Serra
   * @see AbstractSolarLightFluxModel
   * @see LightFluxModel
   * @see "Montenbruck, Oliver, and Gill, Eberhard. Satellite orbits : models, methods, and
   *  * applications. Berlin New York: Springer, 2000."
   * @since 12.2
   */
  public class ConicallyShadowedLightFluxModel extends AbstractSolarLightFluxModel {
  
      /**
       * Default max. check interval for eclipse detection.
       */
      private static final double CONICAL_ECLIPSE_MAX_CHECK = 60;
  
      /**
       * Default threshold for eclipse detection.
       */
      private static final double CONICAL_ECLIPSE_THRESHOLD = 1e-7;
  
      /** Occulted body radius. */
      private final double occultedBodyRadius;
  
      /** Cached date. */
      private AbsoluteDate lastDate;
  
      /** Cached frame. */
      private Frame propagationFrame;
  
      /** Cached position. */
      private Vector3D lastPosition;
  
      /**
       * Constructor.
       * @param kRef reference flux
       * @param occultedBodyRadius radius of occulted body (light source)
       * @param occultedBody position provider for light source
       * @param occultingBodyRadius radius of central, occulting body
       * @param eventDetectionSettings user-defined detection settings for eclipses (if ill-tuned, events might be missed or performance might drop)
       */
      public ConicallyShadowedLightFluxModel(final double kRef, final double occultedBodyRadius,
                                             final ExtendedPositionProvider occultedBody,
                                             final double occultingBodyRadius, final EventDetectionSettings eventDetectionSettings) {
          super(kRef, occultedBody, occultingBodyRadius, eventDetectionSettings);
          this.occultedBodyRadius = occultedBodyRadius;
      }
  
      /**
       * Constructor with default event detection settings.
       * @param kRef reference flux
       * @param occultedBodyRadius radius of occulted body (light source)
       * @param occultedBody position provider for light source
       * @param occultingBodyRadius radius of central, occulting body
      */
     public ConicallyShadowedLightFluxModel(final double kRef, final double occultedBodyRadius,
                                            final ExtendedPositionProvider occultedBody, final double occultingBodyRadius) {
         this(kRef, occultedBodyRadius, occultedBody, occultingBodyRadius, getDefaultEclipseDetectionSettings());
     }
 
     /**
      * Constructor with default value for reference flux.
      * @param occultedBodyRadius radius of occulted body (light source)
      * @param occultedBody position provider for light source
      * @param occultingBodyRadius radius of central, occulting body
      */
     public ConicallyShadowedLightFluxModel(final double occultedBodyRadius, final ExtendedPositionProvider occultedBody,
                                            final double occultingBodyRadius) {
         super(occultedBody, occultingBodyRadius, getDefaultEclipseDetectionSettings());
         this.occultedBodyRadius = occultedBodyRadius;
     }
 
     /**
      * Define default detection settings for eclipses.
      * @return default settings
      */
     public static EventDetectionSettings getDefaultEclipseDetectionSettings() {
         return new EventDetectionSettings(CONICAL_ECLIPSE_MAX_CHECK, CONICAL_ECLIPSE_THRESHOLD,
                 EventDetectionSettings.DEFAULT_MAX_ITER);
     }
 
     /** {@inheritDoc} */
     @Override
     public void init(final SpacecraftState initialState, final AbsoluteDate targetDate) {
         super.init(initialState, targetDate);
         lastDate = initialState.getDate();
         propagationFrame = initialState.getFrame();
         lastPosition = getOccultedBodyPosition(lastDate, propagationFrame);
     }
 
     /** {@inheritDoc} */
     @Override
     public <T extends CalculusFieldElement<T>> void init(final FieldSpacecraftState<T> initialState,
                                                          final FieldAbsoluteDate<T> targetDate) {
         super.init(initialState, targetDate);
         lastDate = initialState.getDate().toAbsoluteDate();
         propagationFrame = initialState.getFrame();
         lastPosition = getOccultedBodyPosition(initialState.getDate(), propagationFrame).toVector3D();
     }
 
     /**
      * Get occulted body position using cache.
      * @param date date
      * @return occulted body position
      */
     private Vector3D getOccultedBodyPosition(final AbsoluteDate date) {
         if (!lastDate.isEqualTo(date)) {
             lastPosition = getOccultedBodyPosition(date, propagationFrame);
             lastDate = date;
         }
         return lastPosition;
     }
 
     /**
      * Get occulted body position using cache (non-Field and no derivatives case).
      * @param fieldDate date
      * @param <T> field type
      * @return occulted body position
      */
     private <T extends CalculusFieldElement<T>> FieldVector3D<T> getOccultedBodyPosition(final FieldAbsoluteDate<T> fieldDate) {
         if (fieldDate.hasZeroField()) {
             final AbsoluteDate date = fieldDate.toAbsoluteDate();
             if (!lastDate.isEqualTo(date)) {
                 lastPosition = getOccultedBodyPosition(date, propagationFrame);
                 lastDate = date;
             }
             return new FieldVector3D<>(fieldDate.getField(), lastPosition);
         } else {
             return getOccultedBodyPosition(fieldDate, propagationFrame);
         }
     }
 
     /** {@inheritDoc} */
     @Override
     protected double getLightingRatio(final Vector3D position, final Vector3D occultedBodyPosition) {
         final double distanceSun = occultedBodyPosition.getNorm();
         final double squaredDistance = position.getNormSq();
         final Vector3D occultedBodyDirection = occultedBodyPosition.normalize();
         final double s0 = -position.dotProduct(occultedBodyDirection);
         if (s0 > 0.0) {
             final double l = FastMath.sqrt(squaredDistance - s0 * s0);
             final double sinf2 = (occultedBodyRadius - getOccultingBodyRadius()) / distanceSun;
             final double l2 = (s0 * sinf2 - getOccultingBodyRadius()) / FastMath.sqrt(1.0 - sinf2 * sinf2);
             if (FastMath.abs(l2) - l >= 0.0) { // umbra
                 return 0.;
             }
             final double sinf1 = (occultedBodyRadius + getOccultingBodyRadius()) / distanceSun;
             final double l1 = (s0 * sinf1 + getOccultingBodyRadius()) / FastMath.sqrt(1.0 - sinf1 * sinf1);
             if (l1 - l > 0.0) { // penumbra
                 final Vector3D relativePosition = occultedBodyPosition.subtract(position);
                 final double relativeDistance = relativePosition.getNorm();
                 final double a = FastMath.asin(occultedBodyRadius / relativeDistance);
                 final double a2 = a * a;
                 final double r = FastMath.sqrt(squaredDistance);
                 final double b = FastMath.asin(getOccultingBodyRadius() / r);
                 final double c = FastMath.acos(-(relativePosition.dotProduct(position)) / (r * relativeDistance));
                 final double x = (c * c + a2 - b * b) / (2 * c);
                 final double y = FastMath.sqrt(FastMath.max(0., a2 - x * x));
                 final double arcCosXOverA = FastMath.acos(FastMath.max(-1, x / a));
                 final double intermediate = (arcCosXOverA + (b * b * FastMath.acos((c - x) / b) - c * y) / a2) / FastMath.PI;
                 return 1. - intermediate;
             }
         }
         return 1.;
     }
 
     /** {@inheritDoc} */
     @Override
     protected <T extends CalculusFieldElement<T>> T getLightingRatio(final FieldVector3D<T> position,
                                                                      final FieldVector3D<T> occultedBodyPosition) {
         final Field<T> field = position.getX().getField();
         final T distanceSun = occultedBodyPosition.getNorm();
         final T squaredDistance = position.getNormSq();
         final FieldVector3D<T> occultedBodyDirection = occultedBodyPosition.normalize();
         final T s0 = position.dotProduct(occultedBodyDirection).negate();
         if (s0.getReal() > 0.0) {
             final T reciprocalDistanceSun = distanceSun.reciprocal();
             final T sinf2 = reciprocalDistanceSun.multiply(occultedBodyRadius - getOccultingBodyRadius());
             final T l2 = (s0.multiply(sinf2).subtract(getOccultingBodyRadius())).divide(FastMath.sqrt(sinf2.square().negate().add(1)));
             final T l = FastMath.sqrt(squaredDistance.subtract(s0.square()));
             if (FastMath.abs(l2).subtract(l).getReal() >= 0.0) { // umbra
                 return field.getZero();
             }
             final T sinf1 = reciprocalDistanceSun.multiply(occultedBodyRadius + getOccultingBodyRadius());
             final T l1 = (s0.multiply(sinf1).add(getOccultingBodyRadius())).divide(FastMath.sqrt(sinf1.square().negate().add(1)));
             if (l1.subtract(l).getReal() > 0.0) { // penumbra
                 final FieldVector3D<T> relativePosition = occultedBodyPosition.subtract(position);
                 final T relativeDistance = relativePosition.getNorm();
                 final T a = FastMath.asin(relativeDistance.reciprocal().multiply(occultedBodyRadius));
                 final T a2 = a.square();
                 final T r = FastMath.sqrt(squaredDistance);
                 final T b = FastMath.asin(r.reciprocal().multiply(getOccultingBodyRadius()));
                 final T b2 = b.square();
                 final T c = FastMath.acos((relativePosition.dotProduct(position).negate()).divide(r.multiply(relativeDistance)));
                 final T x = (c.square().add(a2).subtract(b2)).divide(c.multiply(2));
                 final T x2 = x.square();
                 final T y = (a2.getReal() - x2.getReal() <= 0) ? s0.getField().getZero() : FastMath.sqrt(a2.subtract(x2));
                 final T arcCosXOverA = (x.getReal() / a.getReal() < -1) ? s0.getPi().negate() : FastMath.acos(x.divide(a));
                 final T intermediate = arcCosXOverA.add(((b2.multiply(FastMath.acos((c.subtract(x)).divide(b))))
                         .subtract(c.multiply(y))).divide(a2));
                 return intermediate.divide(-FastMath.PI).add(1);
             }
         }
         return field.getOne();
     }
 
     /** {@inheritDoc} */
     @Override
     public List<EventDetector> getEclipseConditionsDetector() {
         final List<EventDetector> detectors = new ArrayList<>();
         detectors.add(createUmbraEclipseDetector());
         detectors.add(createPenumbraEclipseDetector());
         return detectors;
     }
 
     /**
      * Method to create a new umbra detector.
      * @return detector
      */
     private InternalEclipseDetector createUmbraEclipseDetector() {
         return new InternalEclipseDetector() {
             @Override
             public double g(final SpacecraftState s) {
                 final Vector3D position = s.getPosition();
                 final Vector3D occultedBodyPosition = getOccultedBodyPosition(s.getDate());
                 final Vector3D occultedBodyDirection = occultedBodyPosition.normalize();
                 final double s0 = -position.dotProduct(occultedBodyDirection);
                 final double distanceSun = occultedBodyPosition.getNorm();
                 final double squaredDistance = position.getNormSq();
                 final double sinf2 = (occultedBodyRadius - getOccultingBodyRadius()) / distanceSun;
                 final double l = FastMath.sqrt(squaredDistance - s0 * s0);
                 final double l2 = (s0 * sinf2 - getOccultingBodyRadius()) / FastMath.sqrt(1.0 - sinf2 * sinf2);
                 return FastMath.abs(l2) / l - 1.;
             }
         };
     }
 
     /**
      * Method to create a new penumbra detector.
      * @return detector
      */
     private InternalEclipseDetector createPenumbraEclipseDetector() {
         return new InternalEclipseDetector() {
             @Override
             public double g(final SpacecraftState s) {
                 final Vector3D position = s.getPosition();
                 final Vector3D occultedBodyPosition = getOccultedBodyPosition(s.getDate());
                 final Vector3D occultedBodyDirection = occultedBodyPosition.normalize();
                 final double s0 = -position.dotProduct(occultedBodyDirection);
                 final double distanceSun = occultedBodyPosition.getNorm();
                 final double squaredDistance = position.getNormSq();
                 final double l = FastMath.sqrt(squaredDistance - s0 * s0);
                 final double sinf1 = (occultedBodyRadius + getOccultingBodyRadius()) / distanceSun;
                 final double l1 = (s0 * sinf1 + getOccultingBodyRadius()) / FastMath.sqrt(1.0 - sinf1 * sinf1);
                 return l1 / l - 1.;
             }
         };
     }
 
     /**
      * Internal class for event detector.
      */
     private abstract class InternalEclipseDetector implements EventDetector {
         /** Event handler. */
         private final ResetDerivativesOnEvent handler;
 
         /**
          * Constructor.
          */
         InternalEclipseDetector() {
             this.handler = new ResetDerivativesOnEvent();
         }
 
         @Override
         public EventDetectionSettings getDetectionSettings() {
             return getEventDetectionSettings();
         }
 
         @Override
         public double getThreshold() {
             return getDetectionSettings().getThreshold();
         }
 
         @Override
         public AdaptableInterval getMaxCheckInterval() {
             return getDetectionSettings().getMaxCheckInterval();
         }
 
         @Override
         public int getMaxIterationCount() {
             return getDetectionSettings().getMaxIterationCount();
         }
 
         @Override
         public EventHandler getHandler() {
             return handler;
         }
     }
 
     /** {@inheritDoc} */
     @Override
     public <T extends CalculusFieldElement<T>> List<FieldEventDetector<T>> getFieldEclipseConditionsDetector(final Field<T> field) {
         final List<FieldEventDetector<T>> detectors = new ArrayList<>();
         final FieldEventDetectionSettings<T> detectionSettings = new FieldEventDetectionSettings<>(field,
             getEventDetectionSettings());
         detectors.add(createFieldUmbraEclipseDetector(detectionSettings));
         detectors.add(createFieldPenumbraEclipseDetector(detectionSettings));
         return detectors;
     }
 
     /**
      * Method to create a new umbra detector. Field version.
      * @param detectionSettings non-Field detection settings
      * @param <T> field type
      * @return detector
      */
     private <T extends CalculusFieldElement<T>> FieldInternalEclipseDetector<T> createFieldUmbraEclipseDetector(final FieldEventDetectionSettings<T> detectionSettings) {
         return new FieldInternalEclipseDetector<T>(detectionSettings) {
             @Override
             public T g(final FieldSpacecraftState<T> s) {
                 final FieldVector3D<T> position = s.getPosition();
                 final FieldVector3D<T> occultedBodyPosition = getOccultedBodyPosition(s.getDate());
                 final FieldVector3D<T> occultedBodyDirection = occultedBodyPosition.normalize();
                 final T s0 = position.dotProduct(occultedBodyDirection).negate();
                 final T distanceSun = occultedBodyPosition.getNorm();
                 final T squaredDistance = position.getNormSq();
                 final T reciprocalDistanceSun = distanceSun.reciprocal();
                 final T sinf2 = reciprocalDistanceSun.multiply(occultedBodyRadius - getOccultingBodyRadius());
                 final T l2 = (s0.multiply(sinf2).subtract(getOccultingBodyRadius())).divide(FastMath.sqrt(sinf2.square().negate().add(1)));
                 final T l = FastMath.sqrt(squaredDistance.subtract(s0.square()));
                 return FastMath.abs(l2).divide(l).subtract(1);
             }
         };
     }
 
     /**
      * Method to create a new penumbra detector. Field version.
      * @param detectionSettings non-Field detection settings
      * @param <T> field type
      * @return detector
      */
     private <T extends CalculusFieldElement<T>> FieldInternalEclipseDetector<T> createFieldPenumbraEclipseDetector(final FieldEventDetectionSettings<T> detectionSettings) {
         return new FieldInternalEclipseDetector<T>(detectionSettings) {
             @Override
             public T g(final FieldSpacecraftState<T> s) {
                 final FieldVector3D<T> position = s.getPosition();
                 final FieldVector3D<T> occultedBodyPosition = getOccultedBodyPosition(s.getDate());
                 final FieldVector3D<T> occultedBodyDirection = occultedBodyPosition.normalize();
                 final T s0 = position.dotProduct(occultedBodyDirection).negate();
                 final T distanceSun = occultedBodyPosition.getNorm();
                 final T squaredDistance = position.getNormSq();
                 final T reciprocalDistanceSun = distanceSun.reciprocal();
                 final T sinf1 = reciprocalDistanceSun.multiply(occultedBodyRadius + getOccultingBodyRadius());
                 final T l1 = (s0.multiply(sinf1).add(getOccultingBodyRadius())).divide(FastMath.sqrt(sinf1.square().negate().add(1)));
                 final T l = FastMath.sqrt(squaredDistance.subtract(s0.square()));
                 return l1.divide(l).subtract(1);
             }
         };
     }
 
     /**
      * Internal class for event detector.
      */
     private abstract static class FieldInternalEclipseDetector<T extends CalculusFieldElement<T>> implements FieldEventDetector<T> {
         /** Event handler. */
         private final FieldResetDerivativesOnEvent<T> handler;
 
         /** Detection settings. */
         private final FieldEventDetectionSettings<T> fieldEventDetectionSettings;
 
         /**
          * Constructor.
          * @param fieldEventDetectionSettings detection settings
          */
         FieldInternalEclipseDetector(final FieldEventDetectionSettings<T> fieldEventDetectionSettings) {
             this.handler = new FieldResetDerivativesOnEvent<>();
             this.fieldEventDetectionSettings = fieldEventDetectionSettings;
         }
 
         @Override
         public FieldEventDetectionSettings<T> getDetectionSettings() {
             return fieldEventDetectionSettings;
         }
 
         @Override
         public T getThreshold() {
             return getDetectionSettings().getThreshold();
         }
 
         @Override
         public FieldAdaptableInterval<T> getMaxCheckInterval() {
             return getDetectionSettings().getMaxCheckInterval();
         }
 
         @Override
         public int getMaxIterationCount() {
             return getDetectionSettings().getMaxIterationCount();
         }
 
         @Override
         public FieldEventHandler<T> getHandler() {
             return handler;
         }
     }
 }