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.forces.drag;
18
19 import org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.analysis.differentiation.DSFactory;
21 import org.hipparchus.analysis.differentiation.DerivativeStructure;
22 import org.hipparchus.analysis.differentiation.Gradient;
23 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
24 import org.hipparchus.geometry.euclidean.threed.Vector3D;
25 import org.orekit.forces.ForceModel;
26 import org.orekit.frames.Frame;
27 import org.orekit.frames.StaticTransform;
28 import org.orekit.models.earth.atmosphere.Atmosphere;
29 import org.orekit.propagation.FieldSpacecraftState;
30 import org.orekit.time.AbsoluteDate;
31 import org.orekit.time.FieldAbsoluteDate;
32 import org.orekit.utils.FieldPVCoordinates;
33
34 import java.util.Arrays;
35
36 /**
37 * Base class for drag force models.
38 * @see DragForce
39 * @author Bryan Cazabonne
40 * @since 10.2
41 */
42 public abstract class AbstractDragForceModel implements ForceModel {
43
44 /** Atmospheric model. */
45 private final Atmosphere atmosphere;
46
47 /**
48 * Flag to use (first-order) finite differences instead of automatic differentiation when computing density derivatives w.r.t. position.
49 */
50 private final boolean useFiniteDifferencesOnDensityWrtPosition;
51
52 /**
53 * Constructor with default value for finite differences flag.
54 * @param atmosphere atmospheric model
55 */
56 protected AbstractDragForceModel(final Atmosphere atmosphere) {
57 this(atmosphere, true);
58 }
59
60 /**
61 * Constructor.
62 * @param atmosphere atmospheric model
63 * @param useFiniteDifferencesOnDensityWrtPosition flag to use finite differences to compute density derivatives w.r.t.
64 * position (is less accurate but can be faster depending on model)
65 * @since 12.1
66 */
67 protected AbstractDragForceModel(final Atmosphere atmosphere, final boolean useFiniteDifferencesOnDensityWrtPosition) {
68 this.atmosphere = atmosphere;
69 this.useFiniteDifferencesOnDensityWrtPosition = useFiniteDifferencesOnDensityWrtPosition;
70 }
71
72 /** Get the atmospheric model.
73 * @return atmosphere model
74 * @since 12.1
75 */
76 public Atmosphere getAtmosphere() {
77 return atmosphere;
78 }
79
80 /** {@inheritDoc} */
81 @Override
82 public boolean dependsOnPositionOnly() {
83 return false;
84 }
85
86 /** Check if a field state corresponds to derivatives with respect to state.
87 * @param state state to check
88 * @param <T> type of the field elements
89 * @return true if state corresponds to derivatives with respect to state
90 */
91 protected <T extends CalculusFieldElement<T>> boolean isDSStateDerivative(final FieldSpacecraftState<T> state) {
92 try {
93 final DerivativeStructure dsMass = (DerivativeStructure) state.getMass();
94 final int o = dsMass.getOrder();
95 final int p = dsMass.getFreeParameters();
96
97 // To be in the desired case:
98 // Order must be 1 (first order derivatives only)
99 // Number of parameters must be 6 (PV), 7 (PV + drag coefficient) or 8 (PV + drag coefficient + lift ratio)
100 if (o != 1 || p != 6 && p != 7 && p != 8) {
101 return false;
102 }
103
104 // Check that the first 6 parameters are position and velocity
105 @SuppressWarnings("unchecked")
106 final FieldPVCoordinates<DerivativeStructure> pv = (FieldPVCoordinates<DerivativeStructure>) state.getPVCoordinates();
107 return isVariable(pv.getPosition().getX(), 0) &&
108 isVariable(pv.getPosition().getY(), 1) &&
109 isVariable(pv.getPosition().getZ(), 2) &&
110 isVariable(pv.getVelocity().getX(), 3) &&
111 isVariable(pv.getVelocity().getY(), 4) &&
112 isVariable(pv.getVelocity().getZ(), 5);
113 } catch (ClassCastException cce) {
114 return false;
115 }
116 }
117
118 /** Check if a field state corresponds to derivatives with respect to state.
119 * @param state state to check
120 * @param <T> type of the field elements
121 * @return true if state corresponds to derivatives with respect to state
122 */
123 protected <T extends CalculusFieldElement<T>> boolean isGradientStateDerivative(final FieldSpacecraftState<T> state) {
124 try {
125 final Gradient gMass = (Gradient) state.getMass();
126 final int p = gMass.getFreeParameters();
127
128 // To be in the desired case:
129 // Number of parameters must be 6 (PV), 7 (PV + drag coefficient) or 8 (PV + drag coefficient + lift ratio)
130 if (p != 6 && p != 7 && p != 8) {
131 return false;
132 }
133
134 // Check that the first 6 parameters are position and velocity
135 @SuppressWarnings("unchecked")
136 final FieldPVCoordinates<Gradient> pv = (FieldPVCoordinates<Gradient>) state.getPVCoordinates();
137 return isVariable(pv.getPosition().getX(), 0) &&
138 isVariable(pv.getPosition().getY(), 1) &&
139 isVariable(pv.getPosition().getZ(), 2) &&
140 isVariable(pv.getVelocity().getX(), 3) &&
141 isVariable(pv.getVelocity().getY(), 4) &&
142 isVariable(pv.getVelocity().getZ(), 5);
143 } catch (ClassCastException cce) {
144 return false;
145 }
146 }
147
148 /**
149 * Evaluate the Field density.
150 * @param s spacecraft state
151 * @return atmospheric density
152 * @param <T> field type
153 * @since 12.1
154 */
155 @SuppressWarnings("unchecked")
156 protected <T extends CalculusFieldElement<T>> T getFieldDensity(final FieldSpacecraftState<T> s) {
157 final FieldAbsoluteDate<T> date = s.getDate();
158 final Frame frame = s.getFrame();
159 final FieldVector3D<T> position = s.getPosition();
160 if (isGradientStateDerivative(s)) {
161 if (useFiniteDifferencesOnDensityWrtPosition) {
162 return (T) this.getGradientDensityWrtStateUsingFiniteDifferences(date.toAbsoluteDate(), frame, (FieldVector3D<Gradient>) position);
163 } else {
164 return (T) this.getGradientDensityWrtState(date.toAbsoluteDate(), frame, (FieldVector3D<Gradient>) position);
165 }
166 } else if (isDSStateDerivative(s)) {
167 if (useFiniteDifferencesOnDensityWrtPosition) {
168 return (T) this.getDSDensityWrtStateUsingFiniteDifferences(date.toAbsoluteDate(), frame, (FieldVector3D<DerivativeStructure>) position);
169 } else {
170 return (T) this.getDSDensityWrtState(date.toAbsoluteDate(), frame, (FieldVector3D<DerivativeStructure>) position);
171 }
172 } else {
173 return atmosphere.getDensity(date, position, frame);
174 }
175 }
176
177 /** Check if a derivative represents a specified variable.
178 * @param ds derivative to check
179 * @param index index of the variable
180 * @return true if the derivative represents a specified variable
181 */
182 protected boolean isVariable(final DerivativeStructure ds, final int index) {
183 final double[] derivatives = ds.getAllDerivatives();
184 boolean check = true;
185 for (int i = 1; i < derivatives.length; ++i) {
186 check &= derivatives[i] == ((index + 1 == i) ? 1.0 : 0.0);
187 }
188 return check;
189 }
190
191 /** Check if a derivative represents a specified variable.
192 * @param g derivative to check
193 * @param index index of the variable
194 * @return true if the derivative represents a specified variable
195 */
196 protected boolean isVariable(final Gradient g, final int index) {
197 final double[] derivatives = g.getGradient();
198 boolean check = true;
199 for (int i = 0; i < derivatives.length; ++i) {
200 check &= derivatives[i] == ((index == i) ? 1.0 : 0.0);
201 }
202 return check;
203 }
204
205 /** Compute density and its derivatives.
206 * Using finite differences for the derivatives.
207 * And doing the actual computation only for the derivatives with respect to position (others are set to 0.).
208 * <p>
209 * From a theoretical point of view, this method computes the same values
210 * as {@link Atmosphere#getDensity(FieldAbsoluteDate, FieldVector3D, Frame)} in the
211 * specific case of {@link DerivativeStructure} with respect to state, so
212 * it is less general. However, it can be faster depending the Field implementation.
213 * </p>
214 * <p>
215 * The derivatives should be computed with respect to position. The input
216 * parameters already take into account the free parameters (6, 7 or 8 depending
217 * on derivation with respect to drag coefficient and lift ratio being considered or not)
218 * and order (always 1). Free parameters at indices 0, 1 and 2 correspond to derivatives
219 * with respect to position. Free parameters at indices 3, 4 and 5 correspond
220 * to derivatives with respect to velocity (these derivatives will remain zero
221 * as the atmospheric density does not depend on velocity). Free parameter
222 * at indexes 6 and 7 (if present) corresponds to derivatives with respect to drag coefficient
223 * and/or lift ratio (one of these or both).
224 * This 2 last derivatives will remain zero as atmospheric density does not depend on them.
225 * </p>
226 * @param date current date
227 * @param frame inertial reference frame for state (both orbit and attitude)
228 * @param position position of spacecraft in inertial frame
229 * @return the density and its derivatives
230 */
231 protected DerivativeStructure getDSDensityWrtStateUsingFiniteDifferences(final AbsoluteDate date,
232 final Frame frame,
233 final FieldVector3D<DerivativeStructure> position) {
234
235 // Retrieve derivation properties for parameter T
236 // It is implied here that T is a DerivativeStructure
237 // With order 1 and 6, 7 or 8 free parameters
238 // This is all checked before in method isStateDerivatives
239 final DSFactory factory = position.getX().getFactory();
240
241 // Build a DerivativeStructure using only derivatives with respect to position
242 final DSFactory factory3 = new DSFactory(3, 1);
243 final FieldVector3D<DerivativeStructure> position3 =
244 new FieldVector3D<>(factory3.variable(0, position.getX().getReal()),
245 factory3.variable(1, position.getY().getReal()),
246 factory3.variable(2, position.getZ().getReal()));
247
248 // Get atmosphere properties in atmosphere own frame
249 final Frame atmFrame = atmosphere.getFrame();
250 final StaticTransform toBody = frame.getStaticTransformTo(atmFrame, date);
251 final FieldVector3D<DerivativeStructure> posBodyDS = toBody.transformPosition(position3);
252 final Vector3D posBody = posBodyDS.toVector3D();
253
254 // Estimate density model by finite differences and composition
255 // Using a delta of 1m
256 final double delta = 1.0;
257 final double x = posBody.getX();
258 final double y = posBody.getY();
259 final double z = posBody.getZ();
260 final double rho0 = atmosphere.getDensity(date, posBody, atmFrame);
261 final double dRhodX = (atmosphere.getDensity(date, new Vector3D(x + delta, y, z), atmFrame) - rho0) / delta;
262 final double dRhodY = (atmosphere.getDensity(date, new Vector3D(x, y + delta, z), atmFrame) - rho0) / delta;
263 final double dRhodZ = (atmosphere.getDensity(date, new Vector3D(x, y, z + delta), atmFrame) - rho0) / delta;
264 final double[] dXdQ = posBodyDS.getX().getAllDerivatives();
265 final double[] dYdQ = posBodyDS.getY().getAllDerivatives();
266 final double[] dZdQ = posBodyDS.getZ().getAllDerivatives();
267
268 // Density with derivatives:
269 // - The value and only the 3 first derivatives (those with respect to spacecraft position) are computed
270 // - Others are set to 0.
271 final int p = factory.getCompiler().getFreeParameters();
272 final double[] rhoAll = new double[p + 1];
273 rhoAll[0] = rho0;
274 for (int i = 1; i < 4; ++i) {
275 rhoAll[i] = dRhodX * dXdQ[i] + dRhodY * dYdQ[i] + dRhodZ * dZdQ[i];
276 }
277
278 return factory.build(rhoAll);
279 }
280
281 /** Compute density and its derivatives.
282 * And doing the actual computation only for the derivatives with respect to position (others are set to 0.).
283 * <p>
284 * The derivatives should be computed with respect to position. The input
285 * parameters already take into account the free parameters (6, 7 or 8 depending
286 * on derivation with respect to drag coefficient and lift ratio being considered or not)
287 * and order (always 1). Free parameters at indices 0, 1 and 2 correspond to derivatives
288 * with respect to position. Free parameters at indices 3, 4 and 5 correspond
289 * to derivatives with respect to velocity (these derivatives will remain zero
290 * as the atmospheric density does not depend on velocity). Free parameter
291 * at indexes 6 and 7 (if present) corresponds to derivatives with respect to drag coefficient
292 * and/or lift ratio (one of these or both).
293 * This 2 last derivatives will remain zero as atmospheric density does not depend on them.
294 * </p>
295 * @param date current date
296 * @param frame inertial reference frame for state (both orbit and attitude)
297 * @param position position of spacecraft in inertial frame
298 * @return the density and its derivatives
299 */
300 protected DerivativeStructure getDSDensityWrtState(final AbsoluteDate date, final Frame frame,
301 final FieldVector3D<DerivativeStructure> position) {
302
303 // Retrieve derivation properties for parameter T
304 // It is implied here that T is a DerivativeStructure
305 // With order 1 and 6, 7 or 8 free parameters
306 // This is all checked before in method isStateDerivatives
307 final DSFactory factory = position.getX().getFactory();
308
309 // Build a DerivativeStructure using only derivatives with respect to position
310 final DSFactory factory3 = new DSFactory(3, 1);
311 final FieldVector3D<DerivativeStructure> position3 =
312 new FieldVector3D<>(factory3.variable(0, position.getX().getReal()),
313 factory3.variable(1, position.getY().getReal()),
314 factory3.variable(2, position.getZ().getReal()));
315
316 // Get atmosphere properties in atmosphere own frame
317 final Frame atmFrame = atmosphere.getFrame();
318 final StaticTransform toBody = frame.getStaticTransformTo(atmFrame, date);
319 final FieldVector3D<DerivativeStructure> posBodyDS = toBody.transformPosition(position3);
320 final FieldAbsoluteDate<DerivativeStructure> fieldDate = new FieldAbsoluteDate<>(position3.getX().getField(), date);
321 final DerivativeStructure density = atmosphere.getDensity(fieldDate, posBodyDS, atmFrame);
322
323 // Density with derivatives:
324 // - The value and only the 3 first derivatives (those with respect to spacecraft position) are computed
325 // - Others are set to 0.
326 final double[] derivatives = Arrays.copyOf(density.getAllDerivatives(), factory.getCompiler().getSize());
327 return factory.build(derivatives);
328 }
329
330 /** Compute density and its derivatives.
331 * Using finite differences for the derivatives.
332 * And doing the actual computation only for the derivatives with respect to position (others are set to 0.).
333 * <p>
334 * From a theoretical point of view, this method computes the same values
335 * as {@link Atmosphere#getDensity(FieldAbsoluteDate, FieldVector3D, Frame)} in the
336 * specific case of {@link Gradient} with respect to state, so
337 * it is less general. However, it can be faster depending the Field implementation.
338 * </p>
339 * <p>
340 * The derivatives should be computed with respect to position. The input
341 * parameters already take into account the free parameters (6, 7 or 8 depending
342 * on derivation with respect to drag coefficient and lift ratio being considered or not)
343 * and order (always 1). Free parameters at indices 0, 1 and 2 correspond to derivatives
344 * with respect to position. Free parameters at indices 3, 4 and 5 correspond
345 * to derivatives with respect to velocity (these derivatives will remain zero
346 * as the atmospheric density does not depend on velocity). Free parameter
347 * at indexes 6 and 7 (if present) corresponds to derivatives with respect to drag coefficient
348 * and/or lift ratio (one of these or both).
349 * This 2 last derivatives will remain zero as atmospheric density does not depend on them.
350 * </p>
351 * @param date current date
352 * @param frame inertial reference frame for state (both orbit and attitude)
353 * @param position position of spacecraft in inertial frame
354 * @return the density and its derivatives
355 */
356 protected Gradient getGradientDensityWrtStateUsingFiniteDifferences(final AbsoluteDate date,
357 final Frame frame,
358 final FieldVector3D<Gradient> position) {
359
360 // Build a Gradient using only derivatives with respect to position
361 final FieldVector3D<Gradient> position3 =
362 new FieldVector3D<>(Gradient.variable(3, 0, position.getX().getReal()),
363 Gradient.variable(3, 1, position.getY().getReal()),
364 Gradient.variable(3, 2, position.getZ().getReal()));
365
366 // Get atmosphere properties in atmosphere own frame
367 final Frame atmFrame = atmosphere.getFrame();
368 final StaticTransform toBody = frame.getStaticTransformTo(atmFrame, date);
369 final FieldVector3D<Gradient> posBodyDS = toBody.transformPosition(position3);
370 final Vector3D posBody = posBodyDS.toVector3D();
371
372 // Estimate density model by finite differences and composition
373 // Using a delta of 1m
374 final double delta = 1.0;
375 final double x = posBody.getX();
376 final double y = posBody.getY();
377 final double z = posBody.getZ();
378 final double rho0 = atmosphere.getDensity(date, posBody, atmFrame);
379 final double dRhodX = (atmosphere.getDensity(date, new Vector3D(x + delta, y, z), atmFrame) - rho0) / delta;
380 final double dRhodY = (atmosphere.getDensity(date, new Vector3D(x, y + delta, z), atmFrame) - rho0) / delta;
381 final double dRhodZ = (atmosphere.getDensity(date, new Vector3D(x, y, z + delta), atmFrame) - rho0) / delta;
382 final double[] dXdQ = posBodyDS.getX().getGradient();
383 final double[] dYdQ = posBodyDS.getY().getGradient();
384 final double[] dZdQ = posBodyDS.getZ().getGradient();
385
386 // Density with derivatives:
387 // - The value and only the 3 first derivatives (those with respect to spacecraft position) are computed
388 // - Others are set to 0.
389 final int p = position.getX().getFreeParameters();
390 final double[] rhoAll = new double[p];
391 for (int i = 0; i < 3; ++i) {
392 rhoAll[i] = dRhodX * dXdQ[i] + dRhodY * dYdQ[i] + dRhodZ * dZdQ[i];
393 }
394
395 return new Gradient(rho0, rhoAll);
396 }
397
398 /** Compute density and its derivatives.
399 * <p>
400 * The derivatives should be computed with respect to position. The input
401 * parameters already take into account the free parameters (6, 7 or 8 depending
402 * on derivation with respect to drag coefficient and lift ratio being considered or not)
403 * and order (always 1). Free parameters at indices 0, 1 and 2 correspond to derivatives
404 * with respect to position. Free parameters at indices 3, 4 and 5 correspond
405 * to derivatives with respect to velocity (these derivatives will remain zero
406 * as the atmospheric density does not depend on velocity). Free parameter
407 * at indexes 6 and 7 (if present) corresponds to derivatives with respect to drag coefficient
408 * and/or lift ratio (one of these or both).
409 * This 2 last derivatives will remain zero as atmospheric density does not depend on them.
410 * </p>
411 * @param date current date
412 * @param frame inertial reference frame for state (both orbit and attitude)
413 * @param position position of spacecraft in inertial frame
414 * @return the density and its derivatives
415 * @since 12.1
416 */
417 protected Gradient getGradientDensityWrtState(final AbsoluteDate date, final Frame frame,
418 final FieldVector3D<Gradient> position) {
419
420 // Build a Gradient using only derivatives with respect to position
421 final int positionDimension = 3;
422 final FieldVector3D<Gradient> position3 =
423 new FieldVector3D<>(Gradient.variable(positionDimension, 0, position.getX().getReal()),
424 Gradient.variable(positionDimension, 1, position.getY().getReal()),
425 Gradient.variable(positionDimension, 2, position.getZ().getReal()));
426
427 // Get atmosphere properties in atmosphere own frame
428 final Frame atmFrame = atmosphere.getFrame();
429 final StaticTransform toBody = frame.getStaticTransformTo(atmFrame, date);
430 final FieldVector3D<Gradient> posBodyGradient = toBody.transformPosition(position3);
431 final FieldAbsoluteDate<Gradient> fieldDate = new FieldAbsoluteDate<>(position3.getX().getField(), date);
432 final Gradient density = atmosphere.getDensity(fieldDate, posBodyGradient, atmFrame);
433
434 // Density with derivatives:
435 // - The value and only the 3 first derivatives (those with respect to spacecraft position) are computed
436 // - Others are set to 0.
437 final double[] derivatives = Arrays.copyOf(density.getGradient(), position.getX().getFreeParameters());
438 return new Gradient(density.getValue(), derivatives);
439 }
440 }