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