ConicallyShadowedLightFluxModel.java

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