AbstractRadiationForceModel.java

  1. /* Copyright 2002-2025 CS GROUP
  2.  * Licensed to CS GROUP (CS) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * CS licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *   http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.orekit.forces.radiation;

  19. import java.lang.reflect.Array;
  20. import java.util.ArrayList;
  21. import java.util.Collections;
  22. import java.util.List;
  23. import java.util.stream.Stream;

  25. import org.hipparchus.CalculusFieldElement;
  26. import org.hipparchus.Field;
  27. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  28. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  29. import org.hipparchus.util.FastMath;
  30. import org.hipparchus.util.MathArrays;
  31. import org.orekit.bodies.OneAxisEllipsoid;
  32. import org.orekit.frames.Frame;
  33. import org.orekit.propagation.events.EclipseDetector;
  34. import org.orekit.propagation.events.EventDetector;
  35. import org.orekit.propagation.events.EventDetectionSettings;
  36. import org.orekit.propagation.events.FieldEclipseDetector;
  37. import org.orekit.propagation.events.FieldEventDetector;
  38. import org.orekit.propagation.events.FieldEventDetectionSettings;
  39. import org.orekit.propagation.events.handlers.FieldResetDerivativesOnEvent;
  40. import org.orekit.propagation.events.handlers.ResetDerivativesOnEvent;
  41. import org.orekit.utils.Constants;
  42. import org.orekit.utils.ExtendedPositionProviderAdapter;
  43. import org.orekit.utils.ExtendedPositionProvider;
  44. import org.orekit.utils.OccultationEngine;

  46. /**
  47.  * Base class for radiation force models.
  48.  * @see SolarRadiationPressure
  49.  * @see ECOM2
  50.  * @since 10.2
  51.  */
  52. public abstract class AbstractRadiationForceModel implements RadiationForceModel {

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

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

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

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

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

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

  74.     /** Eclipse detection settings. */
  75.     private final EventDetectionSettings eclipseDetectionSettings;

  77.     /**
  78.      * Default constructor.
  79.      * Only central body is considered.
  80.      * @param sun Sun model
  81.      * @param centralBody central body shape model (for umbra/penumbra computation)
  82.      * @param eclipseDetectionSettings eclipse detection settings (warning: poor settings will miss eclipses)
  83.      * @since 13.0
  84.      */
  85.     protected AbstractRadiationForceModel(final ExtendedPositionProvider sun, final OneAxisEllipsoid centralBody,
  86.                                           final EventDetectionSettings eclipseDetectionSettings) {
  87.         // in most cases, there will be only Earth, sometimes also Moon so an initial capacity of 2 is appropriate
  88.         occultingBodies = new ArrayList<>(2);
  89.         occultingBodies.add(new OccultationEngine(sun, Constants.SUN_RADIUS, centralBody));
  90.         this.eclipseDetectionSettings = eclipseDetectionSettings;
  91.     }

  93.     /** {@inheritDoc} */
  94.     @Override
  95.     public Stream<EventDetector> getEventDetectors() {
  96.         final EventDetector[] detectors = new EventDetector[2 * occultingBodies.size()];
  97.         for (int i = 0; i < occultingBodies.size(); ++i) {
  98.             final OccultationEngine occulting = occultingBodies.get(i);
  99.             detectors[2 * i]     = new EclipseDetector(occulting).
  100.                                    withUmbra().
  101.                                    withMargin(-ANGULAR_MARGIN).
  102.                                    withDetectionSettings(eclipseDetectionSettings).
  103.                                    withHandler(new ResetDerivativesOnEvent());
  104.             detectors[2 * i + 1] = new EclipseDetector(occulting).
  105.                                    withPenumbra().
  106.                                    withMargin(ANGULAR_MARGIN).
  107.                                    withDetectionSettings(eclipseDetectionSettings).
  108.                                    withHandler(new ResetDerivativesOnEvent());
  109.         }
  110.         // Fusion between Date detector for parameter driver span change and
  111.         // Detector for umbra / penumbra events
  112.         return Stream.concat(Stream.of(detectors), RadiationForceModel.super.getEventDetectors());
  113.     }

  115.     /**
  116.      * Get the default eclipse detection settings.
  117.      * @return detection settings
  118.      * @since 13.0
  119.      */
  120.     public static EventDetectionSettings getDefaultEclipseDetectionSettings() {
  121.         return new EventDetectionSettings(ECLIPSE_MAX_CHECK, ECLIPSE_THRESHOLD, EventDetectionSettings.DEFAULT_MAX_ITER);
  122.     }

  124.     /** {@inheritDoc} */
  125.     @Override
  126.     public <T extends CalculusFieldElement<T>> Stream<FieldEventDetector<T>> getFieldEventDetectors(final Field<T> field) {
  127.         final T zero = field.getZero();
  128.         @SuppressWarnings("unchecked")
  129.         final FieldEventDetector<T>[] detectors = (FieldEventDetector<T>[]) Array.newInstance(FieldEventDetector.class,
  130.                                                                                               2 * occultingBodies.size());
  131.         for (int i = 0; i < occultingBodies.size(); ++i) {
  132.             final OccultationEngine occulting = occultingBodies.get(i);
  133.             detectors[2 * i]     = new FieldEclipseDetector<>(field, occulting).
  134.                                    withUmbra().
  135.                                    withMargin(zero.newInstance(-ANGULAR_MARGIN)).
  136.                                    withDetectionSettings(new FieldEventDetectionSettings<>(field, eclipseDetectionSettings)).
  137.                                    withHandler(new FieldResetDerivativesOnEvent<>());
  138.             detectors[2 * i + 1] = new FieldEclipseDetector<>(field, occulting).
  139.                                    withPenumbra().
  140.                                    withMargin(zero.newInstance(ANGULAR_MARGIN)).
  141.                                    withDetectionSettings(new FieldEventDetectionSettings<>(field, eclipseDetectionSettings)).
  142.                                    withHandler(new FieldResetDerivativesOnEvent<>());
  143.         }
  144.         return Stream.concat(Stream.of(detectors), RadiationForceModel.super.getFieldEventDetectors(field));
  145.     }

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

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

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

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

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

  172.         return angle;
  173.     }

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

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

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

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

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

  202.         return angle;
  203.     }

  205.     /**
  206.      * Add a new occulting body.
  207.      * <p>
  208.      * Central body is already considered, it shall not be added this way.
  209.      * </p>
  210.      * @param provider body PV provider
  211.      * @param radius body mean radius
  212.      * @see #addOccultingBody(OneAxisEllipsoid)
  213.      */
  214.     public void addOccultingBody(final ExtendedPositionProvider provider, final double radius) {

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

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

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

  233.         addOccultingBody(sphericalOccultingBody);

  235.     }

  237.     /**
  238.      * Add a new occulting body.
  239.      * <p>
  240.      * Central body is already considered, it shall not be added this way.
  241.      * </p>
  242.      * @param occulting occulting body to add
  243.      * @see #addOccultingBody(ExtendedPositionProvider, double)
  244.      * @since 12.0
  245.      */
  246.     public void addOccultingBody(final OneAxisEllipsoid occulting) {

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

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

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

  258.     }

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

  273. }
Created with JaCoCo 0.8.12.202403310830