1 /* Copyright 2002-2021 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.AbstractForceModel;
26 import org.orekit.frames.Frame;
27 import org.orekit.frames.Transform;
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 /**
35 * Base class for drag force models.
36 * @see DragForce
37 * @see TimeSpanDragForce
38 * @author Bryan Cazabonne
39 * @since 10.2
40 */
41 public abstract class AbstractDragForceModel extends AbstractForceModel {
42
43 /** Atmospheric model. */
44 private final Atmosphere atmosphere;
45
46 /**
47 * Constructor.
48 * @param atmosphere atmospheric model
49 */
50 protected AbstractDragForceModel(final Atmosphere atmosphere) {
51 this.atmosphere = atmosphere;
52 }
53
54 /** {@inheritDoc} */
55 @Override
56 public boolean dependsOnPositionOnly() {
57 return false;
58 }
59
60 /** Check if a field state corresponds to derivatives with respect to state.
61 * @param state state to check
62 * @param <T> type of the field elements
63 * @return true if state corresponds to derivatives with respect to state
64 */
65 protected <T extends CalculusFieldElement<T>> boolean isDSStateDerivative(final FieldSpacecraftState<T> state) {
66 try {
67 final DerivativeStructure dsMass = (DerivativeStructure) state.getMass();
68 final int o = dsMass.getOrder();
69 final int p = dsMass.getFreeParameters();
70
71 // To be in the desired case:
72 // Order must be 1 (first order derivatives only)
73 // Number of parameters must be 6 (PV), 7 (PV + drag coefficient) or 8 (PV + drag coefficient + lift ratio)
74 if (o != 1 || p != 6 && p != 7 && p != 8) {
75 return false;
76 }
77
78 // Check that the first 6 parameters are position and velocity
79 @SuppressWarnings("unchecked")
80 final FieldPVCoordinates<DerivativeStructure> pv = (FieldPVCoordinates<DerivativeStructure>) state.getPVCoordinates();
81 return isVariable(pv.getPosition().getX(), 0) &&
82 isVariable(pv.getPosition().getY(), 1) &&
83 isVariable(pv.getPosition().getZ(), 2) &&
84 isVariable(pv.getVelocity().getX(), 3) &&
85 isVariable(pv.getVelocity().getY(), 4) &&
86 isVariable(pv.getVelocity().getZ(), 5);
87 } catch (ClassCastException cce) {
88 return false;
89 }
90 }
91
92 /** Check if a field state corresponds to derivatives with respect to state.
93 * @param state state to check
94 * @param <T> type of the field elements
95 * @return true if state corresponds to derivatives with respect to state
96 */
97 protected <T extends CalculusFieldElement<T>> boolean isGradientStateDerivative(final FieldSpacecraftState<T> state) {
98 try {
99 final Gradient gMass = (Gradient) state.getMass();
100 final int p = gMass.getFreeParameters();
101
102 // To be in the desired case:
103 // Number of parameters must be 6 (PV), 7 (PV + drag coefficient) or 8 (PV + drag coefficient + lift ratio)
104 if (p != 6 && p != 7 && p != 8) {
105 return false;
106 }
107
108 // Check that the first 6 parameters are position and velocity
109 @SuppressWarnings("unchecked")
110 final FieldPVCoordinates<Gradient> pv = (FieldPVCoordinates<Gradient>) state.getPVCoordinates();
111 return isVariable(pv.getPosition().getX(), 0) &&
112 isVariable(pv.getPosition().getY(), 1) &&
113 isVariable(pv.getPosition().getZ(), 2) &&
114 isVariable(pv.getVelocity().getX(), 3) &&
115 isVariable(pv.getVelocity().getY(), 4) &&
116 isVariable(pv.getVelocity().getZ(), 5);
117 } catch (ClassCastException cce) {
118 return false;
119 }
120 }
121
122 /** Check if a derivative represents a specified variable.
123 * @param ds derivative to check
124 * @param index index of the variable
125 * @return true if the derivative represents a specified variable
126 */
127 protected boolean isVariable(final DerivativeStructure ds, final int index) {
128 final double[] derivatives = ds.getAllDerivatives();
129 boolean check = true;
130 for (int i = 1; i < derivatives.length; ++i) {
131 check &= derivatives[i] == ((index + 1 == i) ? 1.0 : 0.0);
132 }
133 return check;
134 }
135
136 /** Check if a derivative represents a specified variable.
137 * @param g derivative to check
138 * @param index index of the variable
139 * @return true if the derivative represents a specified variable
140 */
141 protected boolean isVariable(final Gradient g, final int index) {
142 final double[] derivatives = g.getGradient();
143 boolean check = true;
144 for (int i = 0; i < derivatives.length; ++i) {
145 check &= derivatives[i] == ((index == i) ? 1.0 : 0.0);
146 }
147 return check;
148 }
149
150 /** Compute density and its derivatives.
151 * Using finite differences for the derivatives.
152 * And doing the actual computation only for the derivatives with respect to position (others are set to 0.).
153 * <p>
154 * From a theoretical point of view, this method computes the same values
155 * as {@link Atmosphere#getDensity(FieldAbsoluteDate, FieldVector3D, Frame)} in the
156 * specific case of {@link DerivativeStructure} with respect to state, so
157 * it is less general. However, it is *much* faster in this important case.
158 * </p>
159 * <p>
160 * The derivatives should be computed with respect to position. The input
161 * parameters already take into account the free parameters (6, 7 or 8 depending
162 * on derivation with respect to drag coefficient and lift ratio being considered or not)
163 * and order (always 1). Free parameters at indices 0, 1 and 2 correspond to derivatives
164 * with respect to position. Free parameters at indices 3, 4 and 5 correspond
165 * to derivatives with respect to velocity (these derivatives will remain zero
166 * as the atmospheric density does not depend on velocity). Free parameter
167 * at indexes 6 and 7 (if present) corresponds to derivatives with respect to drag coefficient
168 * and/or lift ratio (one of these or both).
169 * This 2 last derivatives will remain zero as atmospheric density does not depend on them.
170 * </p>
171 * @param date current date
172 * @param frame inertial reference frame for state (both orbit and attitude)
173 * @param position position of spacecraft in inertial frame
174 * @return the density and its derivatives
175 */
176 protected DerivativeStructure getDSDensityWrtStateUsingFiniteDifferences(final AbsoluteDate date,
177 final Frame frame,
178 final FieldVector3D<DerivativeStructure> position) {
179
180 // Retrieve derivation properties for parameter T
181 // It is implied here that T is a DerivativeStructure
182 // With order 1 and 6, 7 or 8 free parameters
183 // This is all checked before in method isStateDerivatives
184 final DSFactory factory = position.getX().getFactory();
185
186 // Build a DerivativeStructure using only derivatives with respect to position
187 final DSFactory factory3 = new DSFactory(3, 1);
188 final FieldVector3D<DerivativeStructure> position3 =
189 new FieldVector3D<>(factory3.variable(0, position.getX().getReal()),
190 factory3.variable(1, position.getY().getReal()),
191 factory3.variable(2, position.getZ().getReal()));
192
193 // Get atmosphere properties in atmosphere own frame
194 final Frame atmFrame = atmosphere.getFrame();
195 final Transform toBody = frame.getTransformTo(atmFrame, date);
196 final FieldVector3D<DerivativeStructure> posBodyDS = toBody.transformPosition(position3);
197 final Vector3D posBody = posBodyDS.toVector3D();
198
199 // Estimate density model by finite differences and composition
200 // Using a delta of 1m
201 final double delta = 1.0;
202 final double x = posBody.getX();
203 final double y = posBody.getY();
204 final double z = posBody.getZ();
205 final double rho0 = atmosphere.getDensity(date, posBody, atmFrame);
206 final double dRhodX = (atmosphere.getDensity(date, new Vector3D(x + delta, y, z), atmFrame) - rho0) / delta;
207 final double dRhodY = (atmosphere.getDensity(date, new Vector3D(x, y + delta, z), atmFrame) - rho0) / delta;
208 final double dRhodZ = (atmosphere.getDensity(date, new Vector3D(x, y, z + delta), atmFrame) - rho0) / delta;
209 final double[] dXdQ = posBodyDS.getX().getAllDerivatives();
210 final double[] dYdQ = posBodyDS.getY().getAllDerivatives();
211 final double[] dZdQ = posBodyDS.getZ().getAllDerivatives();
212
213 // Density with derivatives:
214 // - The value and only the 3 first derivatives (those with respect to spacecraft position) are computed
215 // - Others are set to 0.
216 final int p = factory.getCompiler().getFreeParameters();
217 final double[] rhoAll = new double[p + 1];
218 rhoAll[0] = rho0;
219 for (int i = 1; i < 4; ++i) {
220 rhoAll[i] = dRhodX * dXdQ[i] + dRhodY * dYdQ[i] + dRhodZ * dZdQ[i];
221 }
222
223 return factory.build(rhoAll);
224 }
225
226 /** Compute density and its derivatives.
227 * Using finite differences for the derivatives.
228 * And doing the actual computation only for the derivatives with respect to position (others are set to 0.).
229 * <p>
230 * From a theoretical point of view, this method computes the same values
231 * as {@link Atmosphere#getDensity(FieldAbsoluteDate, FieldVector3D, Frame)} in the
232 * specific case of {@link Gradient} with respect to state, so
233 * it is less general. However, it is *much* faster in this important case.
234 * </p>
235 * <p>
236 * The derivatives should be computed with respect to position. The input
237 * parameters already take into account the free parameters (6, 7 or 8 depending
238 * on derivation with respect to drag coefficient and lift ratio being considered or not)
239 * and order (always 1). Free parameters at indices 0, 1 and 2 correspond to derivatives
240 * with respect to position. Free parameters at indices 3, 4 and 5 correspond
241 * to derivatives with respect to velocity (these derivatives will remain zero
242 * as the atmospheric density does not depend on velocity). Free parameter
243 * at indexes 6 and 7 (if present) corresponds to derivatives with respect to drag coefficient
244 * and/or lift ratio (one of these or both).
245 * This 2 last derivatives will remain zero as atmospheric density does not depend on them.
246 * </p>
247 * @param date current date
248 * @param frame inertial reference frame for state (both orbit and attitude)
249 * @param position position of spacecraft in inertial frame
250 * @return the density and its derivatives
251 */
252 protected Gradient getGradientDensityWrtStateUsingFiniteDifferences(final AbsoluteDate date,
253 final Frame frame,
254 final FieldVector3D<Gradient> position) {
255
256 // Build a Gradient using only derivatives with respect to position
257 final FieldVector3D<Gradient> position3 =
258 new FieldVector3D<>(Gradient.variable(3, 0, position.getX().getReal()),
259 Gradient.variable(3, 1, position.getY().getReal()),
260 Gradient.variable(3, 2, position.getZ().getReal()));
261
262 // Get atmosphere properties in atmosphere own frame
263 final Frame atmFrame = atmosphere.getFrame();
264 final Transform toBody = frame.getTransformTo(atmFrame, date);
265 final FieldVector3D<Gradient> posBodyDS = toBody.transformPosition(position3);
266 final Vector3D posBody = posBodyDS.toVector3D();
267
268 // Estimate density model by finite differences and composition
269 // Using a delta of 1m
270 final double delta = 1.0;
271 final double x = posBody.getX();
272 final double y = posBody.getY();
273 final double z = posBody.getZ();
274 final double rho0 = atmosphere.getDensity(date, posBody, atmFrame);
275 final double dRhodX = (atmosphere.getDensity(date, new Vector3D(x + delta, y, z), atmFrame) - rho0) / delta;
276 final double dRhodY = (atmosphere.getDensity(date, new Vector3D(x, y + delta, z), atmFrame) - rho0) / delta;
277 final double dRhodZ = (atmosphere.getDensity(date, new Vector3D(x, y, z + delta), atmFrame) - rho0) / delta;
278 final double[] dXdQ = posBodyDS.getX().getGradient();
279 final double[] dYdQ = posBodyDS.getY().getGradient();
280 final double[] dZdQ = posBodyDS.getZ().getGradient();
281
282 // Density with derivatives:
283 // - The value and only the 3 first derivatives (those with respect to spacecraft position) are computed
284 // - Others are set to 0.
285 final int p = position.getX().getFreeParameters();
286 final double[] rhoAll = new double[p];
287 for (int i = 0; i < 3; ++i) {
288 rhoAll[i] = dRhodX * dXdQ[i] + dRhodY * dYdQ[i] + dRhodZ * dZdQ[i];
289 }
290
291 return new Gradient(rho0, rhoAll);
292 }
293
294 }