BrouwerLyddanePropagatorBuilder.java

/* Copyright 2002-2026 CS GROUP
 * Licensed to CS GROUP (CS) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * CS licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.orekit.propagation.conversion;

import java.util.Collections;
import org.hipparchus.util.FastMath;
import org.orekit.attitudes.AttitudeProvider;
import org.orekit.attitudes.FrameAlignedProvider;
import org.orekit.forces.gravity.potential.GravityFieldFactory;
import org.orekit.forces.gravity.potential.TideSystem;
import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
import org.orekit.orbits.AbstractOrbitFactory;
import org.orekit.orbits.Orbit;
import org.orekit.propagation.Propagator;
import org.orekit.propagation.analytical.BrouwerLyddanePropagator;
import org.orekit.propagation.analytical.tle.TLE;
import org.orekit.utils.ParameterDriver;
import org.orekit.utils.ParameterDriversList;

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

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

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

    /** Build a new instance.
     * @param factory factory for initial orbit
     * @param provider for un-normalized zonal coefficients
     * @param M2 value of empirical drag coefficient in rad/s².
     *        If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
     * @see #BrouwerLyddanePropagatorBuilder(AbstractOrbitFactory,
     * UnnormalizedSphericalHarmonicsProvider, AttitudeProvider, double)
     * @since 14.0
     */
    public BrouwerLyddanePropagatorBuilder(final AbstractOrbitFactory<? extends Orbit> factory,
                                           final UnnormalizedSphericalHarmonicsProvider provider,
                                           final double M2) {
        this(factory, provider, FrameAlignedProvider.of(factory.getFrame()), M2);
    }

    /** Build a new instance.
     * @param factory factory for initial orbit
     * @param referenceRadius reference radius of the Earth for the potential model (m)
     * @param mu central attraction coefficient (m³/s²)
     * @param tideSystem tide system
     * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
     * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
     * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
     * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
     * @param M2 value of empirical drag coefficient in rad/s².
     *        If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
     * @see #BrouwerLyddanePropagatorBuilder(AbstractOrbitFactory,
     * UnnormalizedSphericalHarmonicsProvider, AttitudeProvider, double)
     * @since 14.0
     */
    public BrouwerLyddanePropagatorBuilder(final AbstractOrbitFactory<? extends Orbit> factory,
                                           final double referenceRadius,
                                           final double mu,
                                           final TideSystem tideSystem,
                                           final double c20,
                                           final double c30,
                                           final double c40,
                                           final double c50,
                                           final double M2) {
        this(factory,
             GravityFieldFactory.getUnnormalizedProvider(referenceRadius, mu, tideSystem,
                                                         new double[][] {
                                                             {
                                                                 0
                                                             }, {
                                                                 0
                                                             }, {
                                                                 c20
                                                             }, {
                                                                 c30
                                                             }, {
                                                                 c40
                                                             }, {
                                                                 c50
                                                             }
                                                         }, new double[][] {
                                                             {
                                                                 0
                                                             }, {
                                                                 0
                                                             }, {
                                                                 0
                                                             }, {
                                                                 0
                                                             }, {
                                                                 0
                                                             }, {
                                                                 0
                                                             }
                                                         }),
                M2);
    }

    /** Build a new instance.
     * @param factory factory for initial orbit
     * @param provider for un-normalized zonal coefficients
     * @param attitudeProvider attitude law to use
     * @param M2 value of empirical drag coefficient in rad/s².
     *        If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
     */
    public BrouwerLyddanePropagatorBuilder(final AbstractOrbitFactory<? extends Orbit> factory,
                                           final UnnormalizedSphericalHarmonicsProvider provider,
                                           final AttitudeProvider attitudeProvider,
                                           final double M2) {
        super((AbstractOrbitFactory<Orbit>) factory, true, attitudeProvider, Propagator.DEFAULT_MASS);
        factory.setMu(provider.getMu());
        this.provider = provider;
        // initialize M2 driver
        final ParameterDriver M2Driver = new ParameterDriver(BrouwerLyddanePropagator.M2_NAME, M2, SCALE,
                                                             Double.NEGATIVE_INFINITY,
                                                             Double.POSITIVE_INFINITY);
        addPropagationParameters(Collections.singletonList(M2Driver));
    }

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

        // Use cloned builder to unlink orbital drivers
        final BrouwerLyddanePropagatorBuilder builder =
            new BrouwerLyddanePropagatorBuilder((AbstractOrbitFactory<? extends Orbit>) clonedBuilder.getOrbitalParameterFactory().clone(),
                                                clonedBuilder.provider,
                                                clonedBuilder.getAttitudeProvider(),
                                                clonedBuilder.getM2Value());

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

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

        // Return cloned builder
        return builder;
    }

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

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

        // Initialize propagator
        final BrouwerLyddanePropagator propagator =
            new BrouwerLyddanePropagator(getOrbitalParameterFactory().createFromDrivers(),
                                         getAttitudeProvider(), getMass(),
            provider, newM2);
        propagator.getParametersDrivers().getFirst().setSelected(isSelected);
        getImpulseManeuvers().forEach(propagator::addEventDetector);

        // Return
        return propagator;

    }

    /**
     * Get the value of the M2 parameter.
     * <p>
     *  M2 represents the combination of all un-modeled secular along-track effects
     *  (e.g. drag). It is usually fitted during an orbit determination.
     * </p>
     * @return the value of the M2 parameter
     * @since 12.2
     */
    public double getM2Value() {
        double m2 = 0.0;
        for (final ParameterDriver driver : getPropagationParametersDrivers().getDrivers()) {
            if (BrouwerLyddanePropagator.M2_NAME.equals(driver.getName())) {
                m2 = driver.getValue();
            }
        }
        return m2;
    }
}