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.forces;
18
19 import java.util.ArrayList;
20 import java.util.Collections;
21 import java.util.List;
22 import java.util.stream.Collectors;
23
24 import org.hipparchus.CalculusFieldElement;
25 import org.hipparchus.Field;
26 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
27 import org.hipparchus.geometry.euclidean.threed.Vector3D;
28 import org.hipparchus.util.FastMath;
29 import org.hipparchus.util.Precision;
30 import org.orekit.errors.OrekitException;
31 import org.orekit.errors.OrekitInternalError;
32 import org.orekit.forces.drag.DragSensitive;
33 import org.orekit.forces.radiation.RadiationSensitive;
34 import org.orekit.propagation.FieldSpacecraftState;
35 import org.orekit.propagation.SpacecraftState;
36 import org.orekit.utils.ExtendedPositionProvider;
37 import org.orekit.utils.ParameterDriver;
38
39 /** Class representing the features of a classical satellite with a convex body shape.
40 * <p>
41 * The body can be either a simple parallelepipedic box aligned with
42 * spacecraft axes or a set of panels defined by their area and normal vector.
43 * Some panels may be moving to model solar arrays (or antennas that could
44 * point anywhere). This should handle accurately most spacecraft shapes. This model
45 * does not take cast shadows into account.
46 * </p>
47 * <p>
48 * The lift component of the drag force can be optionally considered. It should
49 * probably only be used for reentry computation, with much denser atmosphere
50 * than in regular orbit propagation. The lift component is computed using a
51 * ratio of molecules that experience specular reflection instead of diffuse
52 * reflection (absorption followed by outgassing at negligible velocity).
53 * Without lift (i.e. when the lift ratio is set to 0), drag force is along
54 * atmosphere relative velocity. With lift (i.e. when the lift ratio is set to any
55 * value between 0 and 1), the drag force depends on both relative velocity direction
56 * and panels normal orientation. For a single panel, if the relative velocity is
57 * head-on (i.e. aligned with the panel normal), the force will be in the same
58 * direction with and without lift, but the magnitude with lift ratio set to 1.0 will
59 * be twice the magnitude with lift ratio set to 0.0 (because atmosphere molecules
60 * bounces backward at same velocity in case of specular reflection).
61 * </p>
62 * <p>
63 * Each {@link Panel panel} has its own set of radiation and drag coefficients. In
64 * orbit determination context, it would not be possible to estimate each panel
65 * individually, therefore {@link #getDragParametersDrivers()} returns a single
66 * {@link ParameterDriver parameter driver} representing a {@link DragSensitive#GLOBAL_DRAG_FACTOR
67 * global drag multiplicative factor} that applies to all panels drag coefficients
68 * and the {@link #getRadiationParametersDrivers()} returns a single
69 * {@link ParameterDriver parameter driver} representing a
70 * {@link RadiationSensitive#GLOBAL_RADIATION_FACTOR global radiation multiplicative factor}
71 * that applies to all panels radiation coefficients.
72 * </p>
73 *
74 * @author Luc Maisonobe
75 * @author Pascal Parraud
76 */
77 public class BoxAndSolarArraySpacecraft implements RadiationSensitive, DragSensitive {
78
79 /** Parameters scaling factor.
80 * <p>
81 * We use a power of 2 to avoid numeric noise introduction
82 * in the multiplications/divisions sequences.
83 * </p>
84 */
85 private final double SCALE = FastMath.scalb(1.0, -3);
86
87 /** Driver for drag multiplicative factor parameter. */
88 private final ParameterDriver dragFactorParameterDriver;
89
90 /** Driver for radiation pressure multiplicative factor parameter. */
91 private final ParameterDriver radiationFactorParameterDriver;
92
93 /** Panels composing the spacecraft. */
94 private final List<Panel> panels;
95
96 /** Build a spacecraft model.
97 * @param panels panels composing the body, solar arrays and antennas
98 * (only the panels with strictly positive area will be stored)
99 * @since 12.0
100 */
101 public BoxAndSolarArraySpacecraft(final List<Panel> panels) {
102
103 try {
104 dragFactorParameterDriver = new ParameterDriver(DragSensitive.GLOBAL_DRAG_FACTOR,
105 1.0, SCALE, 0.0, Double.POSITIVE_INFINITY);
106 radiationFactorParameterDriver = new ParameterDriver(RadiationSensitive.GLOBAL_RADIATION_FACTOR,
107 1.0, SCALE, 0.0, Double.POSITIVE_INFINITY);
108 } catch (OrekitException oe) {
109 // this should never happen
110 throw new OrekitInternalError(oe);
111 }
112
113 // remove spurious panels
114 this.panels = panels.stream().filter(p -> p.getArea() > 0).collect(Collectors.toList());
115
116 }
117
118 /** Build a spacecraft model with best lighting of solar array.
119 * <p>
120 * Solar arrays orientation will be such that at each time the Sun direction
121 * will always be in the solar array meridian plane defined by solar array
122 * rotation axis and solar array normal vector.
123 * </p>
124 * @param xLength length of the body along its X axis (m)
125 * @param yLength length of the body along its Y axis (m)
126 * @param zLength length of the body along its Z axis (m)
127 * @param sun sun model
128 * @param solarArrayArea area of the solar array (m²)
129 * @param solarArrayAxis solar array rotation axis in satellite frame
130 * @param dragCoeff drag coefficient (used only for drag)
131 * @param liftRatio lift ratio (proportion between 0 and 1 of atmosphere modecules
132 * that will experience specular reflection when hitting spacecraft instead
133 * of experiencing diffuse reflection, hence producing lift)
134 * @param absorptionCoeff absorption coefficient between 0.0 an 1.0
135 * (used only for radiation pressure)
136 * @param reflectionCoeff specular reflection coefficient between 0.0 an 1.0
137 * (used only for radiation pressure)
138 * @since 12.0
139 */
140 public BoxAndSolarArraySpacecraft(final double xLength, final double yLength, final double zLength,
141 final ExtendedPositionProvider sun,
142 final double solarArrayArea, final Vector3D solarArrayAxis,
143 final double dragCoeff, final double liftRatio,
144 final double absorptionCoeff, final double reflectionCoeff) {
145 this(buildPanels(xLength, yLength, zLength,
146 sun, solarArrayArea, solarArrayAxis,
147 dragCoeff, liftRatio, absorptionCoeff, reflectionCoeff));
148 }
149
150 /** Get the panels composing the body.
151 * @return unmodifiable view of the panels composing the body
152 * @since 12.0
153 */
154 public List<Panel> getPanels() {
155 return Collections.unmodifiableList(panels);
156 }
157
158 /** {@inheritDoc} */
159 @Override
160 public List<ParameterDriver> getDragParametersDrivers() {
161 return Collections.singletonList(dragFactorParameterDriver);
162 }
163
164 /** {@inheritDoc} */
165 @Override
166 public List<ParameterDriver> getRadiationParametersDrivers() {
167 return Collections.singletonList(radiationFactorParameterDriver);
168 }
169
170 /** {@inheritDoc} */
171 @Override
172 public Vector3D dragAcceleration(final SpacecraftState state,
173 final double density, final Vector3D relativeVelocity,
174 final double[] parameters) {
175
176 final double dragFactor = parameters[0];
177
178 // relative velocity in spacecraft frame
179 final double vNorm2 = relativeVelocity.getNormSq();
180 final double vNorm = FastMath.sqrt(vNorm2);
181 final Vector3D vDir = state.getAttitude().getRotation().applyTo(relativeVelocity.scalarMultiply(1.0 / vNorm));
182 final double coeff = density * dragFactor * vNorm2 / (2.0 * state.getMass());
183
184 // panels contribution
185 Vector3D acceleration = Vector3D.ZERO;
186 for (final Panel panel : panels) {
187 Vector3D normal = panel.getNormal(state);
188 double dot = Vector3D.dotProduct(normal, vDir);
189 if (panel.isDoubleSided() && dot > 0) {
190 // the flux comes from the back side
191 normal = normal.negate();
192 dot = -dot;
193 }
194 if (dot < 0) {
195 // the panel intercepts the incoming flux
196 final double f = coeff * panel.getDrag() * panel.getArea() * dot;
197 final double liftRatio = panel.getLiftRatio();
198 acceleration = new Vector3D(1, acceleration,
199 (1 - liftRatio) * FastMath.abs(f), vDir,
200 liftRatio * f * 2, normal);
201 }
202 }
203
204 // convert back to inertial frame
205 return state.getAttitude().getRotation().applyInverseTo(acceleration);
206
207 }
208
209 /** {@inheritDoc} */
210 @Override
211 public <T extends CalculusFieldElement<T>> FieldVector3D<T>
212 dragAcceleration(final FieldSpacecraftState<T> state,
213 final T density, final FieldVector3D<T> relativeVelocity,
214 final T[] parameters) {
215
216 final Field<T> field = state.getDate().getField();
217 final T dragFactor = parameters[0];
218
219 // relative velocity in spacecraft frame
220 final T vNorm2 = relativeVelocity.getNormSq();
221 final T vNorm = FastMath.sqrt(vNorm2);
222 final FieldVector3D<T> vDir = state.getAttitude().getRotation().applyTo(relativeVelocity.scalarMultiply(vNorm.reciprocal()));
223 final T coeff = density.multiply(dragFactor).multiply(vNorm2).divide(state.getMass().multiply(2.0));
224
225 // panels contribution
226 FieldVector3D<T> acceleration = FieldVector3D.getZero(field);
227 for (final Panel panel : panels) {
228 FieldVector3D<T> normal = panel.getNormal(state);
229 T dot = FieldVector3D.dotProduct(normal, vDir);
230 if (panel.isDoubleSided() && dot.getReal() > 0) {
231 // the flux comes from the back side
232 normal = normal.negate();
233 dot = dot.negate();
234 }
235 if (panel.isDoubleSided() || dot.getReal() < 0) {
236 // the panel intercepts the incoming flux
237 final T f = coeff.multiply(panel.getDrag() * panel.getArea()).multiply(dot);
238 final double liftRatio = panel.getLiftRatio();
239 acceleration = new FieldVector3D<>(field.getOne(), acceleration,
240 FastMath.abs(f).multiply(1 - liftRatio), vDir,
241 f.multiply(2 * liftRatio), normal);
242 }
243 }
244
245 // convert back to inertial frame
246 return state.getAttitude().getRotation().applyInverseTo(acceleration);
247
248 }
249
250 /** {@inheritDoc} */
251 @Override
252 public Vector3D radiationPressureAcceleration(final SpacecraftState state,
253 final Vector3D flux,
254 final double[] parameters) {
255
256 if (flux.getNormSq() < Precision.SAFE_MIN) {
257 // null illumination (we are probably in umbra)
258 return Vector3D.ZERO;
259 }
260
261 // radiation flux in spacecraft frame
262 final double radiationFactor = parameters[0];
263 final Vector3D fluxSat = state.getAttitude().getRotation().applyTo(flux).
264 scalarMultiply(radiationFactor);
265
266 // panels contribution
267 Vector3D force = Vector3D.ZERO;
268 for (final Panel panel : panels) {
269 Vector3D normal = panel.getNormal(state);
270 double dot = Vector3D.dotProduct(normal, fluxSat);
271 if (panel.isDoubleSided() && dot > 0) {
272 // the flux comes from the back side
273 normal = normal.negate();
274 dot = -dot;
275 }
276 if (dot < 0) {
277 // the panel intercepts the incoming flux
278
279 final double absorptionCoeff = panel.getAbsorption();
280 final double specularReflectionCoeff = panel.getReflection();
281 final double diffuseReflectionCoeff = 1 - (absorptionCoeff + specularReflectionCoeff);
282 final double psr = fluxSat.getNorm();
283
284 // Vallado's equation 8-44 uses different parameters which are related to our parameters as:
285 // cos (phi) = -dot / (psr * area)
286 // n = panel / area
287 // s = -fluxSat / psr
288 final double cN = 2 * panel.getArea() * dot * (diffuseReflectionCoeff / 3 - specularReflectionCoeff * dot / psr);
289 final double cS = (panel.getArea() * dot / psr) * (specularReflectionCoeff - 1);
290 force = new Vector3D(1, force, cN, normal, cS, fluxSat);
291
292 }
293 }
294
295 // convert to inertial frame
296 return state.getAttitude().getRotation().applyInverseTo(new Vector3D(1.0 / state.getMass(), force));
297
298 }
299
300 /** {@inheritDoc}
301 * <p>This method implements equation 8-44 from David A. Vallado's
302 * Fundamentals of Astrodynamics and Applications, third edition,
303 * 2007, Microcosm Press.</p>
304 */
305 @Override
306 public <T extends CalculusFieldElement<T>> FieldVector3D<T>
307 radiationPressureAcceleration(final FieldSpacecraftState<T> state,
308 final FieldVector3D<T> flux,
309 final T[] parameters) {
310
311 final Field<T> field = state.getDate().getField();
312 if (flux.getNormSq().getReal() < Precision.SAFE_MIN) {
313 // null illumination (we are probably in umbra)
314 return FieldVector3D.getZero(field);
315 }
316
317 // radiation flux in spacecraft frame
318 final T radiationFactor = parameters[0];
319 final FieldVector3D<T> fluxSat = state.getAttitude().getRotation().applyTo(flux).
320 scalarMultiply(radiationFactor);
321
322 // panels contribution
323 FieldVector3D<T> force = FieldVector3D.getZero(field);
324 for (final Panel panel : panels) {
325 FieldVector3D<T> normal = panel.getNormal(state);
326 T dot = FieldVector3D.dotProduct(normal, fluxSat);
327 if (panel.isDoubleSided() && dot.getReal() > 0) {
328 // the flux comes from the back side
329 normal = normal.negate();
330 dot = dot.negate();
331 }
332 if (dot.getReal() < 0) {
333 // the panel intercepts the incoming flux
334
335 final double absorptionCoeff = panel.getAbsorption();
336 final double specularReflectionCoeff = panel.getReflection();
337 final double diffuseReflectionCoeff = 1 - (absorptionCoeff + specularReflectionCoeff);
338 final T psr = fluxSat.getNorm();
339
340 // Vallado's equation 8-44 uses different parameters which are related to our parameters as:
341 // cos (phi) = -dot / (psr * area)
342 // n = panel / area
343 // s = -fluxSat / psr
344 final T cN = dot.multiply(-2 * panel.getArea()).multiply(dot.multiply(specularReflectionCoeff).divide(psr).subtract(diffuseReflectionCoeff / 3));
345 final T cS = dot.multiply(panel.getArea()).multiply(specularReflectionCoeff - 1).divide(psr);
346 force = new FieldVector3D<>(field.getOne(), force, cN, normal, cS, fluxSat);
347 }
348 }
349
350 // convert to inertial frame
351 return state.getAttitude().getRotation().applyInverseTo(new FieldVector3D<>(state.getMass().reciprocal(), force));
352
353 }
354
355 /** Build the panels of a simple parallelepipedic box.
356 * @param xLength length of the body along its X axis (m)
357 * @param yLength length of the body along its Y axis (m)
358 * @param zLength length of the body along its Z axis (m)
359 * @param drag drag coefficient
360 * @param liftRatio drag lift ratio (proportion between 0 and 1 of atmosphere modecules
361 * that will experience specular reflection when hitting spacecraft instead
362 * of experiencing diffuse reflection, hence producing lift)
363 * @param absorption radiation pressure absorption coefficient (between 0 and 1)
364 * @param reflection radiation pressure specular reflection coefficient (between 0 and 1)
365 * @return surface vectors array
366 * @since 12.0
367 */
368 public static List<Panel> buildBox(final double xLength, final double yLength, final double zLength,
369 final double drag, final double liftRatio,
370 final double absorption, final double reflection) {
371
372 final List<Panel> panels = new ArrayList<>(6);
373
374 // spacecraft body, composed of single-sided panels
375 panels.add(new FixedPanel(Vector3D.MINUS_I, yLength * zLength, false, drag, liftRatio, absorption, reflection));
376 panels.add(new FixedPanel(Vector3D.PLUS_I, yLength * zLength, false, drag, liftRatio, absorption, reflection));
377 panels.add(new FixedPanel(Vector3D.MINUS_J, xLength * zLength, false, drag, liftRatio, absorption, reflection));
378 panels.add(new FixedPanel(Vector3D.PLUS_J, xLength * zLength, false, drag, liftRatio, absorption, reflection));
379 panels.add(new FixedPanel(Vector3D.MINUS_K, xLength * yLength, false, drag, liftRatio, absorption, reflection));
380 panels.add(new FixedPanel(Vector3D.PLUS_K, xLength * yLength, false, drag, liftRatio, absorption, reflection));
381
382 return panels;
383
384 }
385
386 /** Build the panels of a simple parallelepiped box plus one solar array panel.
387 * @param xLength length of the body along its X axis (m)
388 * @param yLength length of the body along its Y axis (m)
389 * @param zLength length of the body along its Z axis (m)
390 * @param sun sun model
391 * @param solarArrayArea area of the solar array (m²)
392 * @param solarArrayAxis solar array rotation axis in satellite frame
393 * @param drag drag coefficient
394 * @param liftRatio drag lift ratio (proportion between 0 and 1 of atmosphere modecules
395 * that will experience specular reflection when hitting spacecraft instead
396 * of experiencing diffuse reflection, hence producing lift)
397 * @param absorption radiation pressure absorption coefficient (between 0 and 1)
398 * @param reflection radiation pressure specular reflection coefficient (between 0 and 1)
399 * @return surface vectors array
400 * @since 12.0
401 */
402 public static List<Panel> buildPanels(final double xLength, final double yLength, final double zLength,
403 final ExtendedPositionProvider sun,
404 final double solarArrayArea, final Vector3D solarArrayAxis,
405 final double drag, final double liftRatio,
406 final double absorption, final double reflection) {
407
408 // spacecraft body
409 final List<Panel> panels = buildBox(xLength, yLength, zLength, drag, liftRatio, absorption, reflection);
410
411 // solar array
412 panels.add(new PointingPanel(solarArrayAxis, sun, solarArrayArea, drag, liftRatio, absorption, reflection));
413
414 return panels;
415
416 }
417
418 }