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;
  33. import org.orekit.utils.ParameterDriversList;

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

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

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

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

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

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

  220.     /** Copy constructor.
  221.      * @param builder builder to copy from
  222.      */
  223.     private BrouwerLyddanePropagatorBuilder(final BrouwerLyddanePropagatorBuilder builder) {
  224.         this(builder.createInitialOrbit(), builder.provider, builder.getPositionAngleType(),
  225.              builder.getPositionScale(), builder.getAttitudeProvider(), builder.getM2Value());
  226.     }

  228.     /** {@inheritDoc}. */
  229.     @Override
  230.     public BrouwerLyddanePropagatorBuilder clone() {
  231.         // Call to super clone() method to avoid warning
  232.         final BrouwerLyddanePropagatorBuilder clonedBuilder = (BrouwerLyddanePropagatorBuilder) super.clone();

  234.         // Use copy constructor to unlink orbital drivers
  235.         final BrouwerLyddanePropagatorBuilder builder = new BrouwerLyddanePropagatorBuilder(clonedBuilder);

  237.         // Set mass
  238.         builder.setMass(getMass());

  240.         // Ensure drivers' selection consistency
  241.         final ParameterDriversList propDrivers = clonedBuilder.getPropagationParametersDrivers();
  242.         builder.getPropagationParametersDrivers().getDrivers().
  243.                         forEach(driver -> driver.setSelected(propDrivers.findByName(driver.getName()).isSelected()));

  245.         // Return cloned builder
  246.         return builder;
  247.     }

  249.     /** Override central attraction coefficient.
  250.      * @param templateOrbit template orbit
  251.      * @param provider gravity field provider
  252.      * @param positionAngleType position angle type to use
  253.      * @return orbit with overridden central attraction coefficient
  254.      */
  255.     private static Orbit overrideMu(final Orbit templateOrbit,
  256.                                     final UnnormalizedSphericalHarmonicsProvider provider,
  257.                                     final PositionAngleType positionAngleType) {
  258.         final double[] parameters    = new double[6];
  259.         final double[] parametersDot = parameters.clone();
  260.         templateOrbit.getType().mapOrbitToArray(templateOrbit, positionAngleType, parameters, parametersDot);
  261.         return templateOrbit.getType().mapArrayToOrbit(parameters, parametersDot, positionAngleType,
  262.                                                        templateOrbit.getDate(),
  263.                                                        provider.getMu(),
  264.                                                        templateOrbit.getFrame());
  265.     }

  267.     /** {@inheritDoc} */
  268.     public BrouwerLyddanePropagator buildPropagator(final double[] normalizedParameters) {
  269.         setParameters(normalizedParameters);

  271.         // Update M2 value and selection
  272.         double  newM2      = 0.0;
  273.         boolean isSelected = false;
  274.         for (final ParameterDriver driver : getPropagationParametersDrivers().getDrivers()) {
  275.             if (BrouwerLyddanePropagator.M2_NAME.equals(driver.getName())) {
  276.                 // it is OK as BL m2 parameterDriver has 1 value estimated from -INF to +INF, and
  277.                 // setPeriod method should not be called on this driver (to have several values estimated)
  278.                 newM2      = driver.getValue();
  279.                 isSelected = driver.isSelected();
  280.             }
  281.         }

  283.         // Initialize propagator
  284.         final BrouwerLyddanePropagator propagator = new BrouwerLyddanePropagator(createInitialOrbit(), getAttitudeProvider(), getMass(),
  285.             provider, newM2);
  286.         propagator.getParametersDrivers().get(0).setSelected(isSelected);
  287.         getImpulseManeuvers().forEach(propagator::addEventDetector);

  289.         // Return
  290.         return propagator;

  292.     }

  294.     /**
  295.      * Get the value of the M2 parameter.
  296.      * <p>
  297.      *  M2 represents the combination of all un-modeled secular along-track effects
  298.      *  (e.g. drag). It is usually fitted during an orbit determination.
  299.      * </p>
  300.      * @return the value of the M2 parameter
  301.      * @since 12.2
  302.      */
  303.     public double getM2Value() {
  304.         double m2 = 0.0;
  305.         for (final ParameterDriver driver : getPropagationParametersDrivers().getDrivers()) {
  306.             if (BrouwerLyddanePropagator.M2_NAME.equals(driver.getName())) {
  307.                 m2 = driver.getValue();
  308.             }
  309.         }
  310.         return m2;
  311.     }
  312. }
Created with JaCoCo 0.8.12.202403310830