AbstractRadiationForceModel.java

  1. /* Copyright 2002-2020 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.util.stream.Stream;

  21. import org.hipparchus.Field;
  22. import org.hipparchus.RealFieldElement;
  23. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  24. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  25. import org.hipparchus.ode.events.Action;
  26. import org.hipparchus.util.FastMath;
  27. import org.hipparchus.util.MathArrays;
  28. import org.orekit.errors.OrekitException;
  29. import org.orekit.errors.OrekitMessages;
  30. import org.orekit.forces.AbstractForceModel;
  31. import org.orekit.propagation.FieldSpacecraftState;
  32. import org.orekit.propagation.SpacecraftState;
  33. import org.orekit.propagation.events.AbstractDetector;
  34. import org.orekit.propagation.events.EventDetector;
  35. import org.orekit.propagation.events.FieldAbstractDetector;
  36. import org.orekit.propagation.events.FieldEventDetector;
  37. import org.orekit.propagation.events.handlers.EventHandler;
  38. import org.orekit.propagation.events.handlers.FieldEventHandler;
  39. import org.orekit.utils.Constants;
  40. import org.orekit.utils.ExtendedPVCoordinatesProvider;

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

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

  53.     /** Central body model. */
  54.     private final double equatorialRadius;

  56.     /** Sun model. */
  57.     private final ExtendedPVCoordinatesProvider sun;

  59.     /**
  60.      * Constructor.
  61.      * @param sun Sun model
  62.      * @param equatorialRadius spherical shape model (for umbra/penumbra computation)
  63.      */
  64.     protected AbstractRadiationForceModel(final ExtendedPVCoordinatesProvider sun, final double equatorialRadius) {
  65.         this.sun              = sun;
  66.         this.equatorialRadius = equatorialRadius;
  67.     }

  69.     /** {@inheritDoc} */
  70.     @Override
  71.     public boolean dependsOnPositionOnly() {
  72.         return false;
  73.     }

  75.     /** {@inheritDoc} */
  76.     @Override
  77.     public Stream<EventDetector> getEventsDetectors() {
  78.         return Stream.of(new UmbraDetector(), new PenumbraDetector());
  79.     }

  81.     /** {@inheritDoc} */
  82.     @Override
  83.     public <T extends RealFieldElement<T>> Stream<FieldEventDetector<T>> getFieldEventsDetectors(final Field<T> field) {
  84.         return Stream.of(new FieldUmbraDetector<>(field), new FieldPenumbraDetector<>(field));
  85.     }

  87.     /**
  88.      * Get the useful angles for eclipse computation.
  89.      * @param sunPosition Sun position in the selected frame
  90.      * @param position the satellite's position in the selected frame
  91.      * @return the 3 angles {(satCentral, satSun), Central body apparent radius, Sun apparent radius}
  92.      */
  93.     protected double[] getEclipseAngles(final Vector3D sunPosition, final Vector3D position) {
  94.         final double[] angle = new double[3];

  96.         final Vector3D satSunVector = sunPosition.subtract(position);

  98.         // Sat-Sun / Sat-CentralBody angle
  99.         angle[0] = Vector3D.angle(satSunVector, position.negate());

  101.         // Central body apparent radius
  102.         final double r = position.getNorm();
  103.         if (r <= equatorialRadius) {
  104.             throw new OrekitException(OrekitMessages.TRAJECTORY_INSIDE_BRILLOUIN_SPHERE, r);
  105.         }
  106.         angle[1] = FastMath.asin(equatorialRadius / r);

  108.         // Sun apparent radius
  109.         angle[2] = FastMath.asin(Constants.SUN_RADIUS / satSunVector.getNorm());

  111.         return angle;
  112.     }

  114.     /**
  115.      * Get the useful angles for eclipse computation.
  116.      * @param sunPosition Sun position in the selected frame
  117.      * @param position the satellite's position in the selected frame.
  118.      * @param <T> extends RealFieldElement
  119.      * @return the 3 angles {(satCentral, satSun), Central body apparent radius, Sun apparent radius}
  120.      */
  121.     protected <T extends RealFieldElement<T>> T[] getEclipseAngles(final FieldVector3D<T> sunPosition, final FieldVector3D<T> position) {
  122.         final T[] angle = MathArrays.buildArray(position.getX().getField(), 3);

  124.         final FieldVector3D<T> mP           = position.negate();
  125.         final FieldVector3D<T> satSunVector = mP.add(sunPosition);

  127.         // Sat-Sun / Sat-CentralBody angle
  128.         angle[0] = FieldVector3D.angle(satSunVector, mP);

  130.         // Central body apparent radius
  131.         final T r = position.getNorm();
  132.         if (r.getReal() <= equatorialRadius) {
  133.             throw new OrekitException(OrekitMessages.TRAJECTORY_INSIDE_BRILLOUIN_SPHERE, r);
  134.         }
  135.         angle[1] = r.reciprocal().multiply(equatorialRadius).asin();

  137.         // Sun apparent radius
  138.         angle[2] = satSunVector.getNorm().reciprocal().multiply(Constants.SUN_RADIUS).asin();

  140.         return angle;
  141.     }


  144.     /** This class defines the umbra entry/exit detector. */
  145.     private class UmbraDetector extends AbstractDetector<UmbraDetector> {

  147.         /** Build a new instance. */
  148.         UmbraDetector() {
  149.             super(60.0, 1.0e-3, DEFAULT_MAX_ITER, new EventHandler<UmbraDetector>() {

  151.                 /** {@inheritDoc} */
  152.                 public Action eventOccurred(final SpacecraftState s, final UmbraDetector detector,
  153.                                             final boolean increasing) {
  154.                     return Action.RESET_DERIVATIVES;
  155.                 }

  157.             });
  158.         }

  160.         /** Private constructor with full parameters.
  161.          * <p>
  162.          * This constructor is private as users are expected to use the builder
  163.          * API with the various {@code withXxx()} methods to set up the instance
  164.          * in a readable manner without using a huge amount of parameters.
  165.          * </p>
  166.          * @param maxCheck maximum checking interval (s)
  167.          * @param threshold convergence threshold (s)
  168.          * @param maxIter maximum number of iterations in the event time search
  169.          * @param handler event handler to call at event occurrences
  170.          * @since 6.1
  171.          */
  172.         private UmbraDetector(final double maxCheck, final double threshold,
  173.                               final int maxIter, final EventHandler<? super UmbraDetector> handler) {
  174.             super(maxCheck, threshold, maxIter, handler);
  175.         }

  177.         /** {@inheritDoc} */
  178.         @Override
  179.         protected UmbraDetector create(final double newMaxCheck, final double newThreshold,
  180.                                        final int newMaxIter, final EventHandler<? super UmbraDetector> newHandler) {
  181.             return new UmbraDetector(newMaxCheck, newThreshold, newMaxIter, newHandler);
  182.         }

  184.         /** The G-function is the difference between the Sun-Sat-Central-Body angle and
  185.          * the central body apparent radius.
  186.          * @param s the current state information : date, kinematics, attitude
  187.          * @return value of the g function
  188.          */
  189.         public double g(final SpacecraftState s) {
  190.             final double[] angle = getEclipseAngles(sun.getPVCoordinates(s.getDate(), s.getFrame()).getPosition(),
  191.                                                     s.getPVCoordinates().getPosition());
  192.             return angle[0] - angle[1] + angle[2] - ANGULAR_MARGIN;
  193.         }

  195.     }

  197.     /** This class defines the penumbra entry/exit detector. */
  198.     private class PenumbraDetector extends AbstractDetector<PenumbraDetector> {

  200.         /** Build a new instance. */
  201.         PenumbraDetector() {
  202.             super(60.0, 1.0e-3, DEFAULT_MAX_ITER, new EventHandler<PenumbraDetector>() {

  204.                 /** {@inheritDoc} */
  205.                 public Action eventOccurred(final SpacecraftState s, final PenumbraDetector detector,
  206.                                             final boolean increasing) {
  207.                     return Action.RESET_DERIVATIVES;
  208.                 }

  210.             });
  211.         }

  213.         /** Private constructor with full parameters.
  214.          * <p>
  215.          * This constructor is private as users are expected to use the builder
  216.          * API with the various {@code withXxx()} methods to set up the instance
  217.          * in a readable manner without using a huge amount of parameters.
  218.          * </p>
  219.          * @param maxCheck maximum checking interval (s)
  220.          * @param threshold convergence threshold (s)
  221.          * @param maxIter maximum number of iterations in the event time search
  222.          * @param handler event handler to call at event occurrences
  223.          * @since 6.1
  224.          */
  225.         private PenumbraDetector(final double maxCheck, final double threshold,
  226.                                  final int maxIter, final EventHandler<? super PenumbraDetector> handler) {
  227.             super(maxCheck, threshold, maxIter, handler);
  228.         }

  230.         /** {@inheritDoc} */
  231.         @Override
  232.         protected PenumbraDetector create(final double newMaxCheck, final double newThreshold,
  233.                                           final int newMaxIter, final EventHandler<? super PenumbraDetector> newHandler) {
  234.             return new PenumbraDetector(newMaxCheck, newThreshold, newMaxIter, newHandler);
  235.         }

  237.         /** The G-function is the difference between the Sun-Sat-Central-Body angle and
  238.          * the sum of the central body and Sun's apparent radius.
  239.          * @param s the current state information : date, kinematics, attitude
  240.          * @return value of the g function
  241.          */
  242.         public double g(final SpacecraftState s) {
  243.             final double[] angle = getEclipseAngles(sun.getPVCoordinates(s.getDate(), s.getFrame()).getPosition(),
  244.                                                     s.getPVCoordinates().getPosition());
  245.             return angle[0] - angle[1] - angle[2] + ANGULAR_MARGIN;
  246.         }

  248.     }

  250.     /** This class defines the umbra entry/exit detector. */
  251.     private class FieldUmbraDetector<T extends RealFieldElement<T>>
  252.         extends FieldAbstractDetector<FieldUmbraDetector<T>, T> {

  254.         /** Build a new instance.
  255.          * @param field field to which elements belong
  256.          */
  257.         FieldUmbraDetector(final Field<T> field) {
  258.             super(field.getZero().add(60.0), field.getZero().add(1.0e-3),
  259.                   DEFAULT_MAX_ITER, new FieldEventHandler<FieldUmbraDetector<T>, T>() {

  261.                       /** {@inheritDoc} */
  262.                       public Action eventOccurred(final FieldSpacecraftState<T> s,
  263.                                                   final FieldUmbraDetector<T> detector,
  264.                                                   final boolean increasing) {
  265.                           return Action.RESET_DERIVATIVES;
  266.                       }

  268.                   });
  269.         }

  271.         /** Private constructor with full parameters.
  272.          * <p>
  273.          * This constructor is private as users are expected to use the builder
  274.          * API with the various {@code withXxx()} methods to set up the instance
  275.          * in a readable manner without using a huge amount of parameters.
  276.          * </p>
  277.          * @param maxCheck maximum checking interval (s)
  278.          * @param threshold convergence threshold (s)
  279.          * @param maxIter maximum number of iterations in the event time search
  280.          * @param handler event handler to call at event occurrences
  281.          */
  282.         private FieldUmbraDetector(final T maxCheck, final T threshold,
  283.                                    final int maxIter,
  284.                                    final FieldEventHandler<? super FieldUmbraDetector<T>, T> handler) {
  285.             super(maxCheck, threshold, maxIter, handler);
  286.         }

  288.         /** {@inheritDoc} */
  289.         @Override
  290.         protected FieldUmbraDetector<T> create(final T newMaxCheck, final T newThreshold,
  291.                                                final int newMaxIter,
  292.                                                final FieldEventHandler<? super FieldUmbraDetector<T>, T> newHandler) {
  293.             return new FieldUmbraDetector<>(newMaxCheck, newThreshold, newMaxIter, newHandler);
  294.         }

  296.         /** The G-function is the difference between the Sun-Sat-Central-Body angle and
  297.          * the central body apparent radius.
  298.          * @param s the current state information : date, kinematics, attitude
  299.          * @return value of the g function
  300.          */
  301.         public T g(final FieldSpacecraftState<T> s) {
  302.             final T[] angle = getEclipseAngles(sun.getPVCoordinates(s.getDate(), s.getFrame()).getPosition(),
  303.                                                s.getPVCoordinates().getPosition());
  304.             return angle[0].subtract(angle[1]).add(angle[2]).subtract(ANGULAR_MARGIN);
  305.         }

  307.     }

  309.     /** This class defines the penumbra entry/exit detector. */
  310.     private class FieldPenumbraDetector<T extends RealFieldElement<T>>
  311.           extends FieldAbstractDetector<FieldPenumbraDetector<T>, T> {

  313.         /** Build a new instance.
  314.          * @param field field to which elements belong
  315.          */
  316.         FieldPenumbraDetector(final Field<T> field) {
  317.             super(field.getZero().add(60.0), field.getZero().add(1.0e-3),
  318.                   DEFAULT_MAX_ITER, new FieldEventHandler<FieldPenumbraDetector<T>, T>() {

  320.                       /** {@inheritDoc} */
  321.                       public Action eventOccurred(final FieldSpacecraftState<T> s,
  322.                                                   final FieldPenumbraDetector<T> detector,
  323.                                                   final boolean increasing) {
  324.                           return Action.RESET_DERIVATIVES;
  325.                       }

  327.                   });
  328.         }

  330.         /** Private constructor with full parameters.
  331.          * <p>
  332.          * This constructor is private as users are expected to use the builder
  333.          * API with the various {@code withXxx()} methods to set up the instance
  334.          * in a readable manner without using a huge amount of parameters.
  335.          * </p>
  336.          * @param maxCheck maximum checking interval (s)
  337.          * @param threshold convergence threshold (s)
  338.          * @param maxIter maximum number of iterations in the event time search
  339.          * @param handler event handler to call at event occurrences
  340.          */
  341.         private FieldPenumbraDetector(final T maxCheck, final T threshold,
  342.                                       final int maxIter,
  343.                                       final FieldEventHandler<? super FieldPenumbraDetector<T>, T> handler) {
  344.             super(maxCheck, threshold, maxIter, handler);
  345.         }

  347.         /** {@inheritDoc} */
  348.         @Override
  349.         protected FieldPenumbraDetector<T> create(final T newMaxCheck, final T newThreshold,
  350.                                                   final int newMaxIter,
  351.                                                   final FieldEventHandler<? super FieldPenumbraDetector<T>, T> newHandler) {
  352.             return new FieldPenumbraDetector<>(newMaxCheck, newThreshold, newMaxIter, newHandler);
  353.         }

  355.         /** The G-function is the difference between the Sun-Sat-Central-Body angle and
  356.          * the sum of the central body and Sun's apparent radius.
  357.          * @param s the current state information : date, kinematics, attitude
  358.          * @return value of the g function
  359.          */
  360.         public T g(final FieldSpacecraftState<T> s) {
  361.             final T[] angle = getEclipseAngles(sun.getPVCoordinates(s.getDate(), s.getFrame()).getPosition(),
  362.                                                s.getPVCoordinates().getPosition());
  363.             return angle[0].subtract(angle[1]).subtract(angle[2]).add(ANGULAR_MARGIN);
  364.         }

  366.     }

  368. }
Created with JaCoCo 0.8.6.202009150832