1 /* Copyright 2002-2024 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 java.util.List;
21 import org.hipparchus.util.FastMath;
22 import org.orekit.attitudes.AttitudeProvider;
23 import org.orekit.attitudes.FrameAlignedProvider;
24 import org.orekit.estimation.leastsquares.AbstractBatchLSModel;
25 import org.orekit.estimation.leastsquares.BatchLSModel;
26 import org.orekit.estimation.leastsquares.ModelObserver;
27 import org.orekit.estimation.measurements.ObservedMeasurement;
28 import org.orekit.forces.gravity.potential.GravityFieldFactory;
29 import org.orekit.forces.gravity.potential.TideSystem;
30 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
31 import org.orekit.orbits.Orbit;
32 import org.orekit.orbits.OrbitType;
33 import org.orekit.orbits.PositionAngleType;
34 import org.orekit.propagation.analytical.BrouwerLyddanePropagator;
35 import org.orekit.propagation.analytical.tle.TLE;
36 import org.orekit.utils.ParameterDriver;
37 import org.orekit.utils.ParameterDriversList;
38
39 /** Builder for Brouwer-Lyddane propagator.
40 * <p>
41 * By default, Brouwer-Lyddane model considers only the perturbations due to zonal harmonics.
42 * However, for low Earth orbits, the magnitude of the perturbative acceleration due to
43 * atmospheric drag can be significant. Warren Phipps' 1992 thesis considered the atmospheric
44 * drag by time derivatives of the <i>mean</i> mean anomaly using the catch-all coefficient M2.
45 *
46 * Usually, M2 is adjusted during an orbit determination process and it represents the
47 * combination of all unmodeled secular along-track effects (i.e. not just the atmospheric drag).
48 * The behavior of M2 is closed to the {@link TLE#getBStar()} parameter for the TLE.
49 *
50 * If the value of M2 is equal to {@link BrouwerLyddanePropagator#M2 0.0}, the along-track
51 * secular effects are not considered in the dynamical model. Typical values for M2 are not known.
52 * 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).
53 * The unit of M2 is rad/s².
54 * <p>
55 * To estimate the M2 parameter, it is necessary to call the {@link #getPropagationParametersDrivers()} method
56 * as follow:
57 * <pre>
58 * for (ParameterDriver driver : builder.getPropagationParametersDrivers().getDrivers()) {
59 * if (BrouwerLyddanePropagator.M2_NAME.equals(driver.getName())) {
60 * driver.setSelected(true);
61 * }
62 * }
63 * </pre>
64 * @author Melina Vanel
65 * @author Bryan Cazabonne
66 * @since 11.1
67 */
68 public class BrouwerLyddanePropagatorBuilder extends AbstractPropagatorBuilder {
69
70 /** Parameters scaling factor.
71 * <p>
72 * We use a power of 2 to avoid numeric noise introduction
73 * in the multiplications/divisions sequences.
74 * </p>
75 */
76 private static final double SCALE = FastMath.scalb(1.0, -32);
77
78 /** Provider for un-normalized coefficients. */
79 private final UnnormalizedSphericalHarmonicsProvider provider;
80
81 /** Build a new instance.
82 * <p>
83 * The template orbit is used as a model to {@link
84 * #createInitialOrbit() create initial orbit}. It defines the
85 * inertial frame, the central attraction coefficient, the orbit type, and is also
86 * used together with the {@code positionScale} to convert from the {@link
87 * org.orekit.utils.ParameterDriver#setNormalizedValue(double) normalized} parameters
88 * used by the callers of this builder to the real orbital parameters.
89 * The default attitude provider is aligned with the orbit's inertial frame.
90 * </p>
91 *
92 * @param templateOrbit reference orbit from which real orbits will be built
93 * (note that the mu from this orbit will be overridden with the mu from the
94 * {@code provider})
95 * @param provider for un-normalized zonal coefficients
96 * @param positionAngleType position angle type to use
97 * @param positionScale scaling factor used for orbital parameters normalization
98 * (typically set to the expected standard deviation of the position)
99 * @param M2 value of empirical drag coefficient in rad/s².
100 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
101 * @see #BrouwerLyddanePropagatorBuilder(Orbit,
102 * UnnormalizedSphericalHarmonicsProvider, PositionAngleType, double, AttitudeProvider, double)
103 */
104 public BrouwerLyddanePropagatorBuilder(final Orbit templateOrbit,
105 final UnnormalizedSphericalHarmonicsProvider provider,
106 final PositionAngleType positionAngleType,
107 final double positionScale,
108 final double M2) {
109 this(templateOrbit, provider, positionAngleType, positionScale,
110 FrameAlignedProvider.of(templateOrbit.getFrame()), M2);
111 }
112
113 /** Build a new instance.
114 * <p>
115 * The template orbit is used as a model to {@link
116 * #createInitialOrbit() create initial orbit}. It defines the
117 * inertial frame, the central attraction coefficient, the orbit type, and is also
118 * used together with the {@code positionScale} to convert from the {@link
119 * org.orekit.utils.ParameterDriver#setNormalizedValue(double) normalized} parameters used by the
120 * callers of this builder to the real orbital parameters.
121 * </p>
122 *
123 * @param templateOrbit reference orbit from which real orbits will be built
124 * (note that the mu from this orbit will be overridden with the mu from the
125 * {@code provider})
126 * @param referenceRadius reference radius of the Earth for the potential model (m)
127 * @param mu central attraction coefficient (m³/s²)
128 * @param tideSystem tide system
129 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
130 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
131 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
132 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
133 * @param orbitType orbit type to use
134 * @param positionAngleType position angle type to use
135 * @param positionScale scaling factor used for orbital parameters normalization
136 * (typically set to the expected standard deviation of the position)
137 * @param M2 value of empirical drag coefficient in rad/s².
138 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
139 * @see #BrouwerLyddanePropagatorBuilder(Orbit,
140 * UnnormalizedSphericalHarmonicsProvider, PositionAngleType, double, AttitudeProvider, double)
141 */
142 public BrouwerLyddanePropagatorBuilder(final Orbit templateOrbit,
143 final double referenceRadius,
144 final double mu,
145 final TideSystem tideSystem,
146 final double c20,
147 final double c30,
148 final double c40,
149 final double c50,
150 final OrbitType orbitType,
151 final PositionAngleType positionAngleType,
152 final double positionScale,
153 final double M2) {
154 this(templateOrbit,
155 GravityFieldFactory.getUnnormalizedProvider(referenceRadius, mu, tideSystem,
156 new double[][] {
157 {
158 0
159 }, {
160 0
161 }, {
162 c20
163 }, {
164 c30
165 }, {
166 c40
167 }, {
168 c50
169 }
170 }, new double[][] {
171 {
172 0
173 }, {
174 0
175 }, {
176 0
177 }, {
178 0
179 }, {
180 0
181 }, {
182 0
183 }
184 }),
185 positionAngleType, positionScale, M2);
186 }
187
188 /** Build a new instance.
189 * <p>
190 * The template orbit is used as a model to {@link
191 * #createInitialOrbit() create initial orbit}. It defines the
192 * inertial frame, the central attraction coefficient, the orbit type, and is also
193 * used together with the {@code positionScale} to convert from the {@link
194 * org.orekit.utils.ParameterDriver#setNormalizedValue(double) normalized} parameters used by the
195 * callers of this builder to the real orbital parameters.
196 * </p>
197 * @param templateOrbit reference orbit from which real orbits will be built
198 * (note that the mu from this orbit will be overridden with the mu from the
199 * {@code provider})
200 * @param provider for un-normalized zonal coefficients
201 * @param positionAngleType position angle type to use
202 * @param positionScale scaling factor used for orbital parameters normalization
203 * (typically set to the expected standard deviation of the position)
204 * @param attitudeProvider attitude law to use
205 * @param M2 value of empirical drag coefficient in rad/s².
206 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
207 */
208 public BrouwerLyddanePropagatorBuilder(final Orbit templateOrbit,
209 final UnnormalizedSphericalHarmonicsProvider provider,
210 final PositionAngleType positionAngleType,
211 final double positionScale,
212 final AttitudeProvider attitudeProvider,
213 final double M2) {
214 super(overrideMu(templateOrbit, provider, positionAngleType), positionAngleType, positionScale, true, attitudeProvider);
215 this.provider = provider;
216 // initialize M2 driver
217 final ParameterDriver M2Driver = new ParameterDriver(BrouwerLyddanePropagator.M2_NAME, M2, SCALE,
218 Double.NEGATIVE_INFINITY,
219 Double.POSITIVE_INFINITY);
220 addSupportedParameters(Collections.singletonList(M2Driver));
221 }
222
223 /** Override central attraction coefficient.
224 * @param templateOrbit template orbit
225 * @param provider gravity field provider
226 * @param positionAngleType position angle type to use
227 * @return orbit with overridden central attraction coefficient
228 */
229 private static Orbit overrideMu(final Orbit templateOrbit,
230 final UnnormalizedSphericalHarmonicsProvider provider,
231 final PositionAngleType positionAngleType) {
232 final double[] parameters = new double[6];
233 final double[] parametersDot = templateOrbit.hasDerivatives() ? new double[6] : null;
234 templateOrbit.getType().mapOrbitToArray(templateOrbit, positionAngleType, parameters, parametersDot);
235 return templateOrbit.getType().mapArrayToOrbit(parameters, parametersDot, positionAngleType,
236 templateOrbit.getDate(),
237 provider.getMu(),
238 templateOrbit.getFrame());
239 }
240
241 /** {@inheritDoc} */
242 @Override
243 public BrouwerLyddanePropagatorBuilder copy() {
244
245 // Find M2 value
246 double m2 = 0.0;
247 for (final ParameterDriver driver : getPropagationParametersDrivers().getDrivers()) {
248 if (BrouwerLyddanePropagator.M2_NAME.equals(driver.getName())) {
249 // it is OK as BL m2 parameterDriver has 1 value estimated from -INF to +INF, and
250 // setPeriod method should not be called on this driver (to have several values estimated)
251 m2 = driver.getValue();
252 }
253 }
254
255 return new BrouwerLyddanePropagatorBuilder(createInitialOrbit(), provider, getPositionAngleType(),
256 getPositionScale(), getAttitudeProvider(), m2);
257 }
258
259 /** {@inheritDoc} */
260 public BrouwerLyddanePropagator buildPropagator(final double[] normalizedParameters) {
261 setParameters(normalizedParameters);
262
263 // Update M2 value and selection
264 double newM2 = 0.0;
265 boolean isSelected = false;
266 for (final ParameterDriver driver : getPropagationParametersDrivers().getDrivers()) {
267 if (BrouwerLyddanePropagator.M2_NAME.equals(driver.getName())) {
268 // it is OK as BL m2 parameterDriver has 1 value estimated from -INF to +INF, and
269 // setPeriod method should not be called on this driver (to have several values estimated)
270 newM2 = driver.getValue();
271 isSelected = driver.isSelected();
272 }
273 }
274
275 // Initialize propagator
276 final BrouwerLyddanePropagator propagator = new BrouwerLyddanePropagator(createInitialOrbit(), getAttitudeProvider(), provider, newM2);
277 propagator.getParametersDrivers().get(0).setSelected(isSelected);
278
279 // Return
280 return propagator;
281
282 }
283
284 /** {@inheritDoc} */
285 @Override
286 public AbstractBatchLSModel buildLeastSquaresModel(final PropagatorBuilder[] builders,
287 final List<ObservedMeasurement<?>> measurements,
288 final ParameterDriversList estimatedMeasurementsParameters,
289 final ModelObserver observer) {
290 return new BatchLSModel(builders, measurements, estimatedMeasurementsParameters, observer);
291 }
292
293 }