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