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  import java.util.stream.Stream;
24  
25  import org.hamcrest.MatcherAssert;
26  import org.hipparchus.CalculusFieldElement;
27  import org.hipparchus.Field;
28  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
29  import org.hipparchus.geometry.euclidean.threed.Rotation;
30  import org.hipparchus.geometry.euclidean.threed.Vector3D;
31  import org.hipparchus.util.Binary64;
32  import org.hipparchus.util.Binary64Field;
33  import org.hipparchus.util.FastMath;
34  import org.hipparchus.util.MathArrays;
35  import org.hipparchus.util.Precision;
36  import org.junit.jupiter.api.Assertions;
37  import org.junit.jupiter.api.BeforeEach;
38  import org.junit.jupiter.api.Test;
39  import org.orekit.OrekitMatchers;
40  import org.orekit.Utils;
41  import org.orekit.attitudes.Attitude;
42  import org.orekit.attitudes.LofOffset;
43  import org.orekit.bodies.CelestialBody;
44  import org.orekit.bodies.CelestialBodyFactory;
45  import org.orekit.errors.OrekitException;
46  import org.orekit.forces.drag.DragSensitive;
47  import org.orekit.forces.radiation.RadiationSensitive;
48  import org.orekit.frames.Frame;
49  import org.orekit.frames.FramesFactory;
50  import org.orekit.frames.LOFType;
51  import org.orekit.orbits.CartesianOrbit;
52  import org.orekit.orbits.CircularOrbit;
53  import org.orekit.orbits.Orbit;
54  import org.orekit.orbits.PositionAngleType;
55  import org.orekit.propagation.FieldSpacecraftState;
56  import org.orekit.propagation.Propagator;
57  import org.orekit.propagation.SpacecraftState;
58  import org.orekit.propagation.analytical.EcksteinHechlerPropagator;
59  import org.orekit.time.AbsoluteDate;
60  import org.orekit.time.DateComponents;
61  import org.orekit.time.TimeComponents;
62  import org.orekit.time.TimeScalesFactory;
63  import org.orekit.utils.Constants;
64  import org.orekit.utils.ParameterDriver;
65  import org.orekit.utils.TimeStampedAngularCoordinates;
66  import org.orekit.utils.TimeStampedPVCoordinates;
67  
68  public class BoxAndSolarArraySpacecraftTest {
69  
70      @Test
71      void testCoefficients() {
72          final double drag = 1.1;
73          final double lift = 2.2;
74          final double abso = 3.3;
75          final double refl = 4.4;
76  
77          final double tol = 1.0e-15;
78  
79          // build box
80          List<Panel> box = BoxAndSolarArraySpacecraft.buildBox(1.5, 3.5, 2.5,
81              drag, lift, abso, refl);
82  
83          // check
84          for (Panel panel : box) {
85              Assertions.assertEquals(drag, panel.getDrag(), tol);
86              Assertions.assertEquals(lift, panel.getLiftRatio(), tol);
87              Assertions.assertEquals(refl, panel.getReflection(), tol);
88              Assertions.assertEquals(abso, panel.getAbsorption(), tol);
89          }
90  
91          // build panels
92          List<Panel> panels = BoxAndSolarArraySpacecraft.buildPanels(1.5, 3.5, 2.5,
93              CelestialBodyFactory.getSun(), 0.4, Vector3D.PLUS_J,
94              drag, lift, abso, refl);
95  
96          // check
97          for (Panel panel : panels) {
98              Assertions.assertEquals(drag, panel.getDrag(), tol);
99              Assertions.assertEquals(lift, panel.getLiftRatio(), tol);
100             Assertions.assertEquals(refl, panel.getReflection(), tol);
101             Assertions.assertEquals(abso, panel.getAbsorption(), tol);
102         }
103 
104     }
105 
106     @Test
107     void testParametersDrivers() {
108 
109         CelestialBody sun = CelestialBodyFactory.getSun();
110         List<Panel> cube = new ArrayList<>();
111         cube.add(new FixedPanel(Vector3D.MINUS_I, 3.0, false, 2.0, 0.0, 0.8, 0.1));
112         cube.add(new FixedPanel(Vector3D.PLUS_I,  3.0, false, 2.0, 0.0, 0.8, 0.1));
113         cube.add(new FixedPanel(Vector3D.MINUS_J, 3.0, false, 2.0, 0.0, 0.8, 0.1));
114         cube.add(new FixedPanel(Vector3D.PLUS_J,  3.0, false, 2.0, 0.0, 0.8, 0.1));
115         cube.add(new FixedPanel(Vector3D.MINUS_K, 3.0, false, 2.0, 0.0, 0.8, 0.1));
116         cube.add(new FixedPanel(Vector3D.PLUS_K,  3.0, false, 2.0, 0.0, 0.8, 0.1));
117         List<Panel> boxNoLift = BoxAndSolarArraySpacecraft.buildBox(1.5, 3.5, 2.5, 2.0, 0.0, 0.8, 0.1);
118         List<Panel> boxLift   = BoxAndSolarArraySpacecraft.buildBox(1.5, 3.5, 2.5, 2.0, 0.4, 0.8, 0.1);
119 
120         BoxAndSolarArraySpacecraft s1 =
121                         new BoxAndSolarArraySpacecraft(1.5, 3.5, 2.5, sun, 20.0, Vector3D.PLUS_J, 2.0, 0.0, 0.8, 0.1);
122         Assertions.assertEquals(1, s1.getDragParametersDrivers().size());
123         Assertions.assertEquals(DragSensitive.GLOBAL_DRAG_FACTOR, s1.getDragParametersDrivers().get(0).getName());
124         Assertions.assertEquals(1.0, s1.getDragParametersDrivers().get(0).getValue(), 1.0e-15);
125         Assertions.assertEquals(1, s1.getRadiationParametersDrivers().size());
126         Assertions.assertEquals(RadiationSensitive.GLOBAL_RADIATION_FACTOR,
127                                 s1.getRadiationParametersDrivers().get(0).getName());
128         Assertions.assertEquals(1, s1.getRadiationParametersDrivers().get(0).getValue(), 1.0e-15);
129 
130         BoxAndSolarArraySpacecraft s2 =
131                         new BoxAndSolarArraySpacecraft(1.5, 3.5, 2.5, sun, 20.0, Vector3D.PLUS_J, 2.0, 0.4, 0.8, 0.1);
132         Assertions.assertEquals(1, s2.getDragParametersDrivers().size());
133         Assertions.assertEquals(DragSensitive.GLOBAL_DRAG_FACTOR, s2.getDragParametersDrivers().get(0).getName());
134         Assertions.assertEquals(1.0, s2.getDragParametersDrivers().get(0).getValue(), 1.0e-15);
135         Assertions.assertEquals(1, s2.getRadiationParametersDrivers().size());
136         Assertions.assertEquals(RadiationSensitive.GLOBAL_RADIATION_FACTOR,
137                                 s2.getRadiationParametersDrivers().get(0).getName());
138         Assertions.assertEquals(1.0, s2.getRadiationParametersDrivers().get(0).getValue(), 1.0e-15);
139 
140         PointingPanel pointingNoLift = new PointingPanel(Vector3D.PLUS_J, sun, 20.0, 2.0, 0.0, 0.8, 0.1);
141         BoxAndSolarArraySpacecraft s3 =
142                         new BoxAndSolarArraySpacecraft(Stream.concat(cube.stream(), Stream.of(pointingNoLift)).
143                                                        collect(Collectors.toList()));
144         Assertions.assertEquals(1, s3.getDragParametersDrivers().size());
145         Assertions.assertEquals(DragSensitive.GLOBAL_DRAG_FACTOR, s3.getDragParametersDrivers().get(0).getName());
146         Assertions.assertEquals(1.0, s3.getDragParametersDrivers().get(0).getValue(), 1.0e-15);
147         Assertions.assertEquals(1, s3.getRadiationParametersDrivers().size());
148         Assertions.assertEquals(RadiationSensitive.GLOBAL_RADIATION_FACTOR,
149                                 s3.getRadiationParametersDrivers().get(0).getName());
150         Assertions.assertEquals(1.0, s3.getRadiationParametersDrivers().get(0).getValue(), 1.0e-15);
151 
152         PointingPanel pointingLift = new PointingPanel(Vector3D.PLUS_J, sun, 20.0, 2.0, 0.4, 0.8, 0.1);
153         BoxAndSolarArraySpacecraft s4 =
154                         new BoxAndSolarArraySpacecraft(Stream.concat(cube.stream(), Stream.of(pointingLift)).
155                                                        collect(Collectors.toList()));
156         Assertions.assertEquals(1, s4.getDragParametersDrivers().size());
157         Assertions.assertEquals(DragSensitive.GLOBAL_DRAG_FACTOR, s4.getDragParametersDrivers().get(0).getName());
158         Assertions.assertEquals(1.0, s4.getDragParametersDrivers().get(0).getValue(), 1.0e-15);
159         Assertions.assertEquals(1, s4.getRadiationParametersDrivers().size());
160         Assertions.assertEquals(RadiationSensitive.GLOBAL_RADIATION_FACTOR,
161                                 s4.getRadiationParametersDrivers().get(0).getName());
162         Assertions.assertEquals(1.0, s4.getRadiationParametersDrivers().get(0).getValue(), 1.0e-15);
163 
164         SlewingPanel slewingNoLift = new SlewingPanel(Vector3D.PLUS_J, 7.292e-5, AbsoluteDate.J2000_EPOCH,
165                                                       Vector3D.PLUS_I, 20.0, 2.0, 0.0, 0.8, 0.1);
166         BoxAndSolarArraySpacecraft s5 = new BoxAndSolarArraySpacecraft(Stream.concat(boxNoLift.stream(), Stream.of(slewingNoLift)).
167                                                                        collect(Collectors.toList()));
168         Assertions.assertEquals(1, s5.getDragParametersDrivers().size());
169         Assertions.assertEquals(DragSensitive.GLOBAL_DRAG_FACTOR, s5.getDragParametersDrivers().get(0).getName());
170         Assertions.assertEquals(1.0, s5.getDragParametersDrivers().get(0).getValue(), 1.0e-15);
171         Assertions.assertEquals(1, s5.getRadiationParametersDrivers().size());
172         Assertions.assertEquals(RadiationSensitive.GLOBAL_RADIATION_FACTOR,
173                 s5.getRadiationParametersDrivers().get(0).getName());
174         Assertions.assertEquals(1.0, s5.getRadiationParametersDrivers().get(0).getValue(), 1.0e-15);
175 
176         SlewingPanel slewingLift = new SlewingPanel(Vector3D.PLUS_J, 7.292e-5, AbsoluteDate.J2000_EPOCH,
177                                                     Vector3D.PLUS_I, 20.0, 2.0, 0.4, 0.8, 0.1);
178         BoxAndSolarArraySpacecraft s6 = new BoxAndSolarArraySpacecraft(Stream.concat(boxLift.stream(), Stream.of(slewingLift)).
179                                                                        collect(Collectors.toList()));
180         Assertions.assertEquals(1, s6.getDragParametersDrivers().size());
181         Assertions.assertEquals(DragSensitive.GLOBAL_DRAG_FACTOR, s6.getDragParametersDrivers().get(0).getName());
182         Assertions.assertEquals(1.0, s6.getDragParametersDrivers().get(0).getValue(), 1.0e-15);
183         Assertions.assertEquals(1, s6.getRadiationParametersDrivers().size());
184         Assertions.assertEquals(RadiationSensitive.GLOBAL_RADIATION_FACTOR,
185                 s6.getRadiationParametersDrivers().get(0).getName());
186         Assertions.assertEquals(1.0, s6.getRadiationParametersDrivers().get(0).getValue(), 1.0e-15);
187 
188         BoxAndSolarArraySpacecraft s7 =
189                         new BoxAndSolarArraySpacecraft(Stream.concat(cube.stream(), Stream.of(slewingNoLift)).
190                                                        collect(Collectors.toList()));
191         Assertions.assertEquals(1, s7.getDragParametersDrivers().size());
192         Assertions.assertEquals(DragSensitive.GLOBAL_DRAG_FACTOR, s7.getDragParametersDrivers().get(0).getName());
193         Assertions.assertEquals(1.0, s7.getDragParametersDrivers().get(0).getValue(), 1.0e-15);
194         Assertions.assertEquals(1, s7.getRadiationParametersDrivers().size());
195         Assertions.assertEquals(RadiationSensitive.GLOBAL_RADIATION_FACTOR,
196                                 s7.getRadiationParametersDrivers().get(0).getName());
197         Assertions.assertEquals(1.0, s7.getRadiationParametersDrivers().get(0).getValue(), 1.0e-15);
198 
199         BoxAndSolarArraySpacecraft s8 =
200                         new BoxAndSolarArraySpacecraft(Stream.concat(cube.stream(), Stream.of(slewingLift)).
201                                                        collect(Collectors.toList()));
202         Assertions.assertEquals(1, s8.getDragParametersDrivers().size());
203         Assertions.assertEquals(DragSensitive.GLOBAL_DRAG_FACTOR, s8.getDragParametersDrivers().get(0).getName());
204         Assertions.assertEquals(1.0, s8.getDragParametersDrivers().get(0).getValue(), 1.0e-15);
205         Assertions.assertEquals(1, s8.getRadiationParametersDrivers().size());
206         Assertions.assertEquals(RadiationSensitive.GLOBAL_RADIATION_FACTOR,
207                 s8.getRadiationParametersDrivers().get(0).getName());
208         Assertions.assertEquals(1.0, s8.getRadiationParametersDrivers().get(0).getValue(), 1.0e-15);
209 
210     }
211 
212     @Test
213     void testNoLiftWithoutReflection() {
214 
215         AbsoluteDate initialDate = propagator.getInitialState().getDate();
216         CelestialBody sun = CelestialBodyFactory.getSun();
217         BoxAndSolarArraySpacecraft s =
218             new BoxAndSolarArraySpacecraft(1.5, 3.5, 2.5, sun, 20.0, Vector3D.PLUS_J, 1.0, 0.0, 1.0, 0.0);
219 
220         Vector3D earthRot = new Vector3D(0.0, 0.0, 7.292115e-4);
221         for (double dt = 0; dt < 4000; dt += 60) {
222 
223             AbsoluteDate date = initialDate.shiftedBy(dt);
224             SpacecraftState state = propagator.propagate(date);
225 
226             // simple Earth fixed atmosphere
227             Vector3D p = state.getPosition();
228             Vector3D v = state.getVelocity();
229             Vector3D vAtm = Vector3D.crossProduct(earthRot, p);
230             Vector3D relativeVelocity = vAtm.subtract(v);
231 
232             Vector3D drag = s.dragAcceleration(state, 0.001, relativeVelocity, getDragParameters(s));
233             Assertions.assertEquals(0.0, Vector3D.angle(relativeVelocity, drag), 1.0e-15);
234 
235             Vector3D sunDirection = sun.getPosition(date, state.getFrame()).normalize();
236             Vector3D flux = new Vector3D(-4.56e-6, sunDirection);
237             Vector3D radiation = s.radiationPressureAcceleration(state, flux, getRadiationParameters(s));
238             Assertions.assertEquals(0.0, Vector3D.angle(flux, radiation), 1.0e-9);
239 
240         }
241 
242     }
243 
244     @Test
245     void testOnlyLiftWithoutReflection() {
246 
247         AbsoluteDate initialDate = propagator.getInitialState().getDate();
248         CelestialBody sun = CelestialBodyFactory.getSun();
249         BoxAndSolarArraySpacecraft s =
250             new BoxAndSolarArraySpacecraft(1.5, 3.5, 2.5, sun, 20.0, Vector3D.PLUS_J, 1.0, 1.0, 1.0, 0.0);
251 
252         Vector3D earthRot = new Vector3D(0.0, 0.0, 7.292115e-4);
253         for (double dt = 0; dt < 4000; dt += 60) {
254 
255             AbsoluteDate date = initialDate.shiftedBy(dt);
256             SpacecraftState state = propagator.propagate(date);
257 
258             // simple Earth fixed atmosphere
259             Vector3D p = state.getPosition();
260             Vector3D v = state.getVelocity();
261             Vector3D vAtm = Vector3D.crossProduct(earthRot, p);
262             Vector3D relativeVelocity = vAtm.subtract(v);
263 
264             Vector3D drag = s.dragAcceleration(state, 0.001, relativeVelocity, getDragParameters(s));
265             Assertions.assertTrue(Vector3D.angle(relativeVelocity, drag) > 0.167);
266             Assertions.assertTrue(Vector3D.angle(relativeVelocity, drag) < 0.736);
267 
268             Vector3D sunDirection = sun.getPosition(date, state.getFrame()).normalize();
269             Vector3D flux = new Vector3D(-4.56e-6, sunDirection);
270             Vector3D radiation = s.radiationPressureAcceleration(state, flux, getRadiationParameters(s));
271             Assertions.assertEquals(0.0, Vector3D.angle(flux, radiation), 1.0e-9);
272 
273         }
274 
275     }
276 
277     @Test
278     void testLiftVsNoLift()
279         throws NoSuchFieldException, SecurityException,
280                IllegalArgumentException, IllegalAccessException {
281 
282         // older implementation did not consider lift, so it really worked
283         // only for symmetrical shapes. For testing purposes, we will use a
284         // basic cubic shape without solar arrays and a relative atmosphere
285         // velocity either *exactly* facing a side or *exactly* along a main diagonal
286         List<Panel> facets = new ArrayList<>(7);
287         facets.add(new FixedPanel(Vector3D.MINUS_I, 3.0, false, 1.0, 1.0, 1.0, 0.0));
288         facets.add(new FixedPanel(Vector3D.PLUS_I,  3.0, false, 1.0, 1.0, 1.0, 0.0));
289         facets.add(new FixedPanel(Vector3D.MINUS_J, 3.0, false, 1.0, 1.0, 1.0, 0.0));
290         facets.add(new FixedPanel(Vector3D.PLUS_J,  3.0, false, 1.0, 1.0, 1.0, 0.0));
291         facets.add(new FixedPanel(Vector3D.MINUS_K, 3.0, false, 1.0, 1.0, 1.0, 0.0));
292         facets.add(new FixedPanel(Vector3D.PLUS_K,  3.0, false, 1.0, 1.0, 1.0, 0.0));
293         BoxAndSolarArraySpacecraft cube = new BoxAndSolarArraySpacecraft(facets);
294 
295         AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
296         Frame frame = FramesFactory.getEME2000();
297         Vector3D position = new Vector3D(1234567.8, 9876543.21, 121212.3434);
298         double density = 0.001;
299         SpacecraftState state = new SpacecraftState(new CartesianOrbit(new TimeStampedPVCoordinates(date, position, Vector3D.ZERO),
300                                                                        frame, Constants.EIGEN5C_EARTH_MU),
301                                                     new Attitude(frame,
302                                                                 new TimeStampedAngularCoordinates(date,
303                                                                                                   Rotation.IDENTITY,
304                                                                                                   Vector3D.ZERO,
305                                                                                                   Vector3D.ZERO))).withMass(1000.0);
306 
307         // head-on, there acceleration with lift should be twice acceleration without lift
308         Vector3D headOnVelocity = new Vector3D(2000, 0.0, 0.0);
309         Vector3D newHeadOnDrag  = cube.dragAcceleration(state, density, headOnVelocity, getDragParameters(cube));
310         Vector3D oldHeadOnDrag  = oldDragAcceleration(cube, state, density, headOnVelocity);
311         MatcherAssert.assertThat(newHeadOnDrag, OrekitMatchers.vectorCloseTo(oldHeadOnDrag.scalarMultiply(2), 1));
312 
313         // on an angle, the no lift implementation applies drag to the velocity direction
314         // instead of to the facet normal direction. In the symmetrical case, this implies
315         // it applied a single cos(θ) coefficient (projected surface reduction) instead
316         // of using cos²(θ) (projected surface reduction *and* normal component projection)
317         // and since molecule is reflected backward with the same velocity, this implies a
318         // factor 2 in linear momentum differences
319         Vector3D diagonalVelocity = new Vector3D(2000, 2000, 2000);
320         Vector3D newDiagDrag= cube.dragAcceleration(state, density, diagonalVelocity, getDragParameters(cube));
321         Vector3D oldDiagDrag = oldDragAcceleration(cube, state, density, diagonalVelocity);
322         double oldMissingCoeff = 2.0 / FastMath.sqrt(3.0);
323         Vector3D fixedOldDrag = new Vector3D(oldMissingCoeff, oldDiagDrag);
324         MatcherAssert.assertThat(newDiagDrag, OrekitMatchers.vectorCloseTo(fixedOldDrag, 1));
325 
326     }
327 
328     // this is a slightly adapted version of the pre-9.0 implementation
329     // Beware that this implementation is WRONG
330     private Vector3D oldDragAcceleration(final BoxAndSolarArraySpacecraft boxWithoutSolarArray,
331                                          final SpacecraftState state,
332                                          final double density, final Vector3D relativeVelocity)
333          throws IllegalArgumentException, IllegalAccessException,
334                 NoSuchFieldException, SecurityException {
335 
336         final double dragCoeff = boxWithoutSolarArray.getDragParametersDrivers().get(0).getValue();
337 
338         // relative velocity in spacecraft frame
339         final Vector3D v = state.getAttitude().getRotation().applyTo(relativeVelocity);
340 
341         // body facets contribution
342         double sv = 0;
343         for (final Panel panel : boxWithoutSolarArray.getPanels()) {
344             final double dot = Vector3D.dotProduct(panel.getNormal(state), v);
345             if (dot < 0) {
346                 // the facet intercepts the incoming flux
347                 sv -= panel.getArea() * dot;
348             }
349         }
350 
351         return new Vector3D(sv * density * dragCoeff / (2.0 * state.getMass()), relativeVelocity);
352 
353     }
354 
355     @Test
356     void testPlaneSpecularReflection() {
357 
358         AbsoluteDate initialDate = propagator.getInitialState().getDate();
359         CelestialBody sun = CelestialBodyFactory.getSun();
360         final Panel reflectingSolarArray = new PointingPanel(Vector3D.PLUS_J, sun, 20.0, 0.0, 0.0, 0.0, 1.0);
361         BoxAndSolarArraySpacecraft s = new BoxAndSolarArraySpacecraft(Collections.singletonList(reflectingSolarArray));
362 
363         for (double dt = 0; dt < 4000; dt += 60) {
364 
365             AbsoluteDate date = initialDate.shiftedBy(dt);
366             SpacecraftState state = propagator.propagate(date);
367 
368             Vector3D sunDirection = sun.getPosition(date, state.getFrame()).normalize();
369             Vector3D flux         = new Vector3D(-4.56e-6, sunDirection);
370             Vector3D acceleration = s.radiationPressureAcceleration(state, flux, getRadiationParameters(s));
371             Vector3D normal       = state.getAttitude().getRotation().applyInverseTo(reflectingSolarArray.getNormal(state));
372 
373             // solar array normal is slightly misaligned with Sun direction due to Sun being out of orbital plane
374             Assertions.assertEquals(15.1, FastMath.toDegrees(Vector3D.angle(sunDirection, normal)), 0.11);
375 
376             // radiation pressure is exactly opposed to solar array normal as there is only specular reflection
377             Assertions.assertEquals(180.0, FastMath.toDegrees(Vector3D.angle(acceleration, normal)), 1.0e-3);
378 
379         }
380 
381     }
382 
383     @Test
384     void testPlaneAbsorption() {
385 
386         AbsoluteDate initialDate = propagator.getInitialState().getDate();
387         CelestialBody sun = CelestialBodyFactory.getSun();
388         final Panel absorbingSolarArray = new PointingPanel(Vector3D.PLUS_J, sun, 20.0, 0.0, 0.0, 1.0, 0.0);
389         BoxAndSolarArraySpacecraft s = new BoxAndSolarArraySpacecraft(Collections.singletonList(absorbingSolarArray));
390 
391         for (double dt = 0; dt < 4000; dt += 60) {
392 
393             AbsoluteDate date = initialDate.shiftedBy(dt);
394             SpacecraftState state = propagator.propagate(date);
395 
396             Vector3D sunDirection = sun.getPosition(date, state.getFrame()).normalize();
397             Vector3D flux         = new Vector3D(-4.56e-6, sunDirection);
398             Vector3D acceleration = s.radiationPressureAcceleration(state, flux, getRadiationParameters(s));
399             Vector3D normal       = state.getAttitude().getRotation().applyInverseTo(absorbingSolarArray.getNormal(state));
400 
401             // solar array normal is slightly misaligned with Sun direction due to Sun being out of orbital plane
402             Assertions.assertEquals(15.1, FastMath.toDegrees(Vector3D.angle(sunDirection, normal)), 0.11);
403 
404             // radiation pressure is exactly opposed to Sun direction as there is only absorption
405             Assertions.assertEquals(180.0, FastMath.toDegrees(Vector3D.angle(acceleration, sunDirection)), 1.0e-3);
406 
407         }
408 
409     }
410 
411     /** Test solar array radiation acceleration with zero flux. */
412     @Test
413     void testNullIllumination() {
414         SpacecraftState state = propagator.getInitialState();
415         CelestialBody sun = CelestialBodyFactory.getSun();
416         final Panel absorbingSolarArray = new PointingPanel(Vector3D.PLUS_J, sun, 20.0, 0.0, 0.0, 1.0, 0.0);
417         BoxAndSolarArraySpacecraft s = new BoxAndSolarArraySpacecraft(Collections.singletonList(absorbingSolarArray));
418 
419         // "Field" the inputs using Binary64
420         Field<Binary64> field = Binary64Field.getInstance();
421         Binary64[] srpParam = getRadiationParameters(s, field);
422 
423         FieldSpacecraftState<Binary64> fState = new FieldSpacecraftState<>(field, state);
424         FieldVector3D<Binary64> flux = new FieldVector3D<Binary64>(field.getOne(),
425                         new Vector3D(Precision.SAFE_MIN / 2, Vector3D.PLUS_I));
426 
427 
428         FieldVector3D<Binary64> a = s.radiationPressureAcceleration(fState, flux, srpParam);
429         Assertions.assertEquals(0.0, a.getNorm().getReal(), Double.MIN_VALUE);
430     }
431 
432     /** Test forward/backward acceleration due to solar array radiation pressure. */
433     @Test
434     void testBackwardIllumination() {
435         SpacecraftState state = propagator.getInitialState();
436         CelestialBody sun = CelestialBodyFactory.getSun();
437         final Panel absorbingSolarArray = new PointingPanel(Vector3D.PLUS_J, sun, 20.0, 0.0, 0.0, 1.0, 0.0);
438         BoxAndSolarArraySpacecraft s = new BoxAndSolarArraySpacecraft(Collections.singletonList(absorbingSolarArray));
439 
440         // "Field" the inputs using Binary64
441         Field<Binary64> field = Binary64Field.getInstance();
442         Binary64[] srpParam = getRadiationParameters(s, field);
443 
444         FieldSpacecraftState<Binary64> fState = new FieldSpacecraftState<>(field, state);
445 
446         // Flux equal to SA normal
447         FieldVector3D<Binary64> flux = absorbingSolarArray.getNormal(fState);
448 
449         // Forward flux
450         FieldVector3D<Binary64> aPlus = s.radiationPressureAcceleration(fState, flux, srpParam);
451         // Backward flux
452         FieldVector3D<Binary64> aMinus = s.radiationPressureAcceleration(fState, flux.negate(), srpParam);
453 
454         Assertions.assertEquals(0.0, aPlus.add(aMinus).getNorm().getReal(), Double.MIN_VALUE);
455     }
456 
457     /** Test the functions computing drag and SRP acceleration and giving FieldVector3D outputs.
458      *  By comparing the "double" value with a "Binary64" implementation.
459      */
460     @Test
461     void testFieldAcceleration() {
462 
463         AbsoluteDate initialDate = propagator.getInitialState().getDate();
464         CelestialBody sun = CelestialBodyFactory.getSun();
465 
466         // Assuming simple Earth rotation model, constant density and flux
467         Vector3D earthRot = new Vector3D(0., 0., Constants.GRIM5C1_EARTH_ANGULAR_VELOCITY);
468         double density = 1.e-3;
469         double refFlux = 4.56e-6;
470 
471 
472         // Build a S/C box with non-nil coefficients so that the computation of the acceleration does not
473         // avoid any line of code
474         BoxAndSolarArraySpacecraft s =
475             new BoxAndSolarArraySpacecraft(1., 2., 3., sun, 20.0, Vector3D.PLUS_J, 2.0, 0.3, 0.5, 0.4);
476 
477         for (double dt = 0; dt < 4000; dt += 60) {
478             AbsoluteDate date = initialDate.shiftedBy(dt);
479             SpacecraftState state = propagator.propagate(date);
480 
481             // Data used in acceleration computation
482             Vector3D position = state.getPosition();
483             Vector3D velocity = state.getVelocity();
484             Vector3D vAtm = Vector3D.crossProduct(earthRot, position);
485             Vector3D relativeVelocity = vAtm.subtract(velocity);
486 
487             Frame frame = state.getFrame();
488             Vector3D flux = position.subtract(sun.getPosition(date, frame)).normalize().scalarMultiply(refFlux);
489 
490             // Acceleration in double
491             Vector3D aDrag = s.dragAcceleration(state, density, relativeVelocity, getDragParameters(s));
492             Vector3D aSrp = s.radiationPressureAcceleration(state, flux, getRadiationParameters(s));
493 
494             // "Field" the inputs using Binary64
495             Field<Binary64> field = Binary64Field.getInstance();
496             FieldSpacecraftState<Binary64> fState = new FieldSpacecraftState<>(field, state);
497 
498             FieldVector3D<Binary64> fluxF = new FieldVector3D<Binary64>(field.getOne(), flux);
499             Binary64 densityF = new Binary64(density);
500             FieldVector3D<Binary64> relativeVelocityF = new FieldVector3D<Binary64>(field.getOne(), relativeVelocity);
501 
502 
503             // Acceleration in Binary64
504             FieldVector3D<Binary64> aDragF = s.dragAcceleration(fState, densityF, relativeVelocityF, getDragParameters(s, field));
505             FieldVector3D<Binary64> aSrpF  = s.radiationPressureAcceleration(fState, fluxF, getRadiationParameters(s, field));
506             // Compare double and Binary64 accelerations
507             Assertions.assertEquals(0.0, Vector3D.distance(aDrag, aDragF.toVector3D()), Precision.EPSILON);
508             Assertions.assertEquals(0.0, Vector3D.distance(aSrp,  aSrpF.toVector3D()), Precision.EPSILON);
509         }
510     }
511 
512     /** Get drag parameters as double[]. */
513     private double[] getDragParameters(final BoxAndSolarArraySpacecraft basa) {
514         final List<ParameterDriver> drivers = basa.getDragParametersDrivers();
515         final double[] parameters = new double[drivers.size()];
516         for (int i = 0; i < drivers.size(); ++i) {
517             parameters[i] = drivers.get(i).getValue();
518         }
519         return parameters;
520     }
521 
522     /** Get radiation parameters as double[]. */
523     private double[] getRadiationParameters(final BoxAndSolarArraySpacecraft basa) {
524         final List<ParameterDriver> drivers = basa.getRadiationParametersDrivers();
525         final double[] parameters = new double[drivers.size()];
526         for (int i = 0; i < drivers.size(); ++i) {
527             parameters[i] = drivers.get(i).getValue();
528         }
529         return parameters;
530     }
531 
532     /** Get drag parameters as field[]. */
533     private <T extends CalculusFieldElement<T>> T[] getDragParameters(final BoxAndSolarArraySpacecraft basa,
534                                                                   final Field<T> field) {
535         final List<ParameterDriver> drivers = basa.getDragParametersDrivers();
536         final T[] parameters = MathArrays.buildArray(field, drivers.size());
537         for (int i = 0; i < drivers.size(); ++i) {
538             parameters[i] = field.getZero().add(drivers.get(i).getValue());
539         }
540         return parameters;
541     }
542 
543     /** Get radiation parameters as field[]. */
544     private <T extends CalculusFieldElement<T>> T[] getRadiationParameters(final BoxAndSolarArraySpacecraft basa,
545                                                                   final Field<T> field) {
546         final List<ParameterDriver> drivers = basa.getRadiationParametersDrivers();
547         final T[] parameters = MathArrays.buildArray(field, drivers.size());
548         for (int i = 0; i < drivers.size(); ++i) {
549             parameters[i] = field.getZero().add(drivers.get(i).getValue());
550         }
551         return parameters;
552     }
553 
554     @BeforeEach
555     public void setUp() {
556         try {
557         Utils.setDataRoot("regular-data");
558         mu  = 3.9860047e14;
559         double ae  = 6.378137e6;
560         double c20 = -1.08263e-3;
561         double c30 = 2.54e-6;
562         double c40 = 1.62e-6;
563         double c50 = 2.3e-7;
564         double c60 = -5.5e-7;
565 
566         AbsoluteDate date = new AbsoluteDate(new DateComponents(1970, 7, 1),
567                                              new TimeComponents(13, 59, 27.816),
568                                              TimeScalesFactory.getUTC());
569 
570         // Satellite position as circular parameters, raan chosen to have sun elevation with
571         // respect to orbit plane roughly evolving roughly from 15 to 15.2 degrees in the test range
572         Orbit circ =
573             new CircularOrbit(7178000.0, 0.5e-4, -0.5e-4, FastMath.toRadians(50.), FastMath.toRadians(280),
574                                    FastMath.toRadians(10.0), PositionAngleType.MEAN,
575                                    FramesFactory.getEME2000(), date, mu);
576         propagator =
577             new EcksteinHechlerPropagator(circ,
578                                           new LofOffset(circ.getFrame(), LOFType.LVLH_CCSDS),
579                                           ae, mu, c20, c30, c40, c50, c60);
580         } catch (OrekitException oe) {
581             Assertions.fail(oe.getLocalizedMessage());
582         }
583     }
584 
585     private double mu;
586     private Propagator propagator;
587 
588 }