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;
18
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
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> {
64
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);
72
73 /** Provider for un-normalized coefficients. */
74 private final UnnormalizedSphericalHarmonicsProvider provider;
75
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
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 }
182
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 }
218
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 }
236
237 /** {@inheritDoc} */
238 public BrouwerLyddanePropagator buildPropagator(final double[] normalizedParameters) {
239 setParameters(normalizedParameters);
240
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 }
252
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);
258
259 // Return
260 return propagator;
261
262 }
263
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 }