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