AbstractRadiationForceModel.java

  1. /* Copyright 2002-2023 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;

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

  23. import org.hipparchus.CalculusFieldElement;
  24. import org.hipparchus.Field;
  25. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  26. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  27. import org.hipparchus.ode.events.Action;
  28. import org.hipparchus.util.FastMath;
  29. import org.hipparchus.util.MathArrays;
  30. import org.orekit.bodies.OneAxisEllipsoid;
  31. import org.orekit.forces.ForceModel;
  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.FieldEclipseDetector;
  36. import org.orekit.propagation.events.FieldEventDetector;
  37. import org.orekit.utils.Constants;
  38. import org.orekit.utils.ExtendedPVCoordinatesProvider;
  39. import org.orekit.utils.ExtendedPVCoordinatesProviderAdapter;
  40. import org.orekit.utils.OccultationEngine;

  41. /**
  42.  * Base class for radiation force models.
  43.  * @see SolarRadiationPressure
  44.  * @see ECOM2
  45.  * @since 10.2
  46.  */
  47. public abstract class AbstractRadiationForceModel implements ForceModel {

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

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

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

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

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

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

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

  74.     /** {@inheritDoc} */
  75.     @Override
  76.     public boolean dependsOnPositionOnly() {
  77.         return false;
  78.     }

  79.     /** {@inheritDoc} */
  80.     @Override
  81.     public Stream<EventDetector> getEventDetectors() {
  82.         final EventDetector[] detectors = new EventDetector[2 * occultingBodies.size()];
  83.         for (int i = 0; i < occultingBodies.size(); ++i) {
  84.             final OccultationEngine occulting = occultingBodies.get(i);
  85.             detectors[2 * i]     = new EclipseDetector(occulting).
  86.                                    withUmbra().
  87.                                    withMargin(-ANGULAR_MARGIN).
  88.                                    withMaxCheck(ECLIPSE_MAX_CHECK).
  89.                                    withThreshold(ECLIPSE_THRESHOLD).
  90.                                    withHandler((state, detector, increasing) -> Action.RESET_DERIVATIVES);
  91.             detectors[2 * i + 1] = new EclipseDetector(occulting).
  92.                                    withPenumbra().
  93.                                    withMargin(ANGULAR_MARGIN).
  94.                                    withMaxCheck(ECLIPSE_MAX_CHECK).
  95.                                    withThreshold(ECLIPSE_THRESHOLD).
  96.                                    withHandler((state, detector, increasing) -> Action.RESET_DERIVATIVES);
  97.         }
  98.         // Fusion between Date detector for parameter driver span change and
  99.         // Detector for umbra / penumbra events
  100.         return Stream.concat(Stream.of(detectors), ForceModel.super.getEventDetectors());
  101.     }

  102.     /** {@inheritDoc} */
  103.     @Override
  104.     public <T extends CalculusFieldElement<T>> Stream<FieldEventDetector<T>> getFieldEventDetectors(final Field<T> field) {
  105.         final T zero = field.getZero();
  106.         @SuppressWarnings("unchecked")
  107.         final FieldEventDetector<T>[] detectors = (FieldEventDetector<T>[]) Array.newInstance(FieldEventDetector.class,
  108.                                                                                               2 * occultingBodies.size());
  109.         for (int i = 0; i < occultingBodies.size(); ++i) {
  110.             final OccultationEngine occulting = occultingBodies.get(i);
  111.             detectors[2 * i]     = new FieldEclipseDetector<>(field, occulting).
  112.                                    withUmbra().
  113.                                    withMargin(zero.newInstance(-ANGULAR_MARGIN)).
  114.                                    withMaxCheck(ECLIPSE_MAX_CHECK).
  115.                                    withThreshold(zero.newInstance(ECLIPSE_THRESHOLD)).
  116.                                    withHandler((state, detector, increasing) -> Action.RESET_DERIVATIVES);
  117.             detectors[2 * i + 1] = new FieldEclipseDetector<>(field, occulting).
  118.                                    withPenumbra().
  119.                                    withMargin(zero.newInstance(ANGULAR_MARGIN)).
  120.                                    withMaxCheck(ECLIPSE_MAX_CHECK).
  121.                                    withThreshold(zero.newInstance(ECLIPSE_THRESHOLD)).
  122.                                    withHandler((state, detector, increasing) -> Action.RESET_DERIVATIVES);
  123.         }
  124.         return Stream.concat(Stream.of(detectors), ForceModel.super.getFieldEventDetectors(field));
  125.     }

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

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

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

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

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

  146.         return angle;
  147.     }

  148.     /**
  149.      * Get the useful angles for eclipse computation.
  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.      * @param position the satellite's position in the selected frame
  155.      * @param <T> extends RealFieldElement
  156.      * @return the 3 angles {(satOcculting, satOcculted), Occulting body apparent radius, Occulted body apparent radius}
  157.      */
  158.     protected <T extends CalculusFieldElement<T>> T[] getGeneralEclipseAngles(final FieldVector3D<T> position,
  159.                                                                               final FieldVector3D<T> occultingPosition, final T occultingRadius,
  160.                                                                               final FieldVector3D<T> occultedPosition, final T occultedRadius) {
  161.         final T[] angle = MathArrays.buildArray(position.getX().getField(), 3);

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

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

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

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

  170.         return angle;
  171.     }

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

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

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

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

  196.         addOccultingBody(sphericalOccultingBody);

  197.     }

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

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

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

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

  215.     }

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

  228. }