BrouwerLyddanePropagatorBuilder.java

  1. /* Copyright 2002-2023 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;

  18. import java.util.Collections;
  19. import java.util.List;
  20. import org.hipparchus.util.FastMath;
  21. import org.orekit.attitudes.AttitudeProvider;
  22. import org.orekit.attitudes.FrameAlignedProvider;
  23. import org.orekit.estimation.leastsquares.AbstractBatchLSModel;
  24. import org.orekit.estimation.leastsquares.BatchLSModel;
  25. import org.orekit.estimation.leastsquares.ModelObserver;
  26. import org.orekit.estimation.measurements.ObservedMeasurement;
  27. import org.orekit.forces.gravity.potential.GravityFieldFactory;
  28. import org.orekit.forces.gravity.potential.TideSystem;
  29. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
  30. import org.orekit.orbits.Orbit;
  31. import org.orekit.orbits.OrbitType;
  32. import org.orekit.orbits.PositionAngleType;
  33. import org.orekit.propagation.analytical.BrouwerLyddanePropagator;
  34. import org.orekit.propagation.analytical.tle.TLE;
  35. import org.orekit.utils.ParameterDriver;
  36. import org.orekit.utils.ParameterDriversList;

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

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

  74.     /** Provider for un-normalized coefficients. */
  75.     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.     }

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

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

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

  232.     /** {@inheritDoc} */
  233.     @Override
  234.     public BrouwerLyddanePropagatorBuilder copy() {

  235.         // Find M2 value
  236.         double m2 = 0.0;
  237.         for (final ParameterDriver driver : getPropagationParametersDrivers().getDrivers()) {
  238.             if (BrouwerLyddanePropagator.M2_NAME.equals(driver.getName())) {
  239.                 // it is OK as BL m2 parameterDriver has 1 value estimated from -INF to +INF, and
  240.                 // setPeriod method should not be called on this driver (to have several values estimated)
  241.                 m2 = driver.getValue();
  242.             }
  243.         }

  244.         return new BrouwerLyddanePropagatorBuilder(createInitialOrbit(), provider, getPositionAngleType(),
  245.                                                    getPositionScale(), getAttitudeProvider(), m2);
  246.     }

  247.     /** {@inheritDoc} */
  248.     public BrouwerLyddanePropagator buildPropagator(final double[] normalizedParameters) {
  249.         setParameters(normalizedParameters);

  250.         // Update M2 value and selection
  251.         double  newM2      = 0.0;
  252.         boolean isSelected = false;
  253.         for (final ParameterDriver driver : getPropagationParametersDrivers().getDrivers()) {
  254.             if (BrouwerLyddanePropagator.M2_NAME.equals(driver.getName())) {
  255.                 // it is OK as BL m2 parameterDriver has 1 value estimated from -INF to +INF, and
  256.                 // setPeriod method should not be called on this driver (to have several values estimated)
  257.                 newM2      = driver.getValue();
  258.                 isSelected = driver.isSelected();
  259.             }
  260.         }

  261.         // Initialize propagator
  262.         final BrouwerLyddanePropagator propagator = new BrouwerLyddanePropagator(createInitialOrbit(), getAttitudeProvider(), provider, newM2);
  263.         propagator.getParametersDrivers().get(0).setSelected(isSelected);

  264.         // Return
  265.         return propagator;

  266.     }

  267.     /** {@inheritDoc} */
  268.     @Override
  269.     public AbstractBatchLSModel buildLeastSquaresModel(final PropagatorBuilder[] builders,
  270.                                                        final List<ObservedMeasurement<?>> measurements,
  271.                                                        final ParameterDriversList estimatedMeasurementsParameters,
  272.                                                        final ModelObserver observer) {
  273.         return new BatchLSModel(builders, measurements, estimatedMeasurementsParameters, observer);
  274.     }

  275. }