FixedPointTleGenerationAlgorithm.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.propagation.analytical.tle.generation;

  19. import org.hipparchus.CalculusFieldElement;
  20. import org.hipparchus.Field;
  21. import org.hipparchus.geometry.euclidean.threed.Rotation;
  22. import org.hipparchus.util.FastMath;
  23. import org.hipparchus.util.MathUtils;
  24. import org.orekit.annotation.DefaultDataContext;
  25. import org.orekit.attitudes.FrameAlignedProvider;
  26. import org.orekit.data.DataContext;
  27. import org.orekit.errors.OrekitException;
  28. import org.orekit.errors.OrekitMessages;
  29. import org.orekit.frames.Frame;
  30. import org.orekit.orbits.EquinoctialOrbit;
  31. import org.orekit.orbits.FieldEquinoctialOrbit;
  32. import org.orekit.orbits.FieldKeplerianOrbit;
  33. import org.orekit.orbits.FieldOrbit;
  34. import org.orekit.orbits.KeplerianOrbit;
  35. import org.orekit.orbits.Orbit;
  36. import org.orekit.orbits.OrbitType;
  37. import org.orekit.orbits.PositionAngleType;
  38. import org.orekit.propagation.FieldSpacecraftState;
  39. import org.orekit.propagation.SpacecraftState;
  40. import org.orekit.propagation.analytical.tle.FieldTLE;
  41. import org.orekit.propagation.analytical.tle.FieldTLEPropagator;
  42. import org.orekit.propagation.analytical.tle.TLE;
  43. import org.orekit.propagation.analytical.tle.TLEConstants;
  44. import org.orekit.propagation.analytical.tle.TLEPropagator;
  45. import org.orekit.time.AbsoluteDate;
  46. import org.orekit.time.TimeScale;
  47. import org.orekit.utils.ParameterDriver;

  49. /**
  50.  * Fixed Point method to reverse SGP4 and SDP4 propagation algorithm
  51.  * and generate a usable TLE from a spacecraft state.
  52.  * <p>
  53.  * Using this algorithm, the B* value is not computed. In other words,
  54.  * the B* value from the template TLE is set to the generated one.
  55.  * </p>
  56.  * @author Thomas Paulet
  57.  * @author Bryan Cazabonne
  58.  * @since 12.0
  59.  */
  60. public class FixedPointTleGenerationAlgorithm implements TleGenerationAlgorithm {

  62.     /** Default value for epsilon. */
  63.     public static final double EPSILON_DEFAULT = 1.0e-10;

  65.     /** Default value for maxIterations. */
  66.     public static final int MAX_ITERATIONS_DEFAULT = 100;

  68.     /** Default value for scale. */
  69.     public static final double SCALE_DEFAULT = 1.0;

  71.     /** Used to compute threshold for convergence check. */
  72.     private final double epsilon;

  74.     /** Maximum number of iterations for convergence. */
  75.     private final int maxIterations;

  77.     /** Scale factor of the Fixed Point algorithm. */
  78.     private final double scale;

  80.     /** UTC scale. */
  81.     private final TimeScale utc;

  83.     /** TEME frame. */
  84.     private final Frame teme;

  86.     /**
  87.      * Default constructor.
  88.      * <p>
  89.      * Uses the {@link DataContext#getDefault() default data context}
  90.      * as well as {@link #EPSILON_DEFAULT}, {@link #MAX_ITERATIONS_DEFAULT},
  91.      * {@link #SCALE_DEFAULT} for method convergence.
  92.      * </p>
  93.      */
  94.     @DefaultDataContext
  95.     public FixedPointTleGenerationAlgorithm() {
  96.         this(EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT, SCALE_DEFAULT);
  97.     }

  99.     /**
  100.      * Constructor.
  101.      * <p>
  102.      * Uses the {@link DataContext#getDefault() default data context}.
  103.      * </p>
  104.      * @param epsilon used to compute threshold for convergence check
  105.      * @param maxIterations maximum number of iterations for convergence
  106.      * @param scale scale factor of the Fixed Point algorithm
  107.      */
  108.     @DefaultDataContext
  109.     public FixedPointTleGenerationAlgorithm(final double epsilon, final int maxIterations,
  110.                                             final double scale) {
  111.         this(epsilon, maxIterations, scale,
  112.              DataContext.getDefault().getTimeScales().getUTC(),
  113.              DataContext.getDefault().getFrames().getTEME());
  114.     }

  116.     /**
  117.      * Constructor.
  118.      * @param epsilon used to compute threshold for convergence check
  119.      * @param maxIterations maximum number of iterations for convergence
  120.      * @param scale scale factor of the Fixed Point algorithm
  121.      * @param utc UTC time scale
  122.      * @param teme TEME frame
  123.      */
  124.     public FixedPointTleGenerationAlgorithm(final double epsilon, final int maxIterations,
  125.                                             final double scale, final TimeScale utc,
  126.                                             final Frame teme) {
  127.         this.epsilon       = epsilon;
  128.         this.maxIterations = maxIterations;
  129.         this.scale         = scale;
  130.         this.utc           = utc;
  131.         this.teme          = teme;
  132.     }

  134.     /** {@inheritDoc} */
  135.     @Override
  136.     public TLE generate(final SpacecraftState state, final TLE templateTLE) {

  138.         // Generation epoch
  139.         final AbsoluteDate epoch = state.getDate();

  141.         // gets equinoctial parameters in TEME frame and with TLE gravity parameter from state
  142.         final EquinoctialOrbit equinoctialOrbit = convert(state.getOrbit());
  143.         double sma = equinoctialOrbit.getA();
  144.         double ex  = equinoctialOrbit.getEquinoctialEx();
  145.         double ey  = equinoctialOrbit.getEquinoctialEy();
  146.         double hx  = equinoctialOrbit.getHx();
  147.         double hy  = equinoctialOrbit.getHy();
  148.         double lv  = equinoctialOrbit.getLv();

  150.         // rough initialization of the TLE
  151.         final KeplerianOrbit keplerianOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(equinoctialOrbit);
  152.         TLE current = TleGenerationUtil.newTLE(keplerianOrbit, templateTLE, templateTLE.getBStar(epoch), utc);

  154.         // threshold for each parameter
  155.         final double thrA = epsilon * (1 + sma);
  156.         final double thrE = epsilon * (1 + FastMath.hypot(ex, ey));
  157.         final double thrH = epsilon * (1 + FastMath.hypot(hx, hy));
  158.         final double thrV = epsilon * FastMath.PI;

  160.         int k = 0;
  161.         while (k++ < maxIterations) {

  163.             // recompute the state from the current TLE
  164.             final TLEPropagator propagator = TLEPropagator.selectExtrapolator(current,
  165.                                                                               new FrameAlignedProvider(Rotation.IDENTITY, teme),
  166.                                                                               state.getMass(), teme);
  167.             final Orbit recoveredOrbit = propagator.getInitialState().getOrbit();
  168.             final EquinoctialOrbit recoveredEquiOrbit = (EquinoctialOrbit) OrbitType.EQUINOCTIAL.convertType(recoveredOrbit);

  170.             // adapted parameters residuals
  171.             final double deltaSma = equinoctialOrbit.getA() - recoveredEquiOrbit.getA();
  172.             final double deltaEx  = equinoctialOrbit.getEquinoctialEx() - recoveredEquiOrbit.getEquinoctialEx();
  173.             final double deltaEy  = equinoctialOrbit.getEquinoctialEy() - recoveredEquiOrbit.getEquinoctialEy();
  174.             final double deltaHx  = equinoctialOrbit.getHx() - recoveredEquiOrbit.getHx();
  175.             final double deltaHy  = equinoctialOrbit.getHy() - recoveredEquiOrbit.getHy();
  176.             final double deltaLv  = MathUtils.normalizeAngle(equinoctialOrbit.getLv() - recoveredEquiOrbit.getLv(), 0.0);

  178.             // check convergence
  179.             if (FastMath.abs(deltaSma) < thrA &&
  180.                 FastMath.abs(deltaEx)  < thrE &&
  181.                 FastMath.abs(deltaEy)  < thrE &&
  182.                 FastMath.abs(deltaHx)  < thrH &&
  183.                 FastMath.abs(deltaHy)  < thrH &&
  184.                 FastMath.abs(deltaLv)  < thrV) {

  186.                 // verify if parameters are estimated
  187.                 for (final ParameterDriver templateDrivers : templateTLE.getParametersDrivers()) {
  188.                     if (templateDrivers.isSelected()) {
  189.                         // set to selected for the new TLE
  190.                         current.getParameterDriver(templateDrivers.getName()).setSelected(true);
  191.                     }
  192.                 }

  194.                 // return
  195.                 return current;
  196.             }

  198.             // update state
  199.             sma += scale * deltaSma;
  200.             ex  += scale * deltaEx;
  201.             ey  += scale * deltaEy;
  202.             hx  += scale * deltaHx;
  203.             hy  += scale * deltaHy;
  204.             lv  += scale * deltaLv;
  205.             final EquinoctialOrbit newEquinoctialOrbit =
  206.                                     new EquinoctialOrbit(sma, ex, ey, hx, hy, lv, PositionAngleType.TRUE,
  207.                                                          equinoctialOrbit.getFrame(), equinoctialOrbit.getDate(), equinoctialOrbit.getMu());
  208.             final KeplerianOrbit newKeplerianOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(newEquinoctialOrbit);

  210.             // update TLE
  211.             current = TleGenerationUtil.newTLE(newKeplerianOrbit, templateTLE, templateTLE.getBStar(epoch), utc);

  213.         }

  215.         // unable to generate a TLE
  216.         throw new OrekitException(OrekitMessages.UNABLE_TO_COMPUTE_TLE, k);

  218.     }

  220.     /** {@inheritDoc} */
  221.     @Override
  222.     public <T extends CalculusFieldElement<T>> FieldTLE<T> generate(final FieldSpacecraftState<T> state,
  223.                                                                     final FieldTLE<T> templateTLE) {

  225.         // gets equinoctial parameters in TEME frame and with TLE gravity parameter from state
  226.         final FieldEquinoctialOrbit<T> equinoctialOrbit = convert(state.getOrbit());
  227.         T sma = equinoctialOrbit.getA();
  228.         T ex  = equinoctialOrbit.getEquinoctialEx();
  229.         T ey  = equinoctialOrbit.getEquinoctialEy();
  230.         T hx  = equinoctialOrbit.getHx();
  231.         T hy  = equinoctialOrbit.getHy();
  232.         T lv  = equinoctialOrbit.getLv();

  234.         // rough initialization of the TLE
  235.         final T bStar = state.getA().getField().getZero().newInstance(templateTLE.getBStar());
  236.         final FieldKeplerianOrbit<T> keplerianOrbit = (FieldKeplerianOrbit<T>) OrbitType.KEPLERIAN.convertType(equinoctialOrbit);
  237.         FieldTLE<T> current = TleGenerationUtil.newTLE(keplerianOrbit, templateTLE, bStar, utc);

  239.         // field
  240.         final Field<T> field = state.getDate().getField();

  242.         // threshold for each parameter
  243.         final T thrA = sma.add(1).multiply(epsilon);
  244.         final T thrE = FastMath.hypot(ex, ey).add(1).multiply(epsilon);
  245.         final T thrH = FastMath.hypot(hx, hy).add(1).multiply(epsilon);
  246.         final T thrV = sma.getPi().multiply(epsilon);

  248.         int k = 0;
  249.         while (k++ < maxIterations) {

  251.             // recompute the state from the current TLE
  252.             final FieldTLEPropagator<T> propagator = FieldTLEPropagator.selectExtrapolator(current, new FrameAlignedProvider(Rotation.IDENTITY, teme),
  253.                                                                                            state.getMass(), teme, templateTLE.getParameters(field));
  254.             final FieldOrbit<T> recoveredOrbit = propagator.getInitialState().getOrbit();
  255.             final FieldEquinoctialOrbit<T> recoveredEquinoctialOrbit = (FieldEquinoctialOrbit<T>) OrbitType.EQUINOCTIAL.convertType(recoveredOrbit);

  257.             // adapted parameters residuals
  258.             final T deltaSma = equinoctialOrbit.getA().subtract(recoveredEquinoctialOrbit.getA());
  259.             final T deltaEx  = equinoctialOrbit.getEquinoctialEx().subtract(recoveredEquinoctialOrbit.getEquinoctialEx());
  260.             final T deltaEy  = equinoctialOrbit.getEquinoctialEy().subtract(recoveredEquinoctialOrbit.getEquinoctialEy());
  261.             final T deltaHx  = equinoctialOrbit.getHx().subtract(recoveredEquinoctialOrbit.getHx());
  262.             final T deltaHy  = equinoctialOrbit.getHy().subtract(recoveredEquinoctialOrbit.getHy());
  263.             final T deltaLv  = MathUtils.normalizeAngle(equinoctialOrbit.getLv().subtract(recoveredEquinoctialOrbit.getLv()), field.getZero());

  265.             // check convergence
  266.             if (FastMath.abs(deltaSma.getReal()) < thrA.getReal() &&
  267.                 FastMath.abs(deltaEx.getReal())  < thrE.getReal() &&
  268.                 FastMath.abs(deltaEy.getReal())  < thrE.getReal() &&
  269.                 FastMath.abs(deltaHx.getReal())  < thrH.getReal() &&
  270.                 FastMath.abs(deltaHy.getReal())  < thrH.getReal() &&
  271.                 FastMath.abs(deltaLv.getReal())  < thrV.getReal()) {

  273.                 // return
  274.                 return current;

  276.             }

  278.             // update state
  279.             sma = sma.add(deltaSma.multiply(scale));
  280.             ex  = ex.add(deltaEx.multiply(scale));
  281.             ey  = ey.add(deltaEy.multiply(scale));
  282.             hx  = hx.add(deltaHx.multiply(scale));
  283.             hy  = hy.add(deltaHy.multiply(scale));
  284.             lv  = lv.add(deltaLv.multiply(scale));
  285.             final FieldEquinoctialOrbit<T> newEquinoctialOrbit =
  286.                                     new FieldEquinoctialOrbit<>(sma, ex, ey, hx, hy, lv, PositionAngleType.TRUE,
  287.                                     equinoctialOrbit.getFrame(), equinoctialOrbit.getDate(), equinoctialOrbit.getMu());
  288.             final FieldKeplerianOrbit<T> newKeplerianOrbit = (FieldKeplerianOrbit<T>) OrbitType.KEPLERIAN.convertType(newEquinoctialOrbit);

  290.             // update TLE
  291.             current = TleGenerationUtil.newTLE(newKeplerianOrbit, templateTLE, bStar, utc);

  293.         }

  295.         throw new OrekitException(OrekitMessages.UNABLE_TO_COMPUTE_TLE, k);

  297.     }

  299.     /**
  300.      * Converts an orbit into an equinoctial orbit expressed in TEME frame with the TLE gravity parameter.
  301.      * @param orbitIn the orbit to convert
  302.      * @return the converted orbit, i.e. equinoctial in TEME frame
  303.      */
  304.     private EquinoctialOrbit convert(final Orbit orbitIn) {
  305.         return new EquinoctialOrbit(orbitIn.getPVCoordinates(teme), teme, TLEConstants.MU);
  306.     }

  308.     /**
  309.      * Converts an orbit into an equinoctial orbit expressed in TEME frame with the TLE gravity parameter.
  310.      * @param orbitIn the orbit to convert
  311.      * @param <T> type of the element
  312.      * @return the converted orbit, i.e. equinoctial in TEME frame
  313.      */
  314.     private <T extends CalculusFieldElement<T>> FieldEquinoctialOrbit<T> convert(final FieldOrbit<T> orbitIn) {
  315.         return new FieldEquinoctialOrbit<T>(orbitIn.getPVCoordinates(teme), teme, orbitIn.getMu().newInstance(TLEConstants.MU));
  316.     }

  318. }
Created with JaCoCo 0.8.12.202403310830