1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.control.indirect.adjoint;
18
19 import org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.analysis.differentiation.Gradient;
21 import org.hipparchus.analysis.differentiation.GradientField;
22 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
23 import org.hipparchus.geometry.euclidean.threed.Vector3D;
24 import org.hipparchus.util.Binary64;
25 import org.hipparchus.util.Binary64Field;
26 import org.hipparchus.util.MathArrays;
27 import org.junit.jupiter.api.Assertions;
28 import org.junit.jupiter.api.BeforeEach;
29 import org.junit.jupiter.api.Test;
30 import org.mockito.Mockito;
31 import org.orekit.Utils;
32 import org.orekit.bodies.CelestialBody;
33 import org.orekit.forces.gravity.SingleBodyAbsoluteAttraction;
34 import org.orekit.frames.Frame;
35 import org.orekit.frames.FramesFactory;
36 import org.orekit.propagation.FieldSpacecraftState;
37 import org.orekit.time.AbsoluteDate;
38 import org.orekit.time.DateComponents;
39 import org.orekit.time.FieldAbsoluteDate;
40 import org.orekit.time.TimeScalesFactory;
41 import org.orekit.utils.TimeStampedFieldPVCoordinates;
42
43 class CartesianAdjointSingleBodyTermTest {
44
45 @BeforeEach
46 public void setUp() {
47 Utils.setDataRoot("regular-data");
48 }
49
50 @Test
51 void testGetters() {
52
53 final CelestialBody body = getCelestialBody(Vector3D.MINUS_I);
54 final Frame expectedFrame = FramesFactory.getGCRF();
55 final CartesianAdjointSingleBodyTerm singleBodyTerm = new CartesianAdjointSingleBodyTerm(body.getGM(), body);
56
57 final double mu = singleBodyTerm.getMu();
58
59 Assertions.assertEquals(body.getGM(), mu);
60 }
61
62 @Test
63 @SuppressWarnings("unchecked")
64 void testGetPositionAdjointContribution() {
65
66 final CelestialBody body = getCelestialBody(Vector3D.MINUS_I);
67 final CartesianAdjointSingleBodyTerm singleBodyTerm = new CartesianAdjointSingleBodyTerm(body.getGM(), body);
68 final double[] adjoint = new double[6];
69 final double[] state = new double[6];
70 for (int i = 0; i < adjoint.length; i++) {
71 state[i] = -i+1;
72 adjoint[i] = i;
73 }
74 final AbsoluteDate date = new AbsoluteDate(new DateComponents(2015, 1, 1), TimeScalesFactory.getUTC());
75
76 final double[] contribution = singleBodyTerm.getPositionAdjointContribution(date, state, adjoint,
77 FramesFactory.getGCRF());
78
79 final SingleBodyAbsoluteAttraction singleBodyAbsoluteAttraction = new SingleBodyAbsoluteAttraction(body);
80 final int dimension = 3;
81 final GradientField field = GradientField.getField(dimension);
82 final FieldSpacecraftState<Gradient> mockedState = Mockito.mock(FieldSpacecraftState.class);
83 final FieldVector3D<Gradient> fieldPosition = new FieldVector3D<>(
84 Gradient.variable(dimension, 0, state[0]),
85 Gradient.variable(dimension, 1, state[1]),
86 Gradient.variable(dimension, 2, state[2]));
87 Mockito.when(mockedState.getDate()).thenReturn(new FieldAbsoluteDate<>(field, date));
88 Mockito.when(mockedState.getPosition()).thenReturn(fieldPosition);
89 final Gradient[] fieldMu = MathArrays.buildArray(field, 1);
90 fieldMu[0] = Gradient.constant(dimension, singleBodyTerm.getMu());
91 final FieldVector3D<Gradient> acceleration = singleBodyAbsoluteAttraction.acceleration(mockedState, fieldMu);
92 final double tolerance = 1e-12;
93 Assertions.assertEquals(-contribution[0], acceleration.getX().getGradient()[0] * adjoint[3]
94 + acceleration.getX().getGradient()[1] * adjoint[4] + acceleration.getX().getGradient()[2] * adjoint[5], tolerance);
95 Assertions.assertEquals(-contribution[1], acceleration.getY().getGradient()[0] * adjoint[3]
96 + acceleration.getY().getGradient()[1] * adjoint[4] + acceleration.getY().getGradient()[2] * adjoint[5], tolerance);
97 Assertions.assertEquals(-contribution[2], acceleration.getZ().getGradient()[0] * adjoint[3]
98 + acceleration.getZ().getGradient()[1] * adjoint[4] + acceleration.getZ().getGradient()[2] * adjoint[5], tolerance);
99 }
100
101 @Test
102 void testGetPositionAdjointFieldContribution() {
103
104 final CelestialBody celestialBody = getCelestialBody(Vector3D.ZERO);
105 final CartesianAdjointSingleBodyTerm cartesianAdjointSingleBodyTerm = new CartesianAdjointSingleBodyTerm(celestialBody.getGM(),
106 celestialBody);
107 final Binary64Field field = Binary64Field.getInstance();
108 final Binary64[] fieldAdjoint = MathArrays.buildArray(field, 6);
109 final Binary64[] fieldState = MathArrays.buildArray(field, 6);
110 for (int i = 0; i < fieldAdjoint.length; i++) {
111 fieldState[i] = field.getZero().newInstance(-i+1);
112 fieldAdjoint[i] = field.getZero().newInstance(i);
113 }
114 final Frame frame = celestialBody.getInertiallyOrientedFrame();
115 final FieldAbsoluteDate<Binary64> fieldDate = FieldAbsoluteDate.getArbitraryEpoch(field);
116
117 final Binary64[] fieldContribution = cartesianAdjointSingleBodyTerm.getPositionAdjointFieldContribution(fieldDate,
118 fieldState, fieldAdjoint, frame);
119
120 final CartesianAdjointKeplerianTerm keplerianTerm = new CartesianAdjointKeplerianTerm(celestialBody.getGM());
121 final Binary64[] contribution = keplerianTerm.getPositionAdjointFieldContribution(fieldDate, fieldState,
122 fieldAdjoint, frame);
123 for (int i = 0; i < contribution.length; i++) {
124 Assertions.assertEquals(fieldContribution[i], contribution[i]);
125 }
126 }
127
128 @Test
129 void testGetHamiltonianContribution() {
130
131 final CelestialBody celestialBody = getCelestialBody(Vector3D.ZERO);
132 final CartesianAdjointSingleBodyTerm cartesianAdjointSingleBodyTerm = new CartesianAdjointSingleBodyTerm(celestialBody.getGM(),
133 celestialBody);
134 final AbsoluteDate date = AbsoluteDate.ARBITRARY_EPOCH;
135 final Frame frame = FramesFactory.getGCRF();
136 final double[] adjoint = { -1, 2, -3, 4, -5, 6 };
137 final double[] state = {1, 3, 5, 7, 9, 11};
138
139 final double actualContribution = cartesianAdjointSingleBodyTerm.getHamiltonianContribution(date, state, adjoint, frame);
140
141 final double expectedContribution = new CartesianAdjointKeplerianTerm(celestialBody.getGM()).getHamiltonianContribution(date,
142 state, adjoint, frame);
143 Assertions.assertEquals(expectedContribution, actualContribution);
144 }
145
146 @Test
147 void testGetFieldHamiltonianContribution() {
148
149 final CelestialBody celestialBody = getCelestialBody(Vector3D.ZERO);
150 final CartesianAdjointSingleBodyTerm cartesianAdjointSingleBodyTerm = new CartesianAdjointSingleBodyTerm(celestialBody.getGM(),
151 celestialBody);
152 final Binary64Field field = Binary64Field.getInstance();
153 final Binary64[] fieldAdjoint = MathArrays.buildArray(field, 6);
154 final Binary64[] fieldState = MathArrays.buildArray(field, 6);
155 for (int i = 0; i < fieldAdjoint.length; i++) {
156 fieldState[i] = field.getZero().newInstance(-i + 1);
157 fieldAdjoint[i] = field.getZero().newInstance(i);
158 }
159 final Frame frame = celestialBody.getInertiallyOrientedFrame();
160 final FieldAbsoluteDate<Binary64> fieldDate = FieldAbsoluteDate.getArbitraryEpoch(field);
161
162 final Binary64 fieldContribution = cartesianAdjointSingleBodyTerm.getFieldHamiltonianContribution(fieldDate,
163 fieldState, fieldAdjoint, frame);
164
165 final double[] adjoint = new double[fieldAdjoint.length];
166 final double[] state = adjoint.clone();
167 for (int i = 0; i < fieldAdjoint.length; i++) {
168 state[i] = fieldState[i].getReal();
169 adjoint[i] = fieldAdjoint[i].getReal();
170 }
171 final double contribution = cartesianAdjointSingleBodyTerm.getHamiltonianContribution(fieldDate.toAbsoluteDate(),
172 state, adjoint, frame);
173 Assertions.assertEquals(contribution, fieldContribution.getReal());
174 }
175
176 private static CelestialBody getCelestialBody(final Vector3D position) {
177 return new CelestialBody() {
178 @Override
179 public Frame getInertiallyOrientedFrame() {
180 return FramesFactory.getGCRF();
181 }
182
183 @Override
184 public Frame getBodyOrientedFrame() {
185 return null;
186 }
187
188 @Override
189 public String getName() {
190 return "";
191 }
192
193 @Override
194 public double getGM() {
195 return 1.;
196 }
197
198 @Override
199 public Vector3D getPosition(AbsoluteDate date, Frame frame) {
200 return position;
201 }
202
203 @Override
204 public <T extends CalculusFieldElement<T>> FieldVector3D<T> getPosition(FieldAbsoluteDate<T> date, Frame frame) {
205 return new FieldVector3D<>(date.getField(), getPosition(date.toAbsoluteDate(), frame));
206 }
207
208 @Override
209 public <T extends CalculusFieldElement<T>> TimeStampedFieldPVCoordinates<T> getPVCoordinates(FieldAbsoluteDate<T> date, Frame frame) {
210 return null;
211 }
212 };
213 }
214
215 }