FixedPointConverter.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.conversion.osc2mean;

  19. import org.hipparchus.CalculusFieldElement;
  20. import org.hipparchus.Field;
  21. import org.hipparchus.util.FastMath;
  22. import org.hipparchus.util.MathUtils;
  23. import org.orekit.errors.OrekitException;
  24. import org.orekit.errors.OrekitMessages;
  25. import org.orekit.orbits.EquinoctialOrbit;
  26. import org.orekit.orbits.FieldEquinoctialOrbit;
  27. import org.orekit.orbits.FieldOrbit;
  28. import org.orekit.orbits.Orbit;
  29. import org.orekit.orbits.PositionAngleType;
  30. import org.orekit.time.FieldAbsoluteDate;

  32. /**
  33.  * Class enabling conversion from osculating to mean orbit
  34.  * for a given theory using a fixed-point algorithm.
  35.  *
  36.  * @author Pascal Parraud
  37.  * @since 13.0
  38.  */
  39. public class FixedPointConverter implements OsculatingToMeanConverter {

  41.     /** Default convergence threshold. */
  42.     public static final double DEFAULT_THRESHOLD   = 1e-12;

  44.     /** Default maximum number of iterations. */
  45.     public static final int DEFAULT_MAX_ITERATIONS = 100;

  47.     /** Default damping ratio. */
  48.     public static final double DEFAULT_DAMPING     = 1.;

  50.     /** Mean theory used. */
  51.     private MeanTheory theory;

  53.     /** Convergence threshold. */
  54.     private double threshold;

  56.     /** Maximum number of iterations. */
  57.     private int maxIterations;

  59.     /** Damping ratio. */
  60.     private double damping;

  62.     /** Number of iterations performed. */
  63.     private int iterationsNb;

  65.     /**
  66.      * Default constructor.
  67.      * <p>
  68.      * The mean theory must be set before converting.
  69.      */
  70.     public FixedPointConverter() {
  71.         this(null, DEFAULT_THRESHOLD, DEFAULT_MAX_ITERATIONS, DEFAULT_DAMPING);
  72.     }

  74.     /**
  75.      * Constructor.
  76.      * @param theory mean theory to be used
  77.      */
  78.     public FixedPointConverter(final MeanTheory theory) {
  79.         this(theory, DEFAULT_THRESHOLD, DEFAULT_MAX_ITERATIONS, DEFAULT_DAMPING);
  80.     }

  82.     /**
  83.      * Constructor.
  84.      * <p>
  85.      * The mean theory must be set before converting.
  86.      *
  87.      * @param threshold tolerance for convergence
  88.      * @param maxIterations maximum number of iterations
  89.      * @param damping damping ratio
  90.      */
  91.     public FixedPointConverter(final double threshold,
  92.                                final int maxIterations,
  93.                                final double damping) {
  94.         this(null, threshold, maxIterations, damping);
  95.     }

  97.     /**
  98.      * Constructor.
  99.      * @param theory mean theory to be used
  100.      * @param threshold tolerance for convergence
  101.      * @param maxIterations maximum number of iterations
  102.      * @param damping damping ratio
  103.      */
  104.     public FixedPointConverter(final MeanTheory theory,
  105.                                final double threshold,
  106.                                final int maxIterations,
  107.                                final double damping) {
  108.         setMeanTheory(theory);
  109.         setThreshold(threshold);
  110.         setMaxIterations(maxIterations);
  111.         setDamping(damping);
  112.     }

  114.     /** {@inheritDoc} */
  115.     @Override
  116.     public MeanTheory getMeanTheory() {
  117.         return theory;
  118.     }

  120.     /** {@inheritDoc} */
  121.     @Override
  122.     public void setMeanTheory(final MeanTheory meanTheory) {
  123.         this.theory = meanTheory;
  124.     }

  126.     /**
  127.      * Gets convergence threshold.
  128.      * @return convergence threshold
  129.      */
  130.     public double getThreshold() {
  131.         return threshold;
  132.     }

  134.     /**
  135.      * Sets convergence threshold.
  136.      * @param threshold convergence threshold
  137.      */
  138.     public void setThreshold(final double threshold) {
  139.         this.threshold = threshold;
  140.     }

  142.     /**
  143.      * Gets maximum number of iterations.
  144.      * @return maximum number of iterations
  145.      */
  146.     public int getMaxIterations() {
  147.         return maxIterations;
  148.     }

  150.     /**
  151.      * Sets maximum number of iterations.
  152.      * @param maxIterations maximum number of iterations
  153.      */
  154.     public void setMaxIterations(final int maxIterations) {
  155.         this.maxIterations = maxIterations;
  156.     }

  158.     /**
  159.      * Gets damping ratio.
  160.      * @return damping ratio
  161.      */
  162.     public double getDamping() {
  163.         return damping;
  164.     }

  166.     /**
  167.      * Sets damping ratio.
  168.      * @param damping damping ratio
  169.      */
  170.     public void setDamping(final double damping) {
  171.         this.damping = damping;
  172.     }

  174.     /**
  175.      * Gets the number of iterations performed by the last conversion.
  176.      * @return number of iterations
  177.      */
  178.     public int getIterationsNb() {
  179.         return iterationsNb;
  180.     }

  182.     /** {@inheritDoc}
  183.      *  Uses a fixed-point algorithm.
  184.      */
  185.     @Override
  186.     public Orbit convertToMean(final Orbit osculating) {

  188.         // sanity check
  189.         if (osculating.getA() < theory.getReferenceRadius()) {
  190.             throw new OrekitException(OrekitMessages.TRAJECTORY_INSIDE_BRILLOUIN_SPHERE,
  191.                                       osculating.getA());
  192.         }

  194.         // Get equinoctial osculating parameters
  195.         final Orbit equinoctial = theory.preprocessing(osculating);
  196.         double sma = equinoctial.getA();
  197.         double ex  = equinoctial.getEquinoctialEx();
  198.         double ey  = equinoctial.getEquinoctialEy();
  199.         double hx  = equinoctial.getHx();
  200.         double hy  = equinoctial.getHy();
  201.         double lv  = equinoctial.getLv();

  203.         // Set threshold for each parameter
  204.         final double thresholdA  = threshold * FastMath.abs(sma);
  205.         final double thresholdE  = threshold * (1 + FastMath.hypot(ex, ey));
  206.         final double thresholdH  = threshold * (1 + FastMath.hypot(hx, hy));
  207.         final double thresholdLv = threshold * FastMath.PI;

  209.         // Rough initialization of the mean parameters
  210.         Orbit mean = theory.initialize(equinoctial);

  212.         int i = 0;
  213.         while (i++ < maxIterations) {

  215.             // Update osculating parameters from current mean parameters
  216.             final Orbit updated = theory.meanToOsculating(mean);

  218.             // Updated parameters residuals
  219.             final double deltaA  = equinoctial.getA() - updated.getA();
  220.             final double deltaEx = equinoctial.getEquinoctialEx() - updated.getEquinoctialEx();
  221.             final double deltaEy = equinoctial.getEquinoctialEy() - updated.getEquinoctialEy();
  222.             final double deltaHx = equinoctial.getHx() - updated.getHx();
  223.             final double deltaHy = equinoctial.getHy() - updated.getHy();
  224.             final double deltaLv = MathUtils.normalizeAngle(equinoctial.getLv() - updated.getLv(), 0.0);

  226.             // Check convergence
  227.             if (FastMath.abs(deltaA)  < thresholdA &&
  228.                 FastMath.abs(deltaEx) < thresholdE &&
  229.                 FastMath.abs(deltaEy) < thresholdE &&
  230.                 FastMath.abs(deltaHx) < thresholdH &&
  231.                 FastMath.abs(deltaHy) < thresholdH &&
  232.                 FastMath.abs(deltaLv) < thresholdLv) {
  233.                 // Records number of iterations performed
  234.                 iterationsNb = i;
  235.                 // Returns the mean orbit
  236.                 return theory.postprocessing(osculating, mean);
  237.             }

  239.             // Update mean parameters
  240.             sma += damping * deltaA;
  241.             ex  += damping * deltaEx;
  242.             ey  += damping * deltaEy;
  243.             hx  += damping * deltaHx;
  244.             hy  += damping * deltaHy;
  245.             lv  += damping * deltaLv;

  247.             // Update mean orbit
  248.             mean = new EquinoctialOrbit(sma, ex, ey, hx, hy, lv,
  249.                                         PositionAngleType.TRUE,
  250.                                         equinoctial.getFrame(),
  251.                                         equinoctial.getDate(),
  252.                                         equinoctial.getMu());
  253.         }
  254.         throw new OrekitException(OrekitMessages.UNABLE_TO_COMPUTE_MEAN_PARAMETERS, theory.getTheoryName(), i);
  255.     }

  257.     /** {@inheritDoc}
  258.      *  Uses a fixed-point algorithm.
  259.      */
  260.     @Override
  261.     public <T extends CalculusFieldElement<T>> FieldOrbit<T> convertToMean(final FieldOrbit<T> osculating) {

  263.         // Sanity check
  264.         if (osculating.getA().getReal() < theory.getReferenceRadius()) {
  265.             throw new OrekitException(OrekitMessages.TRAJECTORY_INSIDE_BRILLOUIN_SPHERE,
  266.                                            osculating.getA());
  267.         }

  269.         // Get field
  270.         final FieldAbsoluteDate<T> date = osculating.getDate();
  271.         final Field<T> field = date.getField();
  272.         final T zero = field.getZero();
  273.         final T pi   = zero.getPi();

  275.         // Get equinoctial parameters
  276.         final FieldOrbit<T> equinoctial = theory.preprocessing(osculating);
  277.         T sma = equinoctial.getA();
  278.         T ex  = equinoctial.getEquinoctialEx();
  279.         T ey  = equinoctial.getEquinoctialEy();
  280.         T hx  = equinoctial.getHx();
  281.         T hy  = equinoctial.getHy();
  282.         T lv  = equinoctial.getLv();

  284.         // Set threshold for each parameter
  285.         final T thresholdA  = sma.abs().multiply(threshold);
  286.         final T thresholdE  = FastMath.hypot(ex, ey).add(1).multiply(threshold);
  287.         final T thresholdH  = FastMath.hypot(hx, hy).add(1).multiply(threshold);
  288.         final T thresholdLv = pi.multiply(threshold);

  290.         // Rough initialization of the mean parameters
  291.         FieldOrbit<T> mean = theory.initialize(equinoctial);

  293.         int i = 0;
  294.         while (i++ < maxIterations) {

  296.             // recompute the osculating parameters from the current mean parameters
  297.             final FieldOrbit<T> updated = theory.meanToOsculating(mean);

  299.             // Updated parameters residuals
  300.             final T deltaA  = equinoctial.getA().subtract(updated.getA());
  301.             final T deltaEx = equinoctial.getEquinoctialEx().subtract(updated.getEquinoctialEx());
  302.             final T deltaEy = equinoctial.getEquinoctialEy().subtract(updated.getEquinoctialEy());
  303.             final T deltaHx = equinoctial.getHx().subtract(updated.getHx());
  304.             final T deltaHy = equinoctial.getHy().subtract(updated.getHy());
  305.             final T deltaLv = MathUtils.normalizeAngle(equinoctial.getLv().subtract(updated.getLv()), zero);

  307.             // Check convergence
  308.             if (FastMath.abs(deltaA.getReal())  < thresholdA.getReal() &&
  309.                 FastMath.abs(deltaEx.getReal()) < thresholdE.getReal() &&
  310.                 FastMath.abs(deltaEy.getReal()) < thresholdE.getReal() &&
  311.                 FastMath.abs(deltaHx.getReal()) < thresholdH.getReal() &&
  312.                 FastMath.abs(deltaHy.getReal()) < thresholdH.getReal() &&
  313.                 FastMath.abs(deltaLv.getReal()) < thresholdLv.getReal()) {
  314.                 // Records number of iterations performed
  315.                 iterationsNb = i;
  316.                 // Returns the mean orbit
  317.                 return theory.postprocessing(osculating, mean);
  318.             }

  320.             // Update mean parameters
  321.             sma = sma.add(deltaA.multiply(damping));
  322.             ex  = ex.add(deltaEx.multiply(damping));
  323.             ey  = ey.add(deltaEy.multiply(damping));
  324.             hx  = hx.add(deltaHx.multiply(damping));
  325.             hy  = hy.add(deltaHy.multiply(damping));
  326.             lv  = lv.add(deltaLv.multiply(damping));

  328.             // Update mean orbit
  329.             mean = new FieldEquinoctialOrbit<>(sma, ex, ey, hx, hy, lv,
  330.                                                PositionAngleType.TRUE,
  331.                                                equinoctial.getFrame(),
  332.                                                equinoctial.getDate(),
  333.                                                equinoctial.getMu());
  334.         }
  335.         throw new OrekitException(OrekitMessages.UNABLE_TO_COMPUTE_MEAN_PARAMETERS, theory.getTheoryName(), i);
  336.     }
  337. }
Created with JaCoCo 0.8.12.202403310830