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