1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.forces;
18
19 import org.hamcrest.MatcherAssert;
20 import org.hipparchus.analysis.differentiation.DerivativeStructure;
21 import org.hipparchus.analysis.differentiation.Gradient;
22 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
23 import org.hipparchus.geometry.euclidean.threed.Vector3D;
24 import org.hipparchus.util.Precision;
25 import org.junit.jupiter.api.Assertions;
26 import org.orekit.OrekitMatchers;
27 import org.orekit.attitudes.AttitudeProvider;
28 import org.orekit.propagation.FieldSpacecraftState;
29 import org.orekit.propagation.SpacecraftState;
30
31
32 public abstract class AbstractLegacyForceModelTest extends AbstractForceModelTest {
33
34 protected abstract FieldVector3D<DerivativeStructure> accelerationDerivatives(ForceModel forceModel,
35 FieldSpacecraftState<DerivativeStructure> state);
36
37 protected abstract FieldVector3D<Gradient> accelerationDerivativesGradient(ForceModel forceModel,
38 FieldSpacecraftState<Gradient> state);
39
40 protected void checkStateJacobianVs80Implementation(final SpacecraftState state, final ForceModel forceModel,
41 final AttitudeProvider attitudeProvider,
42 final double checkTolerance, final boolean print) {
43 FieldSpacecraftState<DerivativeStructure> fState = toDS(state, attitudeProvider);
44 FieldVector3D<DerivativeStructure> dsNew = forceModel.acceleration(fState,
45 forceModel.getParameters(fState.getDate().getField(), fState.getDate()));
46 FieldVector3D<DerivativeStructure> dsOld = accelerationDerivatives(forceModel, fState);
47 Vector3D dFdPXRef = new Vector3D(dsOld.getX().getPartialDerivative(1, 0, 0, 0, 0, 0),
48 dsOld.getY().getPartialDerivative(1, 0, 0, 0, 0, 0),
49 dsOld.getZ().getPartialDerivative(1, 0, 0, 0, 0, 0));
50 Vector3D dFdPXRes = new Vector3D(dsNew.getX().getPartialDerivative(1, 0, 0, 0, 0, 0),
51 dsNew.getY().getPartialDerivative(1, 0, 0, 0, 0, 0),
52 dsNew.getZ().getPartialDerivative(1, 0, 0, 0, 0, 0));
53 Vector3D dFdPYRef = new Vector3D(dsOld.getX().getPartialDerivative(0, 1, 0, 0, 0, 0),
54 dsOld.getY().getPartialDerivative(0, 1, 0, 0, 0, 0),
55 dsOld.getZ().getPartialDerivative(0, 1, 0, 0, 0, 0));
56 Vector3D dFdPYRes = new Vector3D(dsNew.getX().getPartialDerivative(0, 1, 0, 0, 0, 0),
57 dsNew.getY().getPartialDerivative(0, 1, 0, 0, 0, 0),
58 dsNew.getZ().getPartialDerivative(0, 1, 0, 0, 0, 0));
59 Vector3D dFdPZRef = new Vector3D(dsOld.getX().getPartialDerivative(0, 0, 1, 0, 0, 0),
60 dsOld.getY().getPartialDerivative(0, 0, 1, 0, 0, 0),
61 dsOld.getZ().getPartialDerivative(0, 0, 1, 0, 0, 0));
62 Vector3D dFdPZRes = new Vector3D(dsNew.getX().getPartialDerivative(0, 0, 1, 0, 0, 0),
63 dsNew.getY().getPartialDerivative(0, 0, 1, 0, 0, 0),
64 dsNew.getZ().getPartialDerivative(0, 0, 1, 0, 0, 0));
65 Vector3D dFdVXRef = new Vector3D(dsOld.getX().getPartialDerivative(0, 0, 0, 1, 0, 0),
66 dsOld.getY().getPartialDerivative(0, 0, 0, 1, 0, 0),
67 dsOld.getZ().getPartialDerivative(0, 0, 0, 1, 0, 0));
68 Vector3D dFdVXRes = new Vector3D(dsNew.getX().getPartialDerivative(0, 0, 0, 1, 0, 0),
69 dsNew.getY().getPartialDerivative(0, 0, 0, 1, 0, 0),
70 dsNew.getZ().getPartialDerivative(0, 0, 0, 1, 0, 0));
71 Vector3D dFdVYRef = new Vector3D(dsOld.getX().getPartialDerivative(0, 0, 0, 0, 1, 0),
72 dsOld.getY().getPartialDerivative(0, 0, 0, 0, 1, 0),
73 dsOld.getZ().getPartialDerivative(0, 0, 0, 0, 1, 0));
74 Vector3D dFdVYRes = new Vector3D(dsNew.getX().getPartialDerivative(0, 0, 0, 0, 1, 0),
75 dsNew.getY().getPartialDerivative(0, 0, 0, 0, 1, 0),
76 dsNew.getZ().getPartialDerivative(0, 0, 0, 0, 1, 0));
77 Vector3D dFdVZRef = new Vector3D(dsOld.getX().getPartialDerivative(0, 0, 0, 0, 0, 1),
78 dsOld.getY().getPartialDerivative(0, 0, 0, 0, 0, 1),
79 dsOld.getZ().getPartialDerivative(0, 0, 0, 0, 0, 1));
80 Vector3D dFdVZRes = new Vector3D(dsNew.getX().getPartialDerivative(0, 0, 0, 0, 0, 1),
81 dsNew.getY().getPartialDerivative(0, 0, 0, 0, 0, 1),
82 dsNew.getZ().getPartialDerivative(0, 0, 0, 0, 0, 1));
83 if (print) {
84 System.out.println("dF/dPX ref norm: " + dFdPXRef.getNorm() + ", abs error: " + Vector3D.distance(dFdPXRef, dFdPXRes) + ", rel error: " + (Vector3D.distance(dFdPXRef, dFdPXRes) / dFdPXRef.getNorm()));
85 System.out.println("dF/dPY ref norm: " + dFdPYRef.getNorm() + ", abs error: " + Vector3D.distance(dFdPYRef, dFdPYRes) + ", rel error: " + (Vector3D.distance(dFdPYRef, dFdPYRes) / dFdPYRef.getNorm()));
86 System.out.println("dF/dPZ ref norm: " + dFdPZRef.getNorm() + ", abs error: " + Vector3D.distance(dFdPZRef, dFdPZRes) + ", rel error: " + (Vector3D.distance(dFdPZRef, dFdPZRes) / dFdPZRef.getNorm()));
87 System.out.println("dF/dVX ref norm: " + dFdVXRef.getNorm() + ", abs error: " + Vector3D.distance(dFdVXRef, dFdVXRes) + ", rel error: " + (Vector3D.distance(dFdVXRef, dFdVXRes) / dFdVXRef.getNorm()));
88 System.out.println("dF/dVY ref norm: " + dFdVYRef.getNorm() + ", abs error: " + Vector3D.distance(dFdVYRef, dFdVYRes) + ", rel error: " + (Vector3D.distance(dFdVYRef, dFdVYRes) / dFdVYRef.getNorm()));
89 System.out.println("dF/dVZ ref norm: " + dFdVZRef.getNorm() + ", abs error: " + Vector3D.distance(dFdVZRef, dFdVZRes) + ", rel error: " + (Vector3D.distance(dFdVZRef, dFdVZRes) / dFdVZRef.getNorm()));
90 }
91 checkdFdP(dFdPXRef, dFdPXRes, checkTolerance);
92 checkdFdP(dFdPYRef, dFdPYRes, checkTolerance);
93 checkdFdP(dFdPZRef, dFdPZRes, checkTolerance);
94 checkdFdP(dFdVXRef, dFdVXRes, checkTolerance);
95 checkdFdP(dFdVYRef, dFdVYRes, checkTolerance);
96 checkdFdP(dFdVZRef, dFdVZRes, checkTolerance);
97
98 }
99
100 protected void checkStateJacobianVs80ImplementationGradient(final SpacecraftState state, final ForceModel forceModel,
101 final AttitudeProvider attitudeProvider,
102 final double checkTolerance, final boolean print)
103 {
104 FieldSpacecraftState<Gradient> fState = toGradient(state, attitudeProvider);
105 FieldVector3D<Gradient> gNew = forceModel.acceleration(fState,
106 forceModel.getParameters(fState.getDate().getField(), fState.getDate()));
107 FieldVector3D<Gradient> gOld = accelerationDerivativesGradient(forceModel, fState);
108 Vector3D dFdPXRef = new Vector3D(gOld.getX().getPartialDerivative(0),
109 gOld.getY().getPartialDerivative(0),
110 gOld.getZ().getPartialDerivative(0));
111 Vector3D dFdPXRes = new Vector3D(gNew.getX().getPartialDerivative(0),
112 gNew.getY().getPartialDerivative(0),
113 gNew.getZ().getPartialDerivative(0));
114 Vector3D dFdPYRef = new Vector3D(gOld.getX().getPartialDerivative(1),
115 gOld.getY().getPartialDerivative(1),
116 gOld.getZ().getPartialDerivative(1));
117 Vector3D dFdPYRes = new Vector3D(gNew.getX().getPartialDerivative(1),
118 gNew.getY().getPartialDerivative(1),
119 gNew.getZ().getPartialDerivative(1));
120 Vector3D dFdPZRef = new Vector3D(gOld.getX().getPartialDerivative(2),
121 gOld.getY().getPartialDerivative(2),
122 gOld.getZ().getPartialDerivative(2));
123 Vector3D dFdPZRes = new Vector3D(gNew.getX().getPartialDerivative(2),
124 gNew.getY().getPartialDerivative(2),
125 gNew.getZ().getPartialDerivative(2));
126 Vector3D dFdVXRef = new Vector3D(gOld.getX().getPartialDerivative(3),
127 gOld.getY().getPartialDerivative(3),
128 gOld.getZ().getPartialDerivative(3));
129 Vector3D dFdVXRes = new Vector3D(gNew.getX().getPartialDerivative(3),
130 gNew.getY().getPartialDerivative(3),
131 gNew.getZ().getPartialDerivative(3));
132 Vector3D dFdVYRef = new Vector3D(gOld.getX().getPartialDerivative(4),
133 gOld.getY().getPartialDerivative(4),
134 gOld.getZ().getPartialDerivative(4));
135 Vector3D dFdVYRes = new Vector3D(gNew.getX().getPartialDerivative(4),
136 gNew.getY().getPartialDerivative(4),
137 gNew.getZ().getPartialDerivative(4));
138 Vector3D dFdVZRef = new Vector3D(gOld.getX().getPartialDerivative(5),
139 gOld.getY().getPartialDerivative(5),
140 gOld.getZ().getPartialDerivative(5));
141 Vector3D dFdVZRes = new Vector3D(gNew.getX().getPartialDerivative(5),
142 gNew.getY().getPartialDerivative(5),
143 gNew.getZ().getPartialDerivative(5));
144 if (print) {
145 System.out.println("dF/dPX ref norm: " + dFdPXRef.getNorm() + ", abs error: " + Vector3D.distance(dFdPXRef, dFdPXRes) + ", rel error: " + (Vector3D.distance(dFdPXRef, dFdPXRes) / dFdPXRef.getNorm()));
146 System.out.println("dF/dPY ref norm: " + dFdPYRef.getNorm() + ", abs error: " + Vector3D.distance(dFdPYRef, dFdPYRes) + ", rel error: " + (Vector3D.distance(dFdPYRef, dFdPYRes) / dFdPYRef.getNorm()));
147 System.out.println("dF/dPZ ref norm: " + dFdPZRef.getNorm() + ", abs error: " + Vector3D.distance(dFdPZRef, dFdPZRes) + ", rel error: " + (Vector3D.distance(dFdPZRef, dFdPZRes) / dFdPZRef.getNorm()));
148 System.out.println("dF/dVX ref norm: " + dFdVXRef.getNorm() + ", abs error: " + Vector3D.distance(dFdVXRef, dFdVXRes) + ", rel error: " + (Vector3D.distance(dFdVXRef, dFdVXRes) / dFdVXRef.getNorm()));
149 System.out.println("dF/dVY ref norm: " + dFdVYRef.getNorm() + ", abs error: " + Vector3D.distance(dFdVYRef, dFdVYRes) + ", rel error: " + (Vector3D.distance(dFdVYRef, dFdVYRes) / dFdVYRef.getNorm()));
150 System.out.println("dF/dVZ ref norm: " + dFdVZRef.getNorm() + ", abs error: " + Vector3D.distance(dFdVZRef, dFdVZRes) + ", rel error: " + (Vector3D.distance(dFdVZRef, dFdVZRes) / dFdVZRef.getNorm()));
151 }
152 checkdFdP(dFdPXRef, dFdPXRes, checkTolerance);
153 checkdFdP(dFdPYRef, dFdPYRes, checkTolerance);
154 checkdFdP(dFdPZRef, dFdPZRes, checkTolerance);
155 checkdFdP(dFdVXRef, dFdVXRes, checkTolerance);
156 checkdFdP(dFdVYRef, dFdVYRes, checkTolerance);
157 checkdFdP(dFdVZRef, dFdVZRes, checkTolerance);
158
159 }
160
161 private void checkdFdP(final Vector3D reference, final Vector3D result, final double checkTolerance) {
162 MatcherAssert.assertThat(result,
163 OrekitMatchers.vectorCloseTo(reference, Precision.SAFE_MIN, checkTolerance));
164 if (reference.getNorm() == 0) {
165
166
167 Assertions.assertEquals(0, result.getNorm(), Precision.SAFE_MIN);
168 } else {
169 Assertions.assertEquals(0, Vector3D.distance(reference, result), checkTolerance * reference.getNorm());
170 }
171 }
172
173 }
174
175