1   /* Copyright 2022-2025 Romain Serra
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.drag;
18  
19  import org.hipparchus.CalculusFieldElement;
20  import org.hipparchus.Field;
21  import org.hipparchus.analysis.differentiation.DSFactory;
22  import org.hipparchus.analysis.differentiation.DerivativeStructure;
23  import org.hipparchus.analysis.differentiation.Gradient;
24  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
25  import org.hipparchus.geometry.euclidean.threed.Vector3D;
26  import org.hipparchus.util.FastMath;
27  import org.junit.jupiter.api.Assertions;
28  import org.junit.jupiter.api.Test;
29  import org.orekit.Utils;
30  import org.orekit.frames.Frame;
31  import org.orekit.frames.FramesFactory;
32  import org.orekit.models.earth.atmosphere.Atmosphere;
33  import org.orekit.propagation.FieldSpacecraftState;
34  import org.orekit.propagation.SpacecraftState;
35  import org.orekit.time.AbsoluteDate;
36  import org.orekit.time.FieldAbsoluteDate;
37  import org.orekit.utils.ParameterDriver;
38  import org.orekit.utils.TimeStampedFieldPVCoordinates;
39  
40  import java.util.Collections;
41  import java.util.List;
42  
43  import static org.junit.jupiter.api.Assertions.assertFalse;
44  import static org.junit.jupiter.api.Assertions.assertTrue;
45  import static org.mockito.Mockito.mock;
46  import static org.mockito.Mockito.when;
47  
48  class AbstractDragForceModelTest {
49  
50      @Test
51      void testGetGradientDensityWrtState() {
52          // GIVEN
53          final TestDragForceModel testDragForceModel = new TestDragForceModel(new TestAtmosphereModel());
54          final Frame frame = FramesFactory.getGCRF();
55          final Vector3D position = new Vector3D(100., 1, 0);
56          final int freeParameters = 3;
57          final FieldVector3D<Gradient> fieldPosition = new FieldVector3D<>(
58                  Gradient.variable(freeParameters, 0, position.getX()),
59                  Gradient.variable(freeParameters, 1, position.getY()),
60                  Gradient.variable(freeParameters, 2, position.getZ()));
61          final AbsoluteDate date = AbsoluteDate.ARBITRARY_EPOCH;
62          // WHEN
63          final Gradient actualGradient = testDragForceModel.getGradientDensityWrtState(date, frame, fieldPosition);
64          // THEN
65          final Gradient expectedGradient = testDragForceModel.getGradientDensityWrtStateUsingFiniteDifferences(date,
66                  frame, fieldPosition);
67          Assertions.assertEquals(expectedGradient.getValue(), actualGradient.getValue());
68          for (int i = 0; i < freeParameters; i++) {
69              Assertions.assertEquals(expectedGradient.getPartialDerivative(i), actualGradient.getPartialDerivative(i), 1e-10);
70          }
71      }
72  
73      @Test
74      void testGetDSDensityWrtState() {
75          // GIVEN
76          final TestDragForceModel testDragForceModel = new TestDragForceModel(new TestAtmosphereModel());
77          final Frame frame = FramesFactory.getGCRF();
78          final Vector3D position = new Vector3D(100., 1, 0);
79          final int freeParameters = 3;
80          final DSFactory factory = new DSFactory(freeParameters, 1);
81          final FieldVector3D<DerivativeStructure> fieldPosition = new FieldVector3D<>(
82                  factory.variable(0, position.getX()),
83                  factory.variable(1, position.getY()),
84                  factory.variable(2, position.getZ()));
85          final AbsoluteDate date = AbsoluteDate.ARBITRARY_EPOCH;
86          // WHEN
87          final DerivativeStructure actualDS = testDragForceModel.getDSDensityWrtState(date, frame, fieldPosition);
88          // THEN
89          final DerivativeStructure expectedDS = testDragForceModel.getDSDensityWrtStateUsingFiniteDifferences(date,
90                  frame, fieldPosition);
91          Assertions.assertEquals(expectedDS.getValue(), actualDS.getValue());
92          final double[] actualDerivatives = actualDS.getAllDerivatives();
93          final double[] expectedDerivatives = expectedDS.getAllDerivatives();
94          for (int i = 1; i < expectedDerivatives.length; i++) {
95              Assertions.assertEquals(expectedDerivatives[i], actualDerivatives[i], 1e-10);
96          }
97      }
98  
99      private static class TestAtmosphereModel implements Atmosphere {
100 
101         @Override
102         public Frame getFrame() {
103             return FramesFactory.getEME2000();
104         }
105 
106         @Override
107         public double getDensity(AbsoluteDate date, Vector3D position, Frame frame) {
108             return FastMath.exp(-position.getNormSq());
109         }
110 
111         @Override
112         public <T extends CalculusFieldElement<T>> T getDensity(FieldAbsoluteDate<T> date, FieldVector3D<T> position, Frame frame) {
113             return FastMath.exp(position.getNormSq().negate());
114         }
115     }
116 
117     private static class TestDragForceModel extends AbstractDragForceModel {
118 
119         protected TestDragForceModel(Atmosphere atmosphere) {
120             super(atmosphere);
121         }
122 
123         @Override
124         public Vector3D acceleration(SpacecraftState s, double[] parameters) {
125             return null;
126         }
127 
128         @Override
129         public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(FieldSpacecraftState<T> s, T[] parameters) {
130             return null;
131         }
132 
133         @Override
134         public List<ParameterDriver> getParametersDrivers() {
135             return Collections.emptyList();
136         }
137     }
138 
139     @Test
140     @SuppressWarnings("unchecked")
141     // This test was added to increase overall conditions coverage in the scope of issue 1453
142     void testIssue1453() {
143         // GIVEN
144         // Load Orekit data
145         Utils.setDataRoot("regular-data");
146 
147         // Create mock field
148         final Field<Gradient> fakeField = new Field<Gradient>() {
149             @Override
150             public Gradient getZero() {
151                 return new Gradient(0, 0);
152             }
153 
154             @Override
155             public Gradient getOne() {
156                 return new Gradient(1, 1);
157             }
158 
159             @Override
160             public Class<Gradient> getRuntimeClass() {
161                 return null;
162             }
163         };
164 
165         // Create mock date
166         final FieldAbsoluteDate<Gradient> fakeFieldAbsoluteDate = new FieldAbsoluteDate<>(fakeField);
167 
168         // Create fake PV
169         final FieldVector3D<Gradient> pos =
170                 new FieldVector3D<>(new FakeGradient(1, 1), new FakeGradient(2, 20), new FakeGradient(3, 30));
171 
172         final FieldVector3D<Gradient> vel =
173                 new FieldVector3D<>(new FakeGradient(4, 40), new FakeGradient(5, 50), new FakeGradient(6, 60));
174 
175         final TimeStampedFieldPVCoordinates<Gradient> fakePV =
176                 new TimeStampedFieldPVCoordinates<>(fakeFieldAbsoluteDate, pos, vel, FieldVector3D.getZero(fakeField));
177 
178         // Create mock mass
179         final Gradient mockMass = mock(Gradient.class);
180         when(mockMass.getFreeParameters()).thenReturn(6);
181 
182         // Create mock spacecraft state
183         final FieldSpacecraftState<Gradient> mockSpacecraftState = mock(FieldSpacecraftState.class);
184         when(mockSpacecraftState.getPVCoordinates()).thenReturn(fakePV);
185         when(mockSpacecraftState.getMass()).thenReturn(mockMass);
186 
187         // Create drag force
188         final DragForce dragForce = new DragForce(null,null);
189 
190         // WHEN & THEN
191         assertFalse(dragForce.isGradientStateDerivative(mockSpacecraftState));
192         ((FakeGradient) pos.getY()).setMutableGradient(new double[] { 0, 1 });
193         assertFalse(dragForce.isGradientStateDerivative(mockSpacecraftState));
194 
195         ((FakeGradient) pos.getZ()).setMutableGradient(new double[] { 0, 0, 1 });
196         assertFalse(dragForce.isGradientStateDerivative(mockSpacecraftState));
197 
198         ((FakeGradient) vel.getX()).setMutableGradient(new double[] { 0, 0, 0, 1 });
199         assertFalse(dragForce.isGradientStateDerivative(mockSpacecraftState));
200 
201         ((FakeGradient) vel.getY()).setMutableGradient(new double[] { 0, 0, 0, 0, 1 });
202         assertFalse(dragForce.isGradientStateDerivative(mockSpacecraftState));
203 
204         ((FakeGradient) vel.getZ()).setMutableGradient(new double[] { 0, 0, 0, 0, 0, 1 });
205         assertTrue(dragForce.isGradientStateDerivative(mockSpacecraftState));
206     }
207 
208     private static class FakeGradient extends Gradient {
209 
210         private final double   mutableValue;
211         private double[] mutableGradient;
212 
213         public FakeGradient(final double value, final double... gradient) {
214             super(value, gradient);
215             this.mutableValue    = value;
216             this.mutableGradient = gradient;
217         }
218 
219         @Override
220         public double getValue() {
221             return mutableValue;
222         }
223 
224         @Override
225         public double[] getGradient() {
226             return mutableGradient;
227         }
228 
229         public void setMutableGradient(final double[] mutableGradient) {
230             this.mutableGradient = mutableGradient;
231         }
232     }
233 
234 }