BrouwerLyddanePropagatorBuilder.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;

  19. import java.util.Collections;
  20. import org.hipparchus.util.FastMath;
  21. import org.orekit.attitudes.AttitudeProvider;
  22. import org.orekit.attitudes.FrameAlignedProvider;
  23. import org.orekit.forces.gravity.potential.GravityFieldFactory;
  24. import org.orekit.forces.gravity.potential.TideSystem;
  25. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
  26. import org.orekit.orbits.Orbit;
  27. import org.orekit.orbits.OrbitType;
  28. import org.orekit.orbits.PositionAngleType;
  29. import org.orekit.propagation.Propagator;
  30. import org.orekit.propagation.analytical.BrouwerLyddanePropagator;
  31. import org.orekit.propagation.analytical.tle.TLE;
  32. import org.orekit.utils.ParameterDriver;

  34. /** Builder for Brouwer-Lyddane propagator.
  35.  * <p>
  36.  * By default, Brouwer-Lyddane model considers only the perturbations due to zonal harmonics.
  37.  * However, for low Earth orbits, the magnitude of the perturbative acceleration due to
  38.  * atmospheric drag can be significant. Warren Phipps' 1992 thesis considered the atmospheric
  39.  * drag by time derivatives of the <i>mean</i> mean anomaly using the catch-all coefficient M2.
  40.  *
  41.  * Usually, M2 is adjusted during an orbit determination process and it represents the
  42.  * combination of all unmodeled secular along-track effects (i.e. not just the atmospheric drag).
  43.  * The behavior of M2 is closed to the {@link TLE#getBStar()} parameter for the TLE.
  44.  *
  45.  * If the value of M2 is equal to {@link BrouwerLyddanePropagator#M2 0.0}, the along-track
  46.  * secular effects are not considered in the dynamical model. Typical values for M2 are not known.
  47.  * It depends on the orbit type. However, the value of M2 must be very small (e.g. between 1.0e-14 and 1.0e-15).
  48.  * The unit of M2 is rad/s².
  49.  * <p>
  50.  * To estimate the M2 parameter, it is necessary to call the {@link #getPropagationParametersDrivers()} method
  51.  * as follow:
  52.  * <pre>
  53.  *  for (ParameterDriver driver : builder.getPropagationParametersDrivers().getDrivers()) {
  54.  *     if (BrouwerLyddanePropagator.M2_NAME.equals(driver.getName())) {
  55.  *        driver.setSelected(true);
  56.  *     }
  57.  *  }
  58.  * </pre>
  59.  * @author Melina Vanel
  60.  * @author Bryan Cazabonne
  61.  * @since 11.1
  62.  */
  63. public class BrouwerLyddanePropagatorBuilder extends AbstractAnalyticalPropagatorBuilder<BrouwerLyddanePropagator> {

  65.     /** Parameters scaling factor.
  66.      * <p>
  67.      * We use a power of 2 to avoid numeric noise introduction
  68.      * in the multiplications/divisions sequences.
  69.      * </p>
  70.      */
  71.     private static final double SCALE = FastMath.scalb(1.0, -32);

  73.     /** Provider for un-normalized coefficients. */
  74.     private final UnnormalizedSphericalHarmonicsProvider provider;

  76.     /** Build a new instance.
  77.      * <p>
  78.      * The template orbit is used as a model to {@link
  79.      * #createInitialOrbit() create initial orbit}. It defines the
  80.      * inertial frame, the central attraction coefficient, the orbit type, and is also
  81.      * used together with the {@code positionScale} to convert from the {@link
  82.      * org.orekit.utils.ParameterDriver#setNormalizedValue(double) normalized} parameters
  83.      * used by the callers of this builder to the real orbital parameters.
  84.      * The default attitude provider is aligned with the orbit's inertial frame.
  85.      * </p>
  86.      *
  87.      * @param templateOrbit reference orbit from which real orbits will be built
  88.      * (note that the mu from this orbit will be overridden with the mu from the
  89.      * {@code provider})
  90.      * @param provider for un-normalized zonal coefficients
  91.      * @param positionAngleType position angle type to use
  92.      * @param positionScale scaling factor used for orbital parameters normalization
  93.      * (typically set to the expected standard deviation of the position)
  94.      * @param M2 value of empirical drag coefficient in rad/s².
  95.      *        If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
  96.      * @see #BrouwerLyddanePropagatorBuilder(Orbit,
  97.      * UnnormalizedSphericalHarmonicsProvider, PositionAngleType, double, AttitudeProvider, double)
  98.      */
  99.     public BrouwerLyddanePropagatorBuilder(final Orbit templateOrbit,
  100.                                            final UnnormalizedSphericalHarmonicsProvider provider,
  101.                                            final PositionAngleType positionAngleType,
  102.                                            final double positionScale,
  103.                                            final double M2) {
  104.         this(templateOrbit, provider, positionAngleType, positionScale,
  105.              FrameAlignedProvider.of(templateOrbit.getFrame()), M2);
  106.     }

  108.     /** Build a new instance.
  109.      * <p>
  110.      * The template orbit is used as a model to {@link
  111.      * #createInitialOrbit() create initial orbit}. It defines the
  112.      * inertial frame, the central attraction coefficient, the orbit type, and is also
  113.      * used together with the {@code positionScale} to convert from the {@link
  114.      * org.orekit.utils.ParameterDriver#setNormalizedValue(double) normalized} parameters used by the
  115.      * callers of this builder to the real orbital parameters.
  116.      * </p>
  117.      *
  118.      * @param templateOrbit reference orbit from which real orbits will be built
  119.      * (note that the mu from this orbit will be overridden with the mu from the
  120.      * {@code provider})
  121.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  122.      * @param mu central attraction coefficient (m³/s²)
  123.      * @param tideSystem tide system
  124.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  125.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  126.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  127.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  128.      * @param orbitType orbit type to use
  129.      * @param positionAngleType position angle type to use
  130.      * @param positionScale scaling factor used for orbital parameters normalization
  131.      * (typically set to the expected standard deviation of the position)
  132.      * @param M2 value of empirical drag coefficient in rad/s².
  133.      *        If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
  134.      * @see #BrouwerLyddanePropagatorBuilder(Orbit,
  135.      * UnnormalizedSphericalHarmonicsProvider, PositionAngleType, double, AttitudeProvider, double)
  136.      */
  137.     public BrouwerLyddanePropagatorBuilder(final Orbit templateOrbit,
  138.                                            final double referenceRadius,
  139.                                            final double mu,
  140.                                            final TideSystem tideSystem,
  141.                                            final double c20,
  142.                                            final double c30,
  143.                                            final double c40,
  144.                                            final double c50,
  145.                                            final OrbitType orbitType,
  146.                                            final PositionAngleType positionAngleType,
  147.                                            final double positionScale,
  148.                                            final double M2) {
  149.         this(templateOrbit,
  150.              GravityFieldFactory.getUnnormalizedProvider(referenceRadius, mu, tideSystem,
  151.                                                          new double[][] {
  152.                                                              {
  153.                                                                  0
  154.                                                              }, {
  155.                                                                  0
  156.                                                              }, {
  157.                                                                  c20
  158.                                                              }, {
  159.                                                                  c30
  160.                                                              }, {
  161.                                                                  c40
  162.                                                              }, {
  163.                                                                  c50
  164.                                                              }
  165.                                                          }, new double[][] {
  166.                                                              {
  167.                                                                  0
  168.                                                              }, {
  169.                                                                  0
  170.                                                              }, {
  171.                                                                  0
  172.                                                              }, {
  173.                                                                  0
  174.                                                              }, {
  175.                                                                  0
  176.                                                              }, {
  177.                                                                  0
  178.                                                              }
  179.                                                          }),
  180.                 positionAngleType, positionScale, M2);
  181.     }

  183.     /** Build a new instance.
  184.      * <p>
  185.      * The template orbit is used as a model to {@link
  186.      * #createInitialOrbit() create initial orbit}. It defines the
  187.      * inertial frame, the central attraction coefficient, the orbit type, and is also
  188.      * used together with the {@code positionScale} to convert from the {@link
  189.      * org.orekit.utils.ParameterDriver#setNormalizedValue(double) normalized} parameters used by the
  190.      * callers of this builder to the real orbital parameters.
  191.      * </p>
  192.      * @param templateOrbit reference orbit from which real orbits will be built
  193.      * (note that the mu from this orbit will be overridden with the mu from the
  194.      * {@code provider})
  195.      * @param provider for un-normalized zonal coefficients
  196.      * @param positionAngleType position angle type to use
  197.      * @param positionScale scaling factor used for orbital parameters normalization
  198.      * (typically set to the expected standard deviation of the position)
  199.      * @param attitudeProvider attitude law to use
  200.      * @param M2 value of empirical drag coefficient in rad/s².
  201.      *        If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
  202.      */
  203.     public BrouwerLyddanePropagatorBuilder(final Orbit templateOrbit,
  204.                                            final UnnormalizedSphericalHarmonicsProvider provider,
  205.                                            final PositionAngleType positionAngleType,
  206.                                            final double positionScale,
  207.                                            final AttitudeProvider attitudeProvider,
  208.                                            final double M2) {
  209.         super(overrideMu(templateOrbit, provider, positionAngleType), positionAngleType, positionScale, true, attitudeProvider,
  210.             Propagator.DEFAULT_MASS);
  211.         this.provider = provider;
  212.         // initialize M2 driver
  213.         final ParameterDriver M2Driver = new ParameterDriver(BrouwerLyddanePropagator.M2_NAME, M2, SCALE,
  214.                                                              Double.NEGATIVE_INFINITY,
  215.                                                              Double.POSITIVE_INFINITY);
  216.         addSupportedParameters(Collections.singletonList(M2Driver));
  217.     }

  219.     /** Override central attraction coefficient.
  220.      * @param templateOrbit template orbit
  221.      * @param provider gravity field provider
  222.      * @param positionAngleType position angle type to use
  223.      * @return orbit with overridden central attraction coefficient
  224.      */
  225.     private static Orbit overrideMu(final Orbit templateOrbit,
  226.                                     final UnnormalizedSphericalHarmonicsProvider provider,
  227.                                     final PositionAngleType positionAngleType) {
  228.         final double[] parameters    = new double[6];
  229.         final double[] parametersDot = parameters.clone();
  230.         templateOrbit.getType().mapOrbitToArray(templateOrbit, positionAngleType, parameters, parametersDot);
  231.         return templateOrbit.getType().mapArrayToOrbit(parameters, parametersDot, positionAngleType,
  232.                                                        templateOrbit.getDate(),
  233.                                                        provider.getMu(),
  234.                                                        templateOrbit.getFrame());
  235.     }

  237.     /** {@inheritDoc} */
  238.     public BrouwerLyddanePropagator buildPropagator(final double[] normalizedParameters) {
  239.         setParameters(normalizedParameters);

  241.         // Update M2 value and selection
  242.         double  newM2      = 0.0;
  243.         boolean isSelected = false;
  244.         for (final ParameterDriver driver : getPropagationParametersDrivers().getDrivers()) {
  245.             if (BrouwerLyddanePropagator.M2_NAME.equals(driver.getName())) {
  246.                 // it is OK as BL m2 parameterDriver has 1 value estimated from -INF to +INF, and
  247.                 // setPeriod method should not be called on this driver (to have several values estimated)
  248.                 newM2      = driver.getValue();
  249.                 isSelected = driver.isSelected();
  250.             }
  251.         }

  253.         // Initialize propagator
  254.         final BrouwerLyddanePropagator propagator = new BrouwerLyddanePropagator(createInitialOrbit(), getAttitudeProvider(), getMass(),
  255.             provider, newM2);
  256.         propagator.getParametersDrivers().get(0).setSelected(isSelected);
  257.         getImpulseManeuvers().forEach(propagator::addEventDetector);

  259.         // Return
  260.         return propagator;

  262.     }

  264.     /**
  265.      * Get the value of the M2 parameter.
  266.      * <p>
  267.      *  M2 represents the combination of all unmodeled secular along-track effects
  268.      *  (e.g. drag). It is usually fitted during an orbit determination.
  269.      * </p>
  270.      * @return the value of the M2 parameter
  271.      * @since 12.2
  272.      */
  273.     public double getM2Value() {
  274.         double m2 = 0.0;
  275.         for (final ParameterDriver driver : getPropagationParametersDrivers().getDrivers()) {
  276.             if (BrouwerLyddanePropagator.M2_NAME.equals(driver.getName())) {
  277.                 m2 = driver.getValue();
  278.             }
  279.         }
  280.         return m2;
  281.     }
  282. }
Created with JaCoCo 0.8.12.202403310830