AbstractRadiationForceModel.java
- /* Copyright 2002-2023 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.forces.ForceModel;
- 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 ForceModel {
- /** 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 boolean dependsOnPositionOnly() {
- return false;
- }
- /** {@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), ForceModel.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), ForceModel.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);
- }
- }