AbstractRadiationForceModel.java

  1. /* Copyright 2002-2024 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.ode.events.Action;
  30. import org.hipparchus.util.FastMath;
  31. import org.hipparchus.util.MathArrays;
  32. import org.orekit.bodies.OneAxisEllipsoid;
  33. import org.orekit.frames.Frame;
  34. import org.orekit.propagation.events.EclipseDetector;
  35. import org.orekit.propagation.events.EventDetector;
  36. import org.orekit.propagation.events.FieldEclipseDetector;
  37. import org.orekit.propagation.events.FieldEventDetector;
  38. import org.orekit.utils.Constants;
  39. import org.orekit.utils.ExtendedPVCoordinatesProvider;
  40. import org.orekit.utils.ExtendedPVCoordinatesProviderAdapter;
  41. import org.orekit.utils.OccultationEngine;

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

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

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

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

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

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

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

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

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

  108.     /** {@inheritDoc} */
  109.     @Override
  110.     public <T extends CalculusFieldElement<T>> Stream<FieldEventDetector<T>> getFieldEventDetectors(final Field<T> field) {
  111.         final T zero = field.getZero();
  112.         @SuppressWarnings("unchecked")
  113.         final FieldEventDetector<T>[] detectors = (FieldEventDetector<T>[]) Array.newInstance(FieldEventDetector.class,
  114.                                                                                               2 * occultingBodies.size());
  115.         for (int i = 0; i < occultingBodies.size(); ++i) {
  116.             final OccultationEngine occulting = occultingBodies.get(i);
  117.             detectors[2 * i]     = new FieldEclipseDetector<>(field, occulting).
  118.                                    withUmbra().
  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.             detectors[2 * i + 1] = new FieldEclipseDetector<>(field, occulting).
  124.                                    withPenumbra().
  125.                                    withMargin(zero.newInstance(ANGULAR_MARGIN)).
  126.                                    withMaxCheck(ECLIPSE_MAX_CHECK).
  127.                                    withThreshold(zero.newInstance(ECLIPSE_THRESHOLD)).
  128.                                    withHandler((state, detector, increasing) -> Action.RESET_DERIVATIVES);
  129.         }
  130.         return Stream.concat(Stream.of(detectors), RadiationForceModel.super.getFieldEventDetectors(field));
  131.     }

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

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

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

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

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

  158.         return angle;
  159.     }

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

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

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

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

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

  188.         return angle;
  189.     }

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

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

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

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

  219.         addOccultingBody(sphericalOccultingBody);

  221.     }

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

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

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

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

  244.     }

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

  259. }
Created with JaCoCo 0.8.12.202403310830