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 import org.orekit.utils.ParameterDriversList;
34
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> {
65
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);
73
74 /** Provider for un-normalized coefficients. */
75 private final UnnormalizedSphericalHarmonicsProvider provider;
76
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 }
108
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 }
183
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 }
219
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 }
227
228 /** {@inheritDoc}. */
229 @Override
230 public BrouwerLyddanePropagatorBuilder clone() {
231 // Call to super clone() method to avoid warning
232 final BrouwerLyddanePropagatorBuilder clonedBuilder = (BrouwerLyddanePropagatorBuilder) super.clone();
233
234 // Use copy constructor to unlink orbital drivers
235 final BrouwerLyddanePropagatorBuilder builder = new BrouwerLyddanePropagatorBuilder(clonedBuilder);
236
237 // Set mass
238 builder.setMass(getMass());
239
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()));
244
245 // Return cloned builder
246 return builder;
247 }
248
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 }
266
267 /** {@inheritDoc} */
268 public BrouwerLyddanePropagator buildPropagator(final double[] normalizedParameters) {
269 setParameters(normalizedParameters);
270
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 }
282
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);
288
289 // Return
290 return propagator;
291
292 }
293
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 }