AbstractRadiationForceModel

/* Copyright 2002-2024 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.lang.reflect.Array;
import java.util.ArrayList;
import java.util.Collections;
import java.util.List;
import java.util.stream.Stream;

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.ode.events.Action;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.MathArrays;
import org.orekit.bodies.OneAxisEllipsoid;
import org.orekit.frames.Frame;
import org.orekit.propagation.events.EclipseDetector;
import org.orekit.propagation.events.EventDetector;
import org.orekit.propagation.events.FieldEclipseDetector;
import org.orekit.propagation.events.FieldEventDetector;
import org.orekit.utils.Constants;
import org.orekit.utils.ExtendedPVCoordinatesProvider;
import org.orekit.utils.ExtendedPVCoordinatesProviderAdapter;
import org.orekit.utils.OccultationEngine;

/**
 * Base class for radiation force models.
 * @see SolarRadiationPressure
 * @see ECOM2
 * @since 10.2
 */
public abstract class AbstractRadiationForceModel implements RadiationForceModel {

    /** Margin to force recompute lighting ratio derivatives when we are really inside penumbra. */
    private static final double ANGULAR_MARGIN = 1.0e-10;

    /** Max check interval for eclipse detectors. */
    private static final double ECLIPSE_MAX_CHECK = 60.0;

    /** Threshold for eclipse detectors. */
    private static final double ECLIPSE_THRESHOLD = 1.0e-7; // this is consistent with ANGULAR_MARGIN = 10⁻¹⁰ rad for LEO

    /** Flatness for spherical bodies. */
    private static final double SPHERICAL_BODY_FLATNESS = 0.0;

    /** Prefix for occulting bodies frames names. */
    private static final String OCCULTING_PREFIX = "occulting-";

    /** Occulting bodies (central body is always the first one).
     * @since 12.0
     */
    private final List<OccultationEngine> occultingBodies;

    /**
     * Default constructor.
     * Only central body is considered.
     * @param sun Sun model
     * @param centralBody central body shape model (for umbra/penumbra computation)
     * @since 12.0
     */
    protected AbstractRadiationForceModel(final ExtendedPVCoordinatesProvider sun, final OneAxisEllipsoid centralBody) {
        // in most cases, there will be only Earth, sometimes also Moon so an initial capacity of 2 is appropriate
        occultingBodies = new ArrayList<>(2);
        occultingBodies.add(new OccultationEngine(sun, Constants.SUN_RADIUS, centralBody));
    }

    /** {@inheritDoc} */
    @Override
    public Stream<EventDetector> getEventDetectors() {
        final EventDetector[] detectors = new EventDetector[2 * occultingBodies.size()];
        for (int i = 0; i < occultingBodies.size(); ++i) {
            final OccultationEngine occulting = occultingBodies.get(i);
            detectors[2 * i]     = new EclipseDetector(occulting).
                                   withUmbra().
                                   withMargin(-ANGULAR_MARGIN).
                                   withMaxCheck(ECLIPSE_MAX_CHECK).
                                   withThreshold(ECLIPSE_THRESHOLD).
                                   withHandler((state, detector, increasing) -> Action.RESET_DERIVATIVES);
            detectors[2 * i + 1] = new EclipseDetector(occulting).
                                   withPenumbra().
                                   withMargin(ANGULAR_MARGIN).
                                   withMaxCheck(ECLIPSE_MAX_CHECK).
                                   withThreshold(ECLIPSE_THRESHOLD).
                                   withHandler((state, detector, increasing) -> Action.RESET_DERIVATIVES);
        }
        // Fusion between Date detector for parameter driver span change and
        // Detector for umbra / penumbra events
        return Stream.concat(Stream.of(detectors), RadiationForceModel.super.getEventDetectors());
    }

    /** {@inheritDoc} */
    @Override
    public <T extends CalculusFieldElement<T>> Stream<FieldEventDetector<T>> getFieldEventDetectors(final Field<T> field) {
        final T zero = field.getZero();
        @SuppressWarnings("unchecked")
        final FieldEventDetector<T>[] detectors = (FieldEventDetector<T>[]) Array.newInstance(FieldEventDetector.class,
                                                                                              2 * occultingBodies.size());
        for (int i = 0; i < occultingBodies.size(); ++i) {
            final OccultationEngine occulting = occultingBodies.get(i);
            detectors[2 * i]     = new FieldEclipseDetector<>(field, occulting).
                                   withUmbra().
                                   withMargin(zero.newInstance(-ANGULAR_MARGIN)).
                                   withMaxCheck(ECLIPSE_MAX_CHECK).
                                   withThreshold(zero.newInstance(ECLIPSE_THRESHOLD)).
                                   withHandler((state, detector, increasing) -> Action.RESET_DERIVATIVES);
            detectors[2 * i + 1] = new FieldEclipseDetector<>(field, occulting).
                                   withPenumbra().
                                   withMargin(zero.newInstance(ANGULAR_MARGIN)).
                                   withMaxCheck(ECLIPSE_MAX_CHECK).
                                   withThreshold(zero.newInstance(ECLIPSE_THRESHOLD)).
                                   withHandler((state, detector, increasing) -> Action.RESET_DERIVATIVES);
        }
        return Stream.concat(Stream.of(detectors), RadiationForceModel.super.getFieldEventDetectors(field));
    }

    /**
     * Get the useful angles for eclipse computation.
     * @param position the satellite's position in the selected frame
     * @param occultingPosition Oculting body position in the selected frame
     * @param occultingRadius Occulting body mean radius
     * @param occultedPosition Occulted body position in the selected frame
     * @param occultedRadius Occulted body mean radius
     * @return the 3 angles {(satOcculting, satOcculted), Occulting body apparent radius, Occulted body apparent radius}
     */
    protected double[] getGeneralEclipseAngles(final Vector3D position, final Vector3D occultingPosition, final double occultingRadius,
                                               final Vector3D occultedPosition, final double occultedRadius) {
        final double[] angle = new double[3];

        final Vector3D satOccultedVector = occultedPosition.subtract(position);
        final Vector3D satOccultingVector = occultingPosition.subtract(position);

        // Sat-Occulted / Sat-Occulting angle
        angle[0] = Vector3D.angle(satOccultedVector, satOccultingVector);

        // Occulting body apparent radius
        angle[1] = FastMath.asin(occultingRadius / satOccultingVector.getNorm());

        // Occulted body apparent radius
        angle[2] = FastMath.asin(occultedRadius / satOccultedVector.getNorm());

        return angle;
    }

    /**
     * Get the useful angles for eclipse computation.
     * @param occultingPosition Oculting body position in the selected frame
     * @param occultingRadius Occulting body mean radius
     * @param occultedPosition Occulted body position in the selected frame
     * @param occultedRadius Occulted body mean radius
     * @param position the satellite's position in the selected frame
     * @param <T> extends RealFieldElement
     * @return the 3 angles {(satOcculting, satOcculted), Occulting body apparent radius, Occulted body apparent radius}
     */
    protected <T extends CalculusFieldElement<T>> T[] getGeneralEclipseAngles(final FieldVector3D<T> position,
                                                                              final FieldVector3D<T> occultingPosition, final T occultingRadius,
                                                                              final FieldVector3D<T> occultedPosition, final T occultedRadius) {
        final T[] angle = MathArrays.buildArray(position.getX().getField(), 3);

        final FieldVector3D<T> satOccultedVector = occultedPosition.subtract(position);
        final FieldVector3D<T> satOccultingVector = occultingPosition.subtract(position);

        // Sat-Occulted / Sat-Occulting angle
        angle[0] = FieldVector3D.angle(satOccultedVector, satOccultingVector);

        // Occulting body apparent radius
        angle[1] = occultingRadius.divide(satOccultingVector.getNorm()).asin();

        // Occulted body apparent radius
        angle[2] = occultedRadius.divide(satOccultedVector.getNorm()).asin();

        return angle;
    }

    /**
     * Add a new occulting body.
     * <p>
     * Central body is already considered, it shall not be added this way.
     * </p>
     * @param provider body PV provider
     * @param radius body mean radius
     * @see #addOccultingBody(OneAxisEllipsoid)
     */
    public void addOccultingBody(final ExtendedPVCoordinatesProvider provider, final double radius) {

        // as parent frame for occulting body frame,
        // we select the first inertial frame in central body hierarchy
        Frame parent = occultingBodies.get(0).getOcculting().getBodyFrame();
        while (!parent.isPseudoInertial()) {
            parent = parent.getParent();
        }

        // as the occulting body will be spherical, we can use an inertially oriented body frame
        final Frame inertiallyOrientedBodyFrame =
                        new ExtendedPVCoordinatesProviderAdapter(parent,
                                                                 provider,
                                                                 OCCULTING_PREFIX + occultingBodies.size());

        // create the spherical occulting body
        final OneAxisEllipsoid sphericalOccultingBody =
                        new OneAxisEllipsoid(radius, SPHERICAL_BODY_FLATNESS, inertiallyOrientedBodyFrame);

        addOccultingBody(sphericalOccultingBody);

    }

    /**
     * Add a new occulting body.
     * <p>
     * Central body is already considered, it shall not be added this way.
     * </p>
     * @param occulting occulting body to add
     * @see #addOccultingBody(ExtendedPVCoordinatesProvider, double)
     * @since 12.0
     */
    public void addOccultingBody(final OneAxisEllipsoid occulting) {

        // retrieve Sun from the central occulting body engine
        final OccultationEngine central = occultingBodies.get(0);

        // create a new occultation engine for the new occulting body
        final OccultationEngine additional =
                        new OccultationEngine(central.getOcculted(), central.getOccultedRadius(), occulting);

        // add it to the list
        occultingBodies.add(additional);

    }

    /**
     * Get all occulting bodies to consider.
     * <p>
     * The list always contains at least one element: the central body
     * which is always the first one in the list.
     * </p>
     * @return immutable list of all occulting bodies to consider, including the central body
     * @since 12.0
     */
    public List<OccultationEngine> getOccultingBodies() {
        return Collections.unmodifiableList(occultingBodies);
    }

}