1   /* Copyright 2002-2023 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      public 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     public 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     public 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.getPVCoordinates().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     public 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.getPVCoordinates().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     public 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)),
306                                                     1000.0);
307 
308         // head-on, there acceleration with lift should be twice acceleration without lift
309         Vector3D headOnVelocity = new Vector3D(2000, 0.0, 0.0);
310         Vector3D newHeadOnDrag  = cube.dragAcceleration(state, density, headOnVelocity, getDragParameters(cube));
311         Vector3D oldHeadOnDrag  = oldDragAcceleration(cube, state, density, headOnVelocity);
312         MatcherAssert.assertThat(newHeadOnDrag, OrekitMatchers.vectorCloseTo(oldHeadOnDrag.scalarMultiply(2), 1));
313 
314         // on an angle, the no lift implementation applies drag to the velocity direction
315         // instead of to the facet normal direction. In the symmetrical case, this implies
316         // it applied a single cos(θ) coefficient (projected surface reduction) instead
317         // of using cos²(θ) (projected surface reduction *and* normal component projection)
318         // and since molecule is reflected backward with the same velocity, this implies a
319         // factor 2 in linear momentum differences
320         Vector3D diagonalVelocity = new Vector3D(2000, 2000, 2000);
321         Vector3D newDiagDrag= cube.dragAcceleration(state, density, diagonalVelocity, getDragParameters(cube));
322         Vector3D oldDiagDrag = oldDragAcceleration(cube, state, density, diagonalVelocity);
323         double oldMissingCoeff = 2.0 / FastMath.sqrt(3.0);
324         Vector3D fixedOldDrag = new Vector3D(oldMissingCoeff, oldDiagDrag);
325         MatcherAssert.assertThat(newDiagDrag, OrekitMatchers.vectorCloseTo(fixedOldDrag, 1));
326 
327     }
328 
329     // this is a slightly adapted version of the pre-9.0 implementation
330     // Beware that this implementation is WRONG
331     private Vector3D oldDragAcceleration(final BoxAndSolarArraySpacecraft boxWithoutSolarArray,
332                                          final SpacecraftState state,
333                                          final double density, final Vector3D relativeVelocity)
334          throws IllegalArgumentException, IllegalAccessException,
335                 NoSuchFieldException, SecurityException {
336 
337         final double dragCoeff = boxWithoutSolarArray.getDragParametersDrivers().get(0).getValue();
338 
339         // relative velocity in spacecraft frame
340         final Vector3D v = state.getAttitude().getRotation().applyTo(relativeVelocity);
341 
342         // body facets contribution
343         double sv = 0;
344         for (final Panel panel : boxWithoutSolarArray.getPanels()) {
345             final double dot = Vector3D.dotProduct(panel.getNormal(state), v);
346             if (dot < 0) {
347                 // the facet intercepts the incoming flux
348                 sv -= panel.getArea() * dot;
349             }
350         }
351 
352         return new Vector3D(sv * density * dragCoeff / (2.0 * state.getMass()), relativeVelocity);
353 
354     }
355 
356     @Test
357     public void testPlaneSpecularReflection() {
358 
359         AbsoluteDate initialDate = propagator.getInitialState().getDate();
360         CelestialBody sun = CelestialBodyFactory.getSun();
361         final Panel reflectingSolarArray = new PointingPanel(Vector3D.PLUS_J, sun, 20.0, 0.0, 0.0, 0.0, 1.0);
362         BoxAndSolarArraySpacecraft s = new BoxAndSolarArraySpacecraft(Collections.singletonList(reflectingSolarArray));
363 
364         for (double dt = 0; dt < 4000; dt += 60) {
365 
366             AbsoluteDate date = initialDate.shiftedBy(dt);
367             SpacecraftState state = propagator.propagate(date);
368 
369             Vector3D sunDirection = sun.getPosition(date, state.getFrame()).normalize();
370             Vector3D flux         = new Vector3D(-4.56e-6, sunDirection);
371             Vector3D acceleration = s.radiationPressureAcceleration(state, flux, getRadiationParameters(s));
372             Vector3D normal       = state.getAttitude().getRotation().applyInverseTo(reflectingSolarArray.getNormal(state));
373 
374             // solar array normal is slightly misaligned with Sun direction due to Sun being out of orbital plane
375             Assertions.assertEquals(15.1, FastMath.toDegrees(Vector3D.angle(sunDirection, normal)), 0.11);
376 
377             // radiation pressure is exactly opposed to solar array normal as there is only specular reflection
378             Assertions.assertEquals(180.0, FastMath.toDegrees(Vector3D.angle(acceleration, normal)), 1.0e-3);
379 
380         }
381 
382     }
383 
384     @Test
385     public void testPlaneAbsorption() {
386 
387         AbsoluteDate initialDate = propagator.getInitialState().getDate();
388         CelestialBody sun = CelestialBodyFactory.getSun();
389         final Panel absorbingSolarArray = new PointingPanel(Vector3D.PLUS_J, sun, 20.0, 0.0, 0.0, 1.0, 0.0);
390         BoxAndSolarArraySpacecraft s = new BoxAndSolarArraySpacecraft(Collections.singletonList(absorbingSolarArray));
391 
392         for (double dt = 0; dt < 4000; dt += 60) {
393 
394             AbsoluteDate date = initialDate.shiftedBy(dt);
395             SpacecraftState state = propagator.propagate(date);
396 
397             Vector3D sunDirection = sun.getPosition(date, state.getFrame()).normalize();
398             Vector3D flux         = new Vector3D(-4.56e-6, sunDirection);
399             Vector3D acceleration = s.radiationPressureAcceleration(state, flux, getRadiationParameters(s));
400             Vector3D normal       = state.getAttitude().getRotation().applyInverseTo(absorbingSolarArray.getNormal(state));
401 
402             // solar array normal is slightly misaligned with Sun direction due to Sun being out of orbital plane
403             Assertions.assertEquals(15.1, FastMath.toDegrees(Vector3D.angle(sunDirection, normal)), 0.11);
404 
405             // radiation pressure is exactly opposed to Sun direction as there is only absorption
406             Assertions.assertEquals(180.0, FastMath.toDegrees(Vector3D.angle(acceleration, sunDirection)), 1.0e-3);
407 
408         }
409 
410     }
411 
412     /** Test solar array radiation acceleration with zero flux. */
413     @Test
414     public void testNullIllumination() {
415         SpacecraftState state = propagator.getInitialState();
416         CelestialBody sun = CelestialBodyFactory.getSun();
417         final Panel absorbingSolarArray = new PointingPanel(Vector3D.PLUS_J, sun, 20.0, 0.0, 0.0, 1.0, 0.0);
418         BoxAndSolarArraySpacecraft s = new BoxAndSolarArraySpacecraft(Collections.singletonList(absorbingSolarArray));
419 
420         // "Field" the inputs using Binary64
421         Field<Binary64> field = Binary64Field.getInstance();
422         Binary64[] srpParam = getRadiationParameters(s, field);
423 
424         FieldSpacecraftState<Binary64> fState = new FieldSpacecraftState<>(field, state);
425         FieldVector3D<Binary64> flux = new FieldVector3D<Binary64>(field.getOne(),
426                         new Vector3D(Precision.SAFE_MIN / 2, Vector3D.PLUS_I));
427 
428 
429         FieldVector3D<Binary64> a = s.radiationPressureAcceleration(fState, flux, srpParam);
430         Assertions.assertEquals(0.0, a.getNorm().getReal(), Double.MIN_VALUE);
431     }
432 
433     /** Test forward/backward acceleration due to solar array radiation pressure. */
434     @Test
435     public void testBackwardIllumination() {
436         SpacecraftState state = propagator.getInitialState();
437         CelestialBody sun = CelestialBodyFactory.getSun();
438         final Panel absorbingSolarArray = new PointingPanel(Vector3D.PLUS_J, sun, 20.0, 0.0, 0.0, 1.0, 0.0);
439         BoxAndSolarArraySpacecraft s = new BoxAndSolarArraySpacecraft(Collections.singletonList(absorbingSolarArray));
440 
441         // "Field" the inputs using Binary64
442         Field<Binary64> field = Binary64Field.getInstance();
443         Binary64[] srpParam = getRadiationParameters(s, field);
444 
445         FieldSpacecraftState<Binary64> fState = new FieldSpacecraftState<>(field, state);
446 
447         // Flux equal to SA normal
448         FieldVector3D<Binary64> flux = absorbingSolarArray.getNormal(fState);
449 
450         // Forward flux
451         FieldVector3D<Binary64> aPlus = s.radiationPressureAcceleration(fState, flux, srpParam);
452         // Backward flux
453         FieldVector3D<Binary64> aMinus = s.radiationPressureAcceleration(fState, flux.negate(), srpParam);
454 
455         Assertions.assertEquals(0.0, aPlus.add(aMinus).getNorm().getReal(), Double.MIN_VALUE);
456     }
457 
458     /** Test the functions computing drag and SRP acceleration and giving FieldVector3D outputs.
459      *  By comparing the "double" value with a "Binary64" implementation.
460      */
461     @Test
462     public void testFieldAcceleration() {
463 
464         AbsoluteDate initialDate = propagator.getInitialState().getDate();
465         CelestialBody sun = CelestialBodyFactory.getSun();
466 
467         // Assuming simple Earth rotation model, constant density and flux
468         Vector3D earthRot = new Vector3D(0., 0., Constants.GRIM5C1_EARTH_ANGULAR_VELOCITY);
469         double density = 1.e-3;
470         double refFlux = 4.56e-6;
471 
472 
473         // Build a S/C box with non-nil coefficients so that the computation of the acceleration does not
474         // avoid any line of code
475         BoxAndSolarArraySpacecraft s =
476             new BoxAndSolarArraySpacecraft(1., 2., 3., sun, 20.0, Vector3D.PLUS_J, 2.0, 0.3, 0.5, 0.4);
477 
478         for (double dt = 0; dt < 4000; dt += 60) {
479             AbsoluteDate date = initialDate.shiftedBy(dt);
480             SpacecraftState state = propagator.propagate(date);
481 
482             // Data used in acceleration computation
483             Vector3D position = state.getPosition();
484             Vector3D velocity = state.getPVCoordinates().getVelocity();
485             Vector3D vAtm = Vector3D.crossProduct(earthRot, position);
486             Vector3D relativeVelocity = vAtm.subtract(velocity);
487 
488             Frame frame = state.getFrame();
489             Vector3D flux = position.subtract(sun.getPosition(date, frame)).normalize().scalarMultiply(refFlux);
490 
491             // Acceleration in double
492             Vector3D aDrag = s.dragAcceleration(state, density, relativeVelocity, getDragParameters(s));
493             Vector3D aSrp = s.radiationPressureAcceleration(state, flux, getRadiationParameters(s));
494 
495             // "Field" the inputs using Binary64
496             Field<Binary64> field = Binary64Field.getInstance();
497             FieldSpacecraftState<Binary64> fState = new FieldSpacecraftState<>(field, state);
498 
499             FieldVector3D<Binary64> fluxF = new FieldVector3D<Binary64>(field.getOne(), flux);
500             Binary64 densityF = new Binary64(density);
501             FieldVector3D<Binary64> relativeVelocityF = new FieldVector3D<Binary64>(field.getOne(), relativeVelocity);
502 
503 
504             // Acceleration in Binary64
505             FieldVector3D<Binary64> aDragF = s.dragAcceleration(fState, densityF, relativeVelocityF, getDragParameters(s, field));
506             FieldVector3D<Binary64> aSrpF  = s.radiationPressureAcceleration(fState, fluxF, getRadiationParameters(s, field));
507             // Compare double and Binary64 accelerations
508             Assertions.assertEquals(0.0, Vector3D.distance(aDrag, aDragF.toVector3D()), Precision.EPSILON);
509             Assertions.assertEquals(0.0, Vector3D.distance(aSrp,  aSrpF.toVector3D()), Precision.EPSILON);
510         }
511     }
512 
513     /** Get drag parameters as double[]. */
514     private double[] getDragParameters(final BoxAndSolarArraySpacecraft basa) {
515         final List<ParameterDriver> drivers = basa.getDragParametersDrivers();
516         final double[] parameters = new double[drivers.size()];
517         for (int i = 0; i < drivers.size(); ++i) {
518             parameters[i] = drivers.get(i).getValue();
519         }
520         return parameters;
521     }
522 
523     /** Get radiation parameters as double[]. */
524     private double[] getRadiationParameters(final BoxAndSolarArraySpacecraft basa) {
525         final List<ParameterDriver> drivers = basa.getRadiationParametersDrivers();
526         final double[] parameters = new double[drivers.size()];
527         for (int i = 0; i < drivers.size(); ++i) {
528             parameters[i] = drivers.get(i).getValue();
529         }
530         return parameters;
531     }
532 
533     /** Get drag parameters as field[]. */
534     private <T extends CalculusFieldElement<T>> T[] getDragParameters(final BoxAndSolarArraySpacecraft basa,
535                                                                   final Field<T> field) {
536         final List<ParameterDriver> drivers = basa.getDragParametersDrivers();
537         final T[] parameters = MathArrays.buildArray(field, drivers.size());
538         for (int i = 0; i < drivers.size(); ++i) {
539             parameters[i] = field.getZero().add(drivers.get(i).getValue());
540         }
541         return parameters;
542     }
543 
544     /** Get radiation parameters as field[]. */
545     private <T extends CalculusFieldElement<T>> T[] getRadiationParameters(final BoxAndSolarArraySpacecraft basa,
546                                                                   final Field<T> field) {
547         final List<ParameterDriver> drivers = basa.getRadiationParametersDrivers();
548         final T[] parameters = MathArrays.buildArray(field, drivers.size());
549         for (int i = 0; i < drivers.size(); ++i) {
550             parameters[i] = field.getZero().add(drivers.get(i).getValue());
551         }
552         return parameters;
553     }
554 
555     @BeforeEach
556     public void setUp() {
557         try {
558         Utils.setDataRoot("regular-data");
559         mu  = 3.9860047e14;
560         double ae  = 6.378137e6;
561         double c20 = -1.08263e-3;
562         double c30 = 2.54e-6;
563         double c40 = 1.62e-6;
564         double c50 = 2.3e-7;
565         double c60 = -5.5e-7;
566 
567         AbsoluteDate date = new AbsoluteDate(new DateComponents(1970, 7, 1),
568                                              new TimeComponents(13, 59, 27.816),
569                                              TimeScalesFactory.getUTC());
570 
571         // Satellite position as circular parameters, raan chosen to have sun elevation with
572         // respect to orbit plane roughly evolving roughly from 15 to 15.2 degrees in the test range
573         Orbit circ =
574             new CircularOrbit(7178000.0, 0.5e-4, -0.5e-4, FastMath.toRadians(50.), FastMath.toRadians(280),
575                                    FastMath.toRadians(10.0), PositionAngleType.MEAN,
576                                    FramesFactory.getEME2000(), date, mu);
577         propagator =
578             new EcksteinHechlerPropagator(circ,
579                                           new LofOffset(circ.getFrame(), LOFType.LVLH_CCSDS),
580                                           ae, mu, c20, c30, c40, c50, c60);
581         } catch (OrekitException oe) {
582             Assertions.fail(oe.getLocalizedMessage());
583         }
584     }
585 
586     private double mu;
587     private Propagator propagator;
588 
589 }