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;
18
19 import java.util.ArrayList;
20 import java.util.List;
21
22 import org.hipparchus.CalculusFieldElement;
23 import org.hipparchus.Field;
24 import org.hipparchus.geometry.euclidean.threed.FieldRotation;
25 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
26 import org.hipparchus.geometry.euclidean.threed.Rotation;
27 import org.hipparchus.geometry.euclidean.threed.Vector3D;
28 import org.hipparchus.util.FastMath;
29 import org.hipparchus.util.Precision;
30 import org.hipparchus.util.SinCos;
31 import org.orekit.errors.OrekitException;
32 import org.orekit.errors.OrekitInternalError;
33 import org.orekit.forces.drag.DragSensitive;
34 import org.orekit.forces.radiation.RadiationSensitive;
35 import org.orekit.frames.Frame;
36 import org.orekit.time.AbsoluteDate;
37 import org.orekit.time.FieldAbsoluteDate;
38 import org.orekit.utils.PVCoordinatesProvider;
39 import org.orekit.utils.ParameterDriver;
40
41 /** Class representing the features of a classical satellite
42 * with a convex body shape and rotating flat solar arrays.
43 * <p>
44 * The body can be either a simple parallelepipedic box aligned with
45 * spacecraft axes or a set of facets defined by their area and normal vector.
46 * This should handle accurately most spacecraft shapes.
47 * </p>
48 * <p>
49 * The solar array rotation with respect to satellite body can be either
50 * the best lighting orientation (i.e. Sun exactly in solar array meridian
51 * plane defined by solar array rotation axis and solar array normal vector)
52 * or a rotation evolving linearly according to a start position and an
53 * angular rate (which can be set to 0 for non-rotating panels, which may
54 * occur in special modes or during contingencies).
55 * </p>
56 * <p>
57 * The lift component of the drag force can be optionally considered. It should
58 * probably only be used for reentry computation, with much denser atmosphere
59 * than in regular orbit propagation. The lift component is computed using a
60 * ratio of molecules that experience specular reflection instead of diffuse
61 * reflection (absorption followed by outgassing at negligible velocity).
62 * Without lift (i.e. when the lift ratio is set to 0), drag force is along
63 * atmosphere relative velocity. With lift (i.e. when the lift ratio is set to any
64 * value between 0 and 1), the drag force depends on both relative velocity direction
65 * and facets normal orientation. For a single panel, if the relative velocity is
66 * head-on (i.e. aligned with the panel normal), the force will be in the same
67 * direction with and without lift, but the magnitude with lift ratio set to 1.0 will
68 * be twice the magnitude with lift ratio set to 0.0 (because atmosphere molecules
69 * bounces backward at same velocity in case of specular reflection).
70 * </p>
71 * <p>
72 * This model does not take cast shadow between body and solar array into account.
73 * </p>
74 *
75 * @author Luc Maisonobe
76 * @author Pascal Parraud
77 */
78 public class BoxAndSolarArraySpacecraft implements RadiationSensitive, DragSensitive {
79
80 /** Parameters scaling factor.
81 * <p>
82 * We use a power of 2 to avoid numeric noise introduction
83 * in the multiplications/divisions sequences.
84 * </p>
85 */
86 private final double SCALE = FastMath.scalb(1.0, -3);
87
88 /** Driver for drag coefficient parameter. */
89 private final ParameterDriver dragParameterDriver;
90
91 /** Driver for lift ratio parameter (may be null is lift is ignored). */
92 private final ParameterDriver liftParameterDriver;
93
94 /** Driver for radiation pressure absorption coefficient parameter. */
95 private final ParameterDriver absorptionParameterDriver;
96
97 /** Driver for radiation pressure reflection coefficient parameter. */
98 private final ParameterDriver reflectionParameterDriver;
99
100 /** Surface vectors for body facets. */
101 private final List<Facet> facets;
102
103 /** Solar array area (m²). */
104 private final double solarArrayArea;
105
106 /** Reference date for linear rotation angle (may be null). */
107 private final AbsoluteDate referenceDate;
108
109 /** Rotation rate for linear rotation angle. */
110 private final double rotationRate;
111
112 /** Solar array reference axis in spacecraft frame (may be null). */
113 private final Vector3D saX;
114
115 /** Solar array third axis in spacecraft frame (may be null). */
116 private final Vector3D saY;
117
118 /** Solar array rotation axis in spacecraft frame. */
119 private final Vector3D saZ;
120
121 /** Sun model. */
122 private final PVCoordinatesProvider sun;
123
124 /** Build a spacecraft model with best lighting of solar array.
125 * <p>
126 * This constructor builds an instance that completely ignores lift
127 * in atmospheric drag (the value of lift coefficient is set to zero,
128 * and there are no {@link ParameterDriver drivers} to change it).
129 * </p>
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 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 */
147 public BoxAndSolarArraySpacecraft(final double xLength, final double yLength,
148 final double zLength,
149 final PVCoordinatesProvider sun, final double solarArrayArea,
150 final Vector3D solarArrayAxis,
151 final double dragCoeff,
152 final double absorptionCoeff,
153 final double reflectionCoeff) {
154 this(simpleBoxFacets(xLength, yLength, zLength), sun, solarArrayArea, solarArrayAxis,
155 dragCoeff, 0.0, false,
156 absorptionCoeff, reflectionCoeff);
157 }
158
159 /** Build a spacecraft model with best lighting of solar array.
160 * <p>
161 * Solar arrays orientation will be such that at each time the Sun direction
162 * will always be in the solar array meridian plane defined by solar array
163 * rotation axis and solar array normal vector.
164 * </p>
165 * @param xLength length of the body along its X axis (m)
166 * @param yLength length of the body along its Y axis (m)
167 * @param zLength length of the body along its Z axis (m)
168 * @param sun sun model
169 * @param solarArrayArea area of the solar array (m²)
170 * @param solarArrayAxis solar array rotation axis in satellite frame
171 * @param dragCoeff drag coefficient (used only for drag)
172 * @param liftRatio lift ratio (proportion between 0 and 1 of atmosphere modecules
173 * that will experience specular reflection when hitting spacecraft instead
174 * of experiencing diffuse reflection, hence producing lift)
175 * @param absorptionCoeff absorption coefficient between 0.0 an 1.0
176 * (used only for radiation pressure)
177 * @param reflectionCoeff specular reflection coefficient between 0.0 an 1.0
178 * (used only for radiation pressure)
179 * @since 9.0
180 */
181 public BoxAndSolarArraySpacecraft(final double xLength, final double yLength,
182 final double zLength,
183 final PVCoordinatesProvider sun, final double solarArrayArea,
184 final Vector3D solarArrayAxis,
185 final double dragCoeff, final double liftRatio,
186 final double absorptionCoeff,
187 final double reflectionCoeff) {
188 this(simpleBoxFacets(xLength, yLength, zLength), sun, solarArrayArea, solarArrayAxis,
189 dragCoeff, liftRatio,
190 absorptionCoeff, reflectionCoeff);
191 }
192
193 /** Build a spacecraft model with best lighting of solar array.
194 * <p>
195 * The spacecraft body is described by an array of surface vectors. Each facet of
196 * the body is described by a vector normal to the facet (pointing outward of the spacecraft)
197 * and whose norm is the surface area in m².
198 * </p>
199 * <p>
200 * Solar arrays orientation will be such that at each time the Sun direction
201 * will always be in the solar array meridian plane defined by solar array
202 * rotation axis and solar array normal vector.
203 * </p>
204 * @param facets body facets (only the facets with strictly positive area will be stored)
205 * @param sun sun model
206 * @param solarArrayArea area of the solar array (m²)
207 * @param solarArrayAxis solar array rotation axis in satellite frame
208 * @param dragCoeff drag coefficient (used only for drag)
209 * @param absorptionCoeff absorption coefficient between 0.0 an 1.0
210 * (used only for radiation pressure)
211 * @param reflectionCoeff specular reflection coefficient between 0.0 an 1.0
212 * (used only for radiation pressure)
213 */
214 public BoxAndSolarArraySpacecraft(final Facet[] facets,
215 final PVCoordinatesProvider sun, final double solarArrayArea,
216 final Vector3D solarArrayAxis,
217 final double dragCoeff,
218 final double absorptionCoeff,
219 final double reflectionCoeff) {
220 this(facets, sun, solarArrayArea, solarArrayAxis,
221 dragCoeff, 0.0, false,
222 absorptionCoeff, reflectionCoeff);
223 }
224
225 /** Build a spacecraft model with best lighting of solar array.
226 * <p>
227 * The spacecraft body is described by an array of surface vectors. Each facet of
228 * the body is described by a vector normal to the facet (pointing outward of the spacecraft)
229 * and whose norm is the surface area in m².
230 * </p>
231 * <p>
232 * Solar arrays orientation will be such that at each time the Sun direction
233 * will always be in the solar array meridian plane defined by solar array
234 * rotation axis and solar array normal vector.
235 * </p>
236 * @param facets body facets (only the facets with strictly positive area will be stored)
237 * @param sun sun model
238 * @param solarArrayArea area of the solar array (m²)
239 * @param solarArrayAxis solar array rotation axis in satellite frame
240 * @param dragCoeff drag coefficient (used only for drag)
241 * @param liftRatio lift ratio (proportion between 0 and 1 of atmosphere modecules
242 * that will experience specular reflection when hitting spacecraft instead
243 * of experiencing diffuse reflection, hence producing lift)
244 * @param absorptionCoeff absorption coefficient between 0.0 an 1.0
245 * (used only for radiation pressure)
246 * @param reflectionCoeff specular reflection coefficient between 0.0 an 1.0
247 * (used only for radiation pressure)
248 * @since 9.0
249 */
250 public BoxAndSolarArraySpacecraft(final Facet[] facets,
251 final PVCoordinatesProvider sun, final double solarArrayArea,
252 final Vector3D solarArrayAxis,
253 final double dragCoeff, final double liftRatio,
254 final double absorptionCoeff,
255 final double reflectionCoeff) {
256 this(facets, sun, solarArrayArea, solarArrayAxis,
257 dragCoeff, liftRatio, true,
258 absorptionCoeff, reflectionCoeff);
259 }
260
261 /** Build a spacecraft model with best lighting of solar array.
262 * <p>
263 * The spacecraft body is described by an array of surface vectors. Each facet of
264 * the body is described by a vector normal to the facet (pointing outward of the spacecraft)
265 * and whose norm is the surface area in m².
266 * </p>
267 * <p>
268 * Solar arrays orientation will be such that at each time the Sun direction
269 * will always be in the solar array meridian plane defined by solar array
270 * rotation axis and solar array normal vector.
271 * </p>
272 * @param facets body facets (only the facets with strictly positive area will be stored)
273 * @param sun sun model
274 * @param solarArrayArea area of the solar array (m²)
275 * @param solarArrayAxis solar array rotation axis in satellite frame
276 * @param dragCoeff drag coefficient (used only for drag)
277 * @param liftRatio lift ratio (proportion between 0 and 1 of atmosphere modecules
278 * that will experience specular reflection when hitting spacecraft instead
279 * of experiencing diffuse reflection, hence producing lift)
280 * @param useLift if true, lift should be used, otherwise it is completely ignored
281 * @param absorptionCoeff absorption coefficient between 0.0 an 1.0
282 * (used only for radiation pressure)
283 * @param reflectionCoeff specular reflection coefficient between 0.0 an 1.0
284 * (used only for radiation pressure)
285 * @since 9.0
286 */
287 private BoxAndSolarArraySpacecraft(final Facet[] facets,
288 final PVCoordinatesProvider sun, final double solarArrayArea,
289 final Vector3D solarArrayAxis,
290 final double dragCoeff, final double liftRatio, final boolean useLift,
291 final double absorptionCoeff,
292 final double reflectionCoeff) {
293
294 // drag
295 dragParameterDriver = buildDragParameterDriver(dragCoeff);
296 liftParameterDriver = useLift ? buildLiftParameterDriver(liftRatio) : null;
297
298 // radiation pressure
299 absorptionParameterDriver = buildAbsorptionParameterDriver(absorptionCoeff);
300 reflectionParameterDriver = buildReflectionParameterDriver(reflectionCoeff);
301
302 this.facets = filter(facets);
303
304 this.sun = sun;
305 this.solarArrayArea = solarArrayArea;
306 this.referenceDate = null;
307 this.rotationRate = 0;
308
309 this.saZ = solarArrayAxis.normalize();
310 this.saY = null;
311 this.saX = null;
312
313 }
314
315 /** Build a spacecraft model with linear rotation of solar array.
316 * <p>
317 * Solar arrays orientation will be a regular rotation from the
318 * reference orientation at reference date and using a constant
319 * rotation rate.
320 * </p>
321 * @param xLength length of the body along its X axis (m)
322 * @param yLength length of the body along its Y axis (m)
323 * @param zLength length of the body along its Z axis (m)
324 * @param sun sun model
325 * @param solarArrayArea area of the solar array (m²)
326 * @param solarArrayAxis solar array rotation axis in satellite frame
327 * @param referenceDate reference date for the solar array rotation
328 * @param referenceNormal direction of the solar array normal at reference date
329 * in spacecraft frame
330 * @param rotationRate rotation rate of the solar array, may be 0 (rad/s)
331 * @param dragCoeff drag coefficient (used only for drag)
332 * @param absorptionCoeff absorption coefficient between 0.0 an 1.0
333 * (used only for radiation pressure)
334 * @param reflectionCoeff specular reflection coefficient between 0.0 an 1.0
335 * (used only for radiation pressure)
336 */
337 public BoxAndSolarArraySpacecraft(final double xLength, final double yLength,
338 final double zLength,
339 final PVCoordinatesProvider sun, final double solarArrayArea,
340 final Vector3D solarArrayAxis,
341 final AbsoluteDate referenceDate,
342 final Vector3D referenceNormal,
343 final double rotationRate,
344 final double dragCoeff,
345 final double absorptionCoeff,
346 final double reflectionCoeff) {
347 this(simpleBoxFacets(xLength, yLength, zLength), sun, solarArrayArea, solarArrayAxis,
348 referenceDate, referenceNormal, rotationRate,
349 dragCoeff, 0.0, false,
350 absorptionCoeff, reflectionCoeff);
351 }
352
353 /** Build a spacecraft model with linear rotation of solar array.
354 * <p>
355 * Solar arrays orientation will be a regular rotation from the
356 * reference orientation at reference date and using a constant
357 * rotation rate.
358 * </p>
359 * @param xLength length of the body along its X axis (m)
360 * @param yLength length of the body along its Y axis (m)
361 * @param zLength length of the body along its Z axis (m)
362 * @param sun sun model
363 * @param solarArrayArea area of the solar array (m²)
364 * @param solarArrayAxis solar array rotation axis in satellite frame
365 * @param referenceDate reference date for the solar array rotation
366 * @param referenceNormal direction of the solar array normal at reference date
367 * in spacecraft frame
368 * @param rotationRate rotation rate of the solar array, may be 0 (rad/s)
369 * @param dragCoeff drag coefficient (used only for drag)
370 * @param liftRatio lift ratio (proportion between 0 and 1 of atmosphere modecules
371 * that will experience specular reflection when hitting spacecraft instead
372 * of experiencing diffuse reflection, hence producing lift)
373 * @param absorptionCoeff absorption coefficient between 0.0 an 1.0
374 * (used only for radiation pressure)
375 * @param reflectionCoeff specular reflection coefficient between 0.0 an 1.0
376 * (used only for radiation pressure)
377 * @since 9.0
378 */
379 public BoxAndSolarArraySpacecraft(final double xLength, final double yLength,
380 final double zLength,
381 final PVCoordinatesProvider sun, final double solarArrayArea,
382 final Vector3D solarArrayAxis,
383 final AbsoluteDate referenceDate,
384 final Vector3D referenceNormal,
385 final double rotationRate,
386 final double dragCoeff, final double liftRatio,
387 final double absorptionCoeff,
388 final double reflectionCoeff) {
389 this(simpleBoxFacets(xLength, yLength, zLength), sun, solarArrayArea, solarArrayAxis,
390 referenceDate, referenceNormal, rotationRate,
391 dragCoeff, liftRatio, true,
392 absorptionCoeff, reflectionCoeff);
393 }
394
395 /** Build a spacecraft model with linear rotation of solar array.
396 * <p>
397 * The spacecraft body is described by an array of surface vectors. Each facet of
398 * the body is described by a vector normal to the facet (pointing outward of the spacecraft)
399 * and whose norm is the surface area in m².
400 * </p>
401 * <p>
402 * Solar arrays orientation will be a regular rotation from the
403 * reference orientation at reference date and using a constant
404 * rotation rate.
405 * </p>
406 * @param facets body facets (only the facets with strictly positive area will be stored)
407 * @param sun sun model
408 * @param solarArrayArea area of the solar array (m²)
409 * @param solarArrayAxis solar array rotation axis in satellite frame
410 * @param referenceDate reference date for the solar array rotation
411 * @param referenceNormal direction of the solar array normal at reference date
412 * in spacecraft frame
413 * @param rotationRate rotation rate of the solar array, may be 0 (rad/s)
414 * @param dragCoeff drag coefficient (used only for drag)
415 * @param absorptionCoeff absorption coefficient between 0.0 an 1.0
416 * (used only for radiation pressure)
417 * @param reflectionCoeff specular reflection coefficient between 0.0 an 1.0
418 * (used only for radiation pressure)
419 */
420 public BoxAndSolarArraySpacecraft(final Facet[] facets,
421 final PVCoordinatesProvider sun, final double solarArrayArea,
422 final Vector3D solarArrayAxis,
423 final AbsoluteDate referenceDate,
424 final Vector3D referenceNormal,
425 final double rotationRate,
426 final double dragCoeff,
427 final double absorptionCoeff,
428 final double reflectionCoeff) {
429 this(facets, sun, solarArrayArea, solarArrayAxis, referenceDate, referenceNormal, rotationRate,
430 dragCoeff, 0.0, false,
431 absorptionCoeff, reflectionCoeff);
432 }
433
434 /** Build a spacecraft model with linear rotation of solar array.
435 * <p>
436 * The spacecraft body is described by an array of surface vectors. Each facet of
437 * the body is described by a vector normal to the facet (pointing outward of the spacecraft)
438 * and whose norm is the surface area in m².
439 * </p>
440 * <p>
441 * Solar arrays orientation will be a regular rotation from the
442 * reference orientation at reference date and using a constant
443 * rotation rate.
444 * </p>
445 * @param facets body facets (only the facets with strictly positive area will be stored)
446 * @param sun sun model
447 * @param solarArrayArea area of the solar array (m²)
448 * @param solarArrayAxis solar array rotation axis in satellite frame
449 * @param referenceDate reference date for the solar array rotation
450 * @param referenceNormal direction of the solar array normal at reference date
451 * in spacecraft frame
452 * @param rotationRate rotation rate of the solar array, may be 0 (rad/s)
453 * @param dragCoeff drag coefficient (used only for drag)
454 * @param liftRatio lift ratio (proportion between 0 and 1 of atmosphere modecules
455 * that will experience specular reflection when hitting spacecraft instead
456 * of experiencing diffuse reflection, hence producing lift)
457 * @param absorptionCoeff absorption coefficient between 0.0 an 1.0
458 * (used only for radiation pressure)
459 * @param reflectionCoeff specular reflection coefficient between 0.0 an 1.0
460 * (used only for radiation pressure)
461 * @since 9.0
462 */
463 public BoxAndSolarArraySpacecraft(final Facet[] facets,
464 final PVCoordinatesProvider sun, final double solarArrayArea,
465 final Vector3D solarArrayAxis,
466 final AbsoluteDate referenceDate,
467 final Vector3D referenceNormal,
468 final double rotationRate,
469 final double dragCoeff, final double liftRatio,
470 final double absorptionCoeff,
471 final double reflectionCoeff) {
472 this(facets, sun, solarArrayArea, solarArrayAxis, referenceDate, referenceNormal, rotationRate,
473 dragCoeff, liftRatio, true,
474 absorptionCoeff, reflectionCoeff);
475 }
476
477 /** Build a spacecraft model with linear rotation of solar array.
478 * <p>
479 * The spacecraft body is described by an array of surface vectors. Each facet of
480 * the body is described by a vector normal to the facet (pointing outward of the spacecraft)
481 * and whose norm is the surface area in m².
482 * </p>
483 * <p>
484 * Solar arrays orientation will be a regular rotation from the
485 * reference orientation at reference date and using a constant
486 * rotation rate.
487 * </p>
488 * @param facets body facets (only the facets with strictly positive area will be stored)
489 * @param sun sun model
490 * @param solarArrayArea area of the solar array (m²)
491 * @param solarArrayAxis solar array rotation axis in satellite frame
492 * @param referenceDate reference date for the solar array rotation
493 * @param referenceNormal direction of the solar array normal at reference date
494 * in spacecraft frame
495 * @param rotationRate rotation rate of the solar array, may be 0 (rad/s)
496 * @param dragCoeff drag coefficient (used only for drag)
497 * @param liftRatio lift ratio (proportion between 0 and 1 of atmosphere modecules
498 * that will experience specular reflection when hitting spacecraft instead
499 * of experiencing diffuse reflection, hence producing lift)
500 * @param useLift if true, lift should be used, otherwise it is completely ignored
501 * @param absorptionCoeff absorption coefficient between 0.0 an 1.0
502 * (used only for radiation pressure)
503 * @param reflectionCoeff specular reflection coefficient between 0.0 an 1.0
504 * (used only for radiation pressure)
505 * @since 9.0
506 */
507 private BoxAndSolarArraySpacecraft(final Facet[] facets,
508 final PVCoordinatesProvider sun, final double solarArrayArea,
509 final Vector3D solarArrayAxis,
510 final AbsoluteDate referenceDate,
511 final Vector3D referenceNormal,
512 final double rotationRate,
513 final double dragCoeff, final double liftRatio, final boolean useLift,
514 final double absorptionCoeff,
515 final double reflectionCoeff) {
516
517 // drag
518 dragParameterDriver = buildDragParameterDriver(dragCoeff);
519 liftParameterDriver = useLift ? buildLiftParameterDriver(liftRatio) : null;
520
521 // radiation pressure
522 absorptionParameterDriver = buildAbsorptionParameterDriver(absorptionCoeff);
523 reflectionParameterDriver = buildReflectionParameterDriver(reflectionCoeff);
524
525 this.facets = filter(facets.clone());
526
527 this.sun = sun;
528 this.solarArrayArea = solarArrayArea;
529 this.referenceDate = referenceDate;
530 this.rotationRate = rotationRate;
531
532 this.saZ = solarArrayAxis.normalize();
533 this.saY = Vector3D.crossProduct(saZ, referenceNormal).normalize();
534 this.saX = Vector3D.crossProduct(saY, saZ);
535
536 }
537
538 /** Build the parameter driver for drag coefficient.
539 * @param coeff drag coefficient
540 * @return parameter driver for drag coefficient
541 * @since 9.0
542 */
543 private ParameterDriver buildDragParameterDriver(final double coeff) {
544 try {
545 return new ParameterDriver(DragSensitive.DRAG_COEFFICIENT,
546 coeff, SCALE, 0.0, Double.POSITIVE_INFINITY);
547 } catch (OrekitException oe) {
548 // this should never happen
549 throw new OrekitInternalError(oe);
550 }
551 }
552
553 /** Build the parameter driver for lift coefficient.
554 * @param coeff lift coefficient
555 * @return parameter driver for lift coefficient
556 * @since 9.0
557 */
558 private ParameterDriver buildLiftParameterDriver(final double coeff) {
559 try {
560 return new ParameterDriver(DragSensitive.LIFT_RATIO, coeff, SCALE, 0.0, 1.0);
561 } catch (OrekitException oe) {
562 // this should never happen
563 throw new OrekitInternalError(oe);
564 }
565 }
566
567 /** Build the parameter driver for absorption coefficient.
568 * @param coeff absorption coefficient
569 * @return parameter driver for absorption coefficient
570 * @since 9.0
571 */
572 private ParameterDriver buildAbsorptionParameterDriver(final double coeff) {
573 try {
574 return new ParameterDriver(RadiationSensitive.ABSORPTION_COEFFICIENT, coeff, SCALE, 0.0, 1.0);
575 } catch (OrekitException oe) {
576 // this should never happen
577 throw new OrekitInternalError(oe);
578 }
579 }
580
581 /** Build the parameter driver for reflection coefficient.
582 * @param coeff absorption coefficient
583 * @return parameter driver for reflection coefficient
584 * @since 9.0
585 */
586 private ParameterDriver buildReflectionParameterDriver(final double coeff) {
587 try {
588 return new ParameterDriver(RadiationSensitive.REFLECTION_COEFFICIENT, coeff, SCALE, 0.0, 1.0);
589 } catch (OrekitException oe) {
590 // this should never happen
591 throw new OrekitInternalError(oe);
592 }
593 }
594
595 /** {@inheritDoc} */
596 @Override
597 public List<ParameterDriver> getDragParametersDrivers() {
598 // Initialize list of drag parameter drivers
599 final List<ParameterDriver> drivers = new ArrayList<>();
600 drivers.add(dragParameterDriver);
601 // Verify if the driver for lift ratio parameter is defined
602 if (liftParameterDriver != null) {
603 drivers.add(liftParameterDriver);
604 }
605 return drivers;
606 }
607
608 /** {@inheritDoc} */
609 @Override
610 public List<ParameterDriver> getRadiationParametersDrivers() {
611 // Initialize list of drag parameter drivers
612 final List<ParameterDriver> drivers = new ArrayList<>();
613 drivers.add(absorptionParameterDriver);
614 drivers.add(reflectionParameterDriver);
615 return drivers;
616 }
617
618 /** Get solar array normal in spacecraft frame.
619 * @param date current date
620 * @param frame inertial reference frame for state (both orbit and attitude)
621 * @param position position of spacecraft in reference frame
622 * @param rotation orientation (attitude) of the spacecraft with respect to reference frame
623 * @return solar array normal in spacecraft frame
624 */
625 public synchronized Vector3D getNormal(final AbsoluteDate date, final Frame frame,
626 final Vector3D position, final Rotation rotation) {
627
628 if (referenceDate != null) {
629 // use a simple rotation at fixed rate
630 final double alpha = rotationRate * date.durationFrom(referenceDate);
631 final SinCos scAlpha = FastMath.sinCos(alpha);
632 return new Vector3D(scAlpha.cos(), saX, scAlpha.sin(), saY);
633 }
634
635 // compute orientation for best lighting
636 final Vector3D sunInert = sun.getPVCoordinates(date, frame).getPosition().subtract(position).normalize();
637 final Vector3D sunSpacecraft = rotation.applyTo(sunInert);
638 final double d = Vector3D.dotProduct(sunSpacecraft, saZ);
639 final double f = 1 - d * d;
640 if (f < Precision.EPSILON) {
641 // extremely rare case: the sun is along solar array rotation axis
642 // (there will not be much output power ...)
643 // we set up an arbitrary normal
644 return saZ.orthogonal();
645 }
646
647 final double s = 1.0 / FastMath.sqrt(f);
648 return new Vector3D(s, sunSpacecraft, -s * d, saZ);
649
650 }
651
652 /** Get solar array normal in spacecraft frame.
653 * @param date current date
654 * @param frame inertial reference frame for state (both orbit and attitude)
655 * @param position position of spacecraft in reference frame
656 * @param rotation orientation (attitude) of the spacecraft with respect to reference frame
657 * @return solar array normal in spacecraft frame
658 * @param <T> type of the field elements
659 */
660 public synchronized <T extends CalculusFieldElement<T>> FieldVector3D<T> getNormal(final FieldAbsoluteDate<T> date,
661 final Frame frame,
662 final FieldVector3D<T> position,
663 final FieldRotation<T> rotation) {
664
665 if (referenceDate != null) {
666 // use a simple rotation at fixed rate
667 final T alpha = date.durationFrom(referenceDate).multiply(rotationRate);
668 return new FieldVector3D<>(alpha.cos(), saX, alpha.sin(), saY);
669 }
670
671 // compute orientation for best lighting
672 final FieldVector3D<T> sunInert = position.subtract(sun.getPVCoordinates(date.toAbsoluteDate(), frame).getPosition()).negate().normalize();
673 final FieldVector3D<T> sunSpacecraft = rotation.applyTo(sunInert);
674 final T d = FieldVector3D.dotProduct(sunSpacecraft, saZ);
675 final T f = d.multiply(d).subtract(1).negate();
676 if (f.getReal() < Precision.EPSILON) {
677 // extremely rare case: the sun is along solar array rotation axis
678 // (there will not be much output power ...)
679 // we set up an arbitrary normal
680 return new FieldVector3D<>(f.getField(), saZ.orthogonal());
681 }
682
683 final T s = f.sqrt().reciprocal();
684 return new FieldVector3D<>(s, sunSpacecraft,
685 s.multiply(d).negate(), new FieldVector3D<>(date.getField(), saZ));
686
687 }
688
689 /** {@inheritDoc} */
690 @Override
691 public Vector3D dragAcceleration(final AbsoluteDate date, final Frame frame, final Vector3D position,
692 final Rotation rotation, final double mass,
693 final double density, final Vector3D relativeVelocity,
694 final double[] parameters) {
695
696 final double dragCoeff = parameters[0];
697 final double liftRatio = liftParameterDriver == null ? 0.0 : parameters[1];
698
699 // relative velocity in spacecraft frame
700 final double vNorm2 = relativeVelocity.getNormSq();
701 final double vNorm = FastMath.sqrt(vNorm2);
702 final Vector3D vDir = rotation.applyTo(relativeVelocity.scalarMultiply(1.0 / vNorm));
703 final double coeff = density * dragCoeff * vNorm2 / (2.0 * mass);
704 final double oMr = 1 - liftRatio;
705
706 // solar array contribution
707 final Vector3D frontNormal = getNormal(date, frame, position, rotation);
708 final double s = coeff * solarArrayArea * Vector3D.dotProduct(frontNormal, vDir);
709 Vector3D acceleration = new Vector3D(oMr * FastMath.abs(s), vDir,
710 liftRatio * s * 2, frontNormal);
711
712 // body facets contribution
713 for (final Facet facet : facets) {
714 final double dot = Vector3D.dotProduct(facet.getNormal(), vDir);
715 if (dot < 0) {
716 // the facet intercepts the incoming flux
717 final double f = coeff * facet.getArea() * dot;
718 acceleration = new Vector3D(1, acceleration,
719 oMr * FastMath.abs(f), vDir,
720 liftRatio * f * 2, facet.getNormal());
721 }
722 }
723
724 // convert back to inertial frame
725 return rotation.applyInverseTo(acceleration);
726
727 }
728
729 /** {@inheritDoc} */
730 @Override
731 public Vector3D radiationPressureAcceleration(final AbsoluteDate date, final Frame frame, final Vector3D position,
732 final Rotation rotation, final double mass, final Vector3D flux,
733 final double[] parameters) {
734
735 if (flux.getNormSq() < Precision.SAFE_MIN) {
736 // null illumination (we are probably in umbra)
737 return Vector3D.ZERO;
738 }
739
740 // radiation flux in spacecraft frame
741 final Vector3D fluxSat = rotation.applyTo(flux);
742
743 // solar array contribution
744 Vector3D normal = getNormal(date, frame, position, rotation);
745 double dot = Vector3D.dotProduct(normal, fluxSat);
746 if (dot > 0) {
747 // the solar array is illuminated backward,
748 // fix signs to compute contribution correctly
749 dot = -dot;
750 normal = normal.negate();
751 }
752 Vector3D force = facetRadiationAcceleration(normal, solarArrayArea, fluxSat, dot, parameters);
753
754 // body facets contribution
755 for (final Facet bodyFacet : facets) {
756 normal = bodyFacet.getNormal();
757 dot = Vector3D.dotProduct(normal, fluxSat);
758 if (dot < 0) {
759 // the facet intercepts the incoming flux
760 force = force.add(facetRadiationAcceleration(normal, bodyFacet.getArea(), fluxSat, dot, parameters));
761 }
762 }
763
764 // convert to inertial frame
765 return rotation.applyInverseTo(new Vector3D(1.0 / mass, force));
766
767 }
768
769 /** {@inheritDoc} */
770 @Override
771 public <T extends CalculusFieldElement<T>> FieldVector3D<T>
772 dragAcceleration(final FieldAbsoluteDate<T> date, final Frame frame,
773 final FieldVector3D<T> position, final FieldRotation<T> rotation,
774 final T mass, final T density, final FieldVector3D<T> relativeVelocity,
775 final T[] parameters) {
776
777 final T dragCoeff = parameters[0];
778 final T liftRatio = liftParameterDriver == null ? dragCoeff.getField().getZero() : parameters[1];
779
780 // relative velocity in spacecraft frame
781 final T vNorm2 = relativeVelocity.getNormSq();
782 final T vNorm = vNorm2.sqrt();
783 final FieldVector3D<T> vDir = rotation.applyTo(relativeVelocity.scalarMultiply(vNorm.reciprocal()));
784 final T coeff = density.multiply(0.5).multiply(dragCoeff).multiply(vNorm2).divide(mass);
785 final T oMr = liftRatio.negate().add(1);
786
787 // solar array facet contribution
788 final FieldVector3D<T> frontNormal = getNormal(date, frame, position, rotation);
789 final T s = coeff.
790 multiply(solarArrayArea).
791 multiply(FieldVector3D.dotProduct(frontNormal, vDir));
792 FieldVector3D<T> acceleration = new FieldVector3D<>(s.abs().multiply(oMr), vDir,
793 s.multiply(liftRatio).multiply(2), frontNormal);
794
795 // body facets contribution
796 final Field<T> field = coeff.getField();
797 for (final Facet facet : facets) {
798 final T dot = FieldVector3D.dotProduct(facet.getNormal(), vDir);
799 if (dot.getReal() < 0) {
800 // the facet intercepts the incoming flux
801 final T f = coeff.multiply(facet.getArea()).multiply(dot);
802 acceleration = new FieldVector3D<>(field.getOne(), acceleration,
803 f.abs().multiply(oMr), vDir,
804 f.multiply(liftRatio).multiply(2), new FieldVector3D<>(field, facet.getNormal()));
805 }
806 }
807
808 // convert back to inertial frame
809 return rotation.applyInverseTo(acceleration);
810
811 }
812
813 /** {@inheritDoc} */
814 @Override
815 public <T extends CalculusFieldElement<T>> FieldVector3D<T>
816 radiationPressureAcceleration(final FieldAbsoluteDate<T> date, final Frame frame,
817 final FieldVector3D<T> position,
818 final FieldRotation<T> rotation, final T mass,
819 final FieldVector3D<T> flux,
820 final T[] parameters) {
821
822 if (flux.getNormSq().getReal() < Precision.SAFE_MIN) {
823 // null illumination (we are probably in umbra)
824 return FieldVector3D.getZero(date.getField());
825 }
826
827 // radiation flux in spacecraft frame
828 final FieldVector3D<T> fluxSat = rotation.applyTo(flux);
829
830 // solar array contribution
831 FieldVector3D<T> normal = getNormal(date, frame, position, rotation);
832 T dot = FieldVector3D.dotProduct(normal, fluxSat);
833 if (dot.getReal() > 0) {
834 // the solar array is illuminated backward,
835 // fix signs to compute contribution correctly
836 dot = dot.negate();
837 normal = normal.negate();
838 }
839 FieldVector3D<T> force = facetRadiationAcceleration(normal, solarArrayArea, fluxSat, dot, parameters);
840
841 // body facets contribution
842 for (final Facet bodyFacet : facets) {
843 normal = new FieldVector3D<>(date.getField(), bodyFacet.getNormal());
844 dot = FieldVector3D.dotProduct(fluxSat, normal);
845 if (dot.getReal() < 0) {
846 // the facet intercepts the incoming flux
847 force = force.add(facetRadiationAcceleration(normal, bodyFacet.getArea(), fluxSat, dot, parameters));
848 }
849 }
850
851 // convert to inertial frame
852 return rotation.applyInverseTo(new FieldVector3D<>(mass.reciprocal(), force));
853
854 }
855
856 /** Compute contribution of one facet to force.
857 * <p>This method implements equation 8-44 from David A. Vallado's
858 * Fundamentals of Astrodynamics and Applications, third edition,
859 * 2007, Microcosm Press.</p>
860 * @param normal facet normal
861 * @param area facet area
862 * @param fluxSat radiation pressure flux in spacecraft frame
863 * @param dot dot product of facet and fluxSat (must be negative)
864 * @param parameters values of the force model parameters
865 * @return contribution of the facet to force in spacecraft frame
866 */
867 private Vector3D facetRadiationAcceleration(final Vector3D normal, final double area, final Vector3D fluxSat,
868 final double dot, final double[] parameters) {
869
870 final double absorptionCoeff = parameters[0];
871 final double specularReflectionCoeff = parameters[1];
872 final double diffuseReflectionCoeff = 1 - (absorptionCoeff + specularReflectionCoeff);
873
874 final double psr = fluxSat.getNorm();
875
876 // Vallado's equation 8-44 uses different parameters which are related to our parameters as:
877 // cos (phi) = -dot / (psr * area)
878 // n = facet / area
879 // s = -fluxSat / psr
880 final double cN = 2 * area * dot * (diffuseReflectionCoeff / 3 - specularReflectionCoeff * dot / psr);
881 final double cS = (area * dot / psr) * (specularReflectionCoeff - 1);
882 return new Vector3D(cN, normal, cS, fluxSat);
883
884 }
885
886 /** Compute contribution of one facet to force.
887 * <p>This method implements equation 8-44 from David A. Vallado's
888 * Fundamentals of Astrodynamics and Applications, third edition,
889 * 2007, Microcosm Press.</p>
890 * @param normal facet normal
891 * @param area facet area
892 * @param fluxSat radiation pressure flux in spacecraft frame
893 * @param dot dot product of facet and fluxSat (must be negative)
894 * @param parameters values of the force model parameters
895 * @param <T> type of the field elements
896 * @return contribution of the facet to force in spacecraft frame
897 */
898 private <T extends CalculusFieldElement<T>> FieldVector3D<T>
899 facetRadiationAcceleration(final FieldVector3D<T> normal, final double area, final FieldVector3D<T> fluxSat,
900 final T dot, final T[] parameters) {
901
902 final T absorptionCoeff = parameters[0];
903 final T specularReflectionCoeff = parameters[1];
904 final T diffuseReflectionCoeff = absorptionCoeff.add(specularReflectionCoeff).negate().add(1);
905
906 final T psr = fluxSat.getNorm();
907
908 // Vallado's equation 8-44 uses different parameters which are related to our parameters as:
909 // cos (phi) = -dot / (psr * area)
910 // n = facet / area
911 // s = -fluxSat / psr
912 final T cN = dot.multiply(-2 * area).multiply(dot.multiply(specularReflectionCoeff).divide(psr).subtract(diffuseReflectionCoeff.divide(3)));
913 final T cS = dot.multiply(area).multiply(specularReflectionCoeff.subtract(1)).divide(psr);
914 return new FieldVector3D<>(cN, normal, cS, fluxSat);
915
916 }
917
918 /** Class representing a single facet of a convex spacecraft body.
919 * <p>Instance of this class are guaranteed to be immutable.</p>
920 * @author Luc Maisonobe
921 */
922 public static class Facet {
923
924 /** Unit Normal vector, pointing outward. */
925 private final Vector3D normal;
926
927 /** Area in m². */
928 private final double area;
929
930 /** Simple constructor.
931 * @param normal vector normal to the facet, pointing outward (will be normalized)
932 * @param area facet area in m²
933 */
934 public Facet(final Vector3D normal, final double area) {
935 this.normal = normal.normalize();
936 this.area = area;
937 }
938
939 /** Get unit normal vector.
940 * @return unit normal vector
941 */
942 public Vector3D getNormal() {
943 return normal;
944 }
945
946 /** Get facet area.
947 * @return facet area
948 */
949 public double getArea() {
950 return area;
951 }
952
953 }
954
955 /** Build the surface vectors for body facets of a simple parallelepipedic box.
956 * @param xLength length of the body along its X axis (m)
957 * @param yLength length of the body along its Y axis (m)
958 * @param zLength length of the body along its Z axis (m)
959 * @return surface vectors array
960 */
961 private static Facet[] simpleBoxFacets(final double xLength, final double yLength, final double zLength) {
962 return new Facet[] {
963 new Facet(Vector3D.MINUS_I, yLength * zLength),
964 new Facet(Vector3D.PLUS_I, yLength * zLength),
965 new Facet(Vector3D.MINUS_J, xLength * zLength),
966 new Facet(Vector3D.PLUS_J, xLength * zLength),
967 new Facet(Vector3D.MINUS_K, xLength * yLength),
968 new Facet(Vector3D.PLUS_K, xLength * yLength)
969 };
970 }
971
972 /** Filter out zero area facets.
973 * @param facets original facets (may include zero area facets)
974 * @return filtered array
975 */
976 private static List<Facet> filter(final Facet[] facets) {
977 final List<Facet> filtered = new ArrayList<Facet>(facets.length);
978 for (Facet facet : facets) {
979 if (facet.getArea() > 0) {
980 filtered.add(facet);
981 }
982 }
983 return filtered;
984 }
985
986 }