PositionAngleDetector.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.propagation.events;

  19. import java.util.function.Function;

  21. import org.hipparchus.analysis.UnivariateFunction;
  22. import org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver;
  23. import org.hipparchus.util.FastMath;
  24. import org.hipparchus.util.MathUtils;
  25. import org.orekit.errors.OrekitIllegalArgumentException;
  26. import org.orekit.errors.OrekitMessages;
  27. import org.orekit.orbits.CircularOrbit;
  28. import org.orekit.orbits.EquinoctialOrbit;
  29. import org.orekit.orbits.KeplerianOrbit;
  30. import org.orekit.orbits.Orbit;
  31. import org.orekit.orbits.OrbitType;
  32. import org.orekit.orbits.PositionAngleType;
  33. import org.orekit.propagation.SpacecraftState;
  34. import org.orekit.propagation.events.handlers.EventHandler;
  35. import org.orekit.propagation.events.handlers.StopOnEvent;
  36. import org.orekit.time.AbsoluteDate;
  37. import org.orekit.utils.TimeSpanMap;

  39. /** Detector for in-orbit position angle.
  40.  * <p>
  41.  * The detector is based on anomaly for {@link OrbitType#KEPLERIAN Keplerian}
  42.  * orbits, latitude argument for {@link OrbitType#CIRCULAR circular} orbits,
  43.  * or longitude argument for {@link OrbitType#EQUINOCTIAL equinoctial} orbits.
  44.  * It does not support {@link OrbitType#CARTESIAN Cartesian} orbits. The
  45.  * angles can be either {@link PositionAngleType#TRUE true}, {@link PositionAngleType#MEAN
  46.  * mean} or {@link PositionAngleType#ECCENTRIC eccentric} angles.
  47.  * </p>
  48.  * @author Luc Maisonobe
  49.  * @since 7.1
  50.  */
  51. public class PositionAngleDetector extends AbstractDetector<PositionAngleDetector> {

  53.     /** Orbit type defining the angle type. */
  54.     private final OrbitType orbitType;

  56.     /** Type of position angle. */
  57.     private final PositionAngleType positionAngleType;

  59.     /** Fixed angle to be crossed. */
  60.     private final double angle;

  62.     /** Position angle extraction function. */
  63.     private final Function<Orbit, Double> positionAngleExtractor;

  65.     /** Estimators for the offset angle, taking care of 2π wrapping and g function continuity. */
  66.     private TimeSpanMap<OffsetEstimator> offsetEstimators;

  68.     /** Build a new detector.
  69.      * <p>The new instance uses default values for maximal checking interval
  70.      * ({@link #DEFAULT_MAX_CHECK}) and convergence threshold ({@link
  71.      * #DEFAULT_THRESHOLD}).</p>
  72.      * @param orbitType orbit type defining the angle type
  73.      * @param positionAngleType type of position angle
  74.      * @param angle fixed angle to be crossed
  75.      * @exception OrekitIllegalArgumentException if orbit type is {@link OrbitType#CARTESIAN}
  76.      */
  77.     public PositionAngleDetector(final OrbitType orbitType, final PositionAngleType positionAngleType,
  78.                                  final double angle)
  79.         throws OrekitIllegalArgumentException {
  80.         this(DEFAULT_MAX_CHECK, DEFAULT_THRESHOLD, orbitType, positionAngleType, angle);
  81.     }

  83.     /** Build a detector.
  84.      * <p> This instance uses by default the {@link StopOnEvent} handler </p>
  85.      * @param maxCheck maximal checking interval (s)
  86.      * @param threshold convergence threshold (s)
  87.      * @param orbitType orbit type defining the angle type
  88.      * @param positionAngleType type of position angle
  89.      * @param angle fixed angle to be crossed
  90.      * @exception OrekitIllegalArgumentException if orbit type is {@link OrbitType#CARTESIAN}
  91.      */
  92.     public PositionAngleDetector(final double maxCheck, final double threshold,
  93.                                  final OrbitType orbitType, final PositionAngleType positionAngleType,
  94.                                  final double angle)
  95.         throws OrekitIllegalArgumentException {
  96.         this(new EventDetectionSettings(maxCheck, threshold, DEFAULT_MAX_ITER), new StopOnEvent(),
  97.              orbitType, positionAngleType, angle);
  98.     }

  100.     /** Protected constructor with full parameters.
  101.      * <p>
  102.      * This constructor is not public as users are expected to use the builder
  103.      * API with the various {@code withXxx()} methods to set up the instance
  104.      * in a readable manner without using a huge amount of parameters.
  105.      * </p>
  106.      * @param detectionSettings event detection settings
  107.      * @param handler event handler to call at event occurrences
  108.      * @param orbitType orbit type defining the angle type
  109.      * @param positionAngleType type of position angle
  110.      * @param angle fixed angle to be crossed
  111.      * @exception OrekitIllegalArgumentException if orbit type is {@link OrbitType#CARTESIAN}
  112.      * @since 13.0
  113.      */
  114.     protected PositionAngleDetector(final EventDetectionSettings detectionSettings, final EventHandler handler,
  115.                                     final OrbitType orbitType, final PositionAngleType positionAngleType,
  116.                                     final double angle)
  117.         throws OrekitIllegalArgumentException {

  119.         super(detectionSettings, handler);

  121.         this.orbitType        = orbitType;
  122.         this.positionAngleType = positionAngleType;
  123.         this.angle            = angle;
  124.         this.offsetEstimators = null;

  126.         switch (orbitType) {
  127.             case KEPLERIAN:
  128.                 positionAngleExtractor = o -> ((KeplerianOrbit) orbitType.convertType(o)).getAnomaly(positionAngleType);
  129.                 break;
  130.             case CIRCULAR:
  131.                 positionAngleExtractor = o -> ((CircularOrbit) orbitType.convertType(o)).getAlpha(positionAngleType);
  132.                 break;
  133.             case EQUINOCTIAL:
  134.                 positionAngleExtractor = o -> ((EquinoctialOrbit) orbitType.convertType(o)).getL(positionAngleType);
  135.                 break;
  136.             default:
  137.                 final String sep = ", ";
  138.                 throw new OrekitIllegalArgumentException(OrekitMessages.ORBIT_TYPE_NOT_ALLOWED,
  139.                                                          orbitType,
  140.                                                          OrbitType.KEPLERIAN   + sep +
  141.                                                          OrbitType.CIRCULAR    + sep +
  142.                                                          OrbitType.EQUINOCTIAL);
  143.         }

  145.     }

  147.     /** {@inheritDoc} */
  148.     @Override
  149.     protected PositionAngleDetector create(final EventDetectionSettings detectionSettings,
  150.                                            final EventHandler newHandler) {
  151.         return new PositionAngleDetector(detectionSettings, newHandler, orbitType, positionAngleType, angle);
  152.     }

  154.     /** Get the orbit type defining the angle type.
  155.      * @return orbit type defining the angle type
  156.      */
  157.     public OrbitType getOrbitType() {
  158.         return orbitType;
  159.     }

  161.     /** Get the type of position angle.
  162.      * @return type of position angle
  163.      */
  164.     public PositionAngleType getPositionAngleType() {
  165.         return positionAngleType;
  166.     }

  168.     /** Get the fixed angle to be crossed (radians).
  169.      * @return fixed angle to be crossed (radians)
  170.      */
  171.     public double getAngle() {
  172.         return angle;
  173.     }

  175.     /** {@inheritDoc} */
  176.     @Override
  177.     public void init(final SpacecraftState s0, final AbsoluteDate t) {
  178.         super.init(s0, t);
  179.         offsetEstimators = new TimeSpanMap<>(new OffsetEstimator(s0.getOrbit(), +1.0));
  180.     }

  182.     /** Compute the value of the detection function.
  183.      * <p>
  184.      * The value is the angle difference between the spacecraft and the fixed
  185.      * angle to be crossed, with some sign tweaks to ensure continuity.
  186.      * These tweaks imply the {@code increasing} flag in events detection becomes
  187.      * irrelevant here! As an example, the angle always increase in a Keplerian
  188.      * orbit, but this g function will increase and decrease so it
  189.      * will cross the zero value once per orbit, in increasing and decreasing
  190.      * directions on alternate orbits..
  191.      * </p>
  192.      * @param s the current state information: date, kinematics, attitude
  193.      * @return angle difference between the spacecraft and the fixed
  194.      * angle, with some sign tweaks to ensure continuity
  195.      */
  196.     public double g(final SpacecraftState s) {

  198.         final Orbit orbit = s.getOrbit();

  200.         // angle difference
  201.         OffsetEstimator estimator = offsetEstimators.get(s.getDate());
  202.         double          delta     = estimator.delta(orbit);

  204.         // we use a value greater than π for handover in order to avoid
  205.         // several switches to be estimated as the calling propagator
  206.         // and Orbit.shiftedBy have different accuracy. It is sufficient
  207.         // to have a handover roughly opposite to the detected position angle
  208.         while (FastMath.abs(delta) >= 3.5) {
  209.             // we are too far away from the current estimator, we need to set up a new one
  210.             // ensuring that we do have a crossing event in the current orbit
  211.             // and we ensure sign continuity with the current estimator

  213.             // find when the previous estimator becomes invalid
  214.             final AbsoluteDate handover = estimator.dateForOffset(FastMath.copySign(FastMath.PI, delta), orbit);

  216.             // perform handover to a new estimator at this date
  217.             estimator = new OffsetEstimator(orbit, delta);
  218.             delta     = estimator.delta(orbit);
  219.             if (isForward()) {
  220.                 offsetEstimators.addValidAfter(estimator, handover.getDate(), false);
  221.             } else {
  222.                 offsetEstimators.addValidBefore(estimator, handover.getDate(), false);
  223.             }

  225.         }

  227.         return delta;

  229.     }

  231.     /** Local class for estimating offset angle, handling 2π wrap-up and sign continuity. */
  232.     private class OffsetEstimator {

  234.         /** Target angle. */
  235.         private final double target;

  237.         /** Sign correction to offset. */
  238.         private final double sign;

  240.         /** Reference angle. */
  241.         private final double r0;

  243.         /** Slope of the linearized model. */
  244.         private final double r1;

  246.         /** Reference date. */
  247.         private final AbsoluteDate t0;

  249.         /** Simple constructor.
  250.          * @param orbit current orbit
  251.          * @param currentSign desired sign of the offset at current orbit time (magnitude is ignored)
  252.          */
  253.         OffsetEstimator(final Orbit orbit, final double currentSign) {
  254.             r0     = positionAngleExtractor.apply(orbit);
  255.             target = MathUtils.normalizeAngle(angle, r0);
  256.             sign   = FastMath.copySign(1.0, (r0 - target) * currentSign);
  257.             r1     = orbit.getKeplerianMeanMotion();
  258.             t0     = orbit.getDate();
  259.         }

  261.         /** Compute offset from reference angle.
  262.          * @param orbit current orbit
  263.          * @return offset between current angle and reference angle
  264.          */
  265.         public double delta(final Orbit orbit) {
  266.             final double rawAngle        = positionAngleExtractor.apply(orbit);
  267.             final double linearReference = r0 + r1 * orbit.getDate().durationFrom(t0);
  268.             final double linearizedAngle = MathUtils.normalizeAngle(rawAngle, linearReference);
  269.             return sign * (linearizedAngle - target);
  270.         }

  272.         /** Find date at which offset reaches specified value.
  273.          * <p>
  274.          * This computation is an approximation because it relies on
  275.          * {@link Orbit#shiftedBy(double)} only.
  276.          * </p>
  277.          * @param offset target value for offset angle
  278.          * @param orbit current orbit
  279.          * @return approximate date at which offset reached specified value
  280.          */
  281.         public AbsoluteDate dateForOffset(final double offset, final Orbit orbit) {

  283.             // bracket the search
  284.             final double period = orbit.getKeplerianPeriod();
  285.             final double delta0 = delta(orbit);
  286.             final double searchInf;
  287.             final double searchSup;
  288.             if ((delta0 - offset) * sign >= 0) {
  289.                 // the date is before current orbit
  290.                 searchInf = -period;
  291.                 searchSup = 0;
  292.             } else {
  293.                 // the date is after current orbit
  294.                 searchInf = 0;
  295.                 searchSup = +period;
  296.             }

  298.             // find the date as an offset from current orbit
  299.             final BracketingNthOrderBrentSolver solver = new BracketingNthOrderBrentSolver(getThreshold(), 5);
  300.             final UnivariateFunction            f      = dt -> delta(orbit.shiftedBy(dt)) - offset;
  301.             final double                        root   = solver.solve(getMaxIterationCount(), f, searchInf, searchSup);

  303.             return orbit.getDate().shiftedBy(root);

  305.         }

  307.     }

  309. }
Created with JaCoCo 0.8.12.202403310830