1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.forces.gravity;
18
19 import java.io.Serializable;
20 import java.lang.reflect.InvocationTargetException;
21
22 import org.hipparchus.Field;
23 import org.hipparchus.analysis.differentiation.DSFactory;
24 import org.hipparchus.analysis.differentiation.DerivativeStructure;
25 import org.hipparchus.analysis.differentiation.Gradient;
26 import org.hipparchus.analysis.differentiation.GradientField;
27 import org.hipparchus.dfp.Dfp;
28 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
29 import org.hipparchus.geometry.euclidean.threed.Rotation;
30 import org.hipparchus.geometry.euclidean.threed.SphericalCoordinates;
31 import org.hipparchus.geometry.euclidean.threed.Vector3D;
32 import org.hipparchus.linear.Array2DRowRealMatrix;
33 import org.hipparchus.linear.MatrixUtils;
34 import org.hipparchus.linear.RealMatrix;
35 import org.hipparchus.ode.AbstractIntegrator;
36 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
37 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
38 import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator;
39 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
40 import org.hipparchus.util.FastMath;
41 import org.junit.jupiter.api.AfterEach;
42 import org.junit.jupiter.api.Assertions;
43 import org.junit.jupiter.api.BeforeEach;
44 import org.junit.jupiter.api.DisplayName;
45 import org.junit.jupiter.api.Test;
46 import org.mockito.Mockito;
47 import org.orekit.Utils;
48 import org.orekit.attitudes.Attitude;
49 import org.orekit.attitudes.LofOffset;
50 import org.orekit.bodies.CelestialBodyFactory;
51 import org.orekit.errors.OrekitException;
52 import org.orekit.forces.AbstractLegacyForceModelTest;
53 import org.orekit.forces.ForceModel;
54 import org.orekit.forces.gravity.potential.GRGSFormatReader;
55 import org.orekit.forces.gravity.potential.GravityFieldFactory;
56 import org.orekit.forces.gravity.potential.ICGEMFormatReader;
57 import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
58 import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider.NormalizedSphericalHarmonics;
59 import org.orekit.forces.gravity.potential.TideSystem;
60 import org.orekit.frames.Frame;
61 import org.orekit.frames.FramesFactory;
62 import org.orekit.frames.LOFType;
63 import org.orekit.frames.StaticTransform;
64 import org.orekit.frames.Transform;
65 import org.orekit.orbits.CartesianOrbit;
66 import org.orekit.orbits.EquinoctialOrbit;
67 import org.orekit.orbits.FieldKeplerianOrbit;
68 import org.orekit.orbits.KeplerianOrbit;
69 import org.orekit.orbits.Orbit;
70 import org.orekit.orbits.OrbitType;
71 import org.orekit.orbits.PositionAngleType;
72 import org.orekit.propagation.*;
73 import org.orekit.propagation.analytical.EcksteinHechlerPropagator;
74 import org.orekit.propagation.numerical.FieldNumericalPropagator;
75 import org.orekit.propagation.numerical.NumericalPropagator;
76 import org.orekit.propagation.sampling.OrekitFixedStepHandler;
77 import org.orekit.time.AbsoluteDate;
78 import org.orekit.time.DateComponents;
79 import org.orekit.time.FieldAbsoluteDate;
80 import org.orekit.time.TimeComponents;
81 import org.orekit.time.TimeScalesFactory;
82 import org.orekit.utils.Constants;
83 import org.orekit.utils.FieldPVCoordinates;
84 import org.orekit.utils.IERSConventions;
85 import org.orekit.utils.PVCoordinates;
86 import org.orekit.utils.PVCoordinatesProvider;
87
88
89 public class HolmesFeatherstoneAttractionModelTest extends AbstractLegacyForceModelTest {
90
91 private static final int SCALING = 930;
92
93 @Override
94 protected FieldVector3D<DerivativeStructure> accelerationDerivatives(final ForceModel forceModel,
95 final FieldSpacecraftState<DerivativeStructure> state) {
96 try {
97 final AbsoluteDate date = state.getDate().toAbsoluteDate();
98 final FieldVector3D<DerivativeStructure> position = state.getPVCoordinates().getPosition();
99 java.lang.reflect.Field bodyFrameField = HolmesFeatherstoneAttractionModel.class.getDeclaredField("bodyFrame");
100 bodyFrameField.setAccessible(true);
101 Frame bodyFrame = (Frame) bodyFrameField.get(forceModel);
102
103
104 final StaticTransform fromBodyFrame = bodyFrame.getStaticTransformTo(state.getFrame(), date);
105 final StaticTransform toBodyFrame = fromBodyFrame.getInverse();
106 final Vector3D positionBody = toBodyFrame.transformPosition(position.toVector3D());
107
108
109 final GradientHessian gh = gradientHessian((HolmesFeatherstoneAttractionModel) forceModel,
110 date, positionBody);
111
112
113 final double[] gInertial = fromBodyFrame.transformVector(new Vector3D(gh.getGradient())).toArray();
114
115
116 final RealMatrix hBody = new Array2DRowRealMatrix(gh.getHessian(), false);
117 final RealMatrix rot = new Array2DRowRealMatrix(toBodyFrame.getRotation().getMatrix());
118 final RealMatrix hInertial = rot.transpose().multiply(hBody).multiply(rot);
119
120
121 final double[] derivatives = new double[1 + state.getMass().getFreeParameters()];
122 final DerivativeStructure[] accDer = new DerivativeStructure[3];
123 for (int i = 0; i < 3; ++i) {
124
125
126 derivatives[0] = gInertial[i];
127
128
129 derivatives[1] = hInertial.getEntry(i, 0);
130 derivatives[2] = hInertial.getEntry(i, 1);
131 derivatives[3] = hInertial.getEntry(i, 2);
132
133
134
135 accDer[i] = state.getMass().getFactory().build(derivatives);
136
137 }
138
139 return new FieldVector3D<>(accDer);
140
141 } catch (IllegalArgumentException | IllegalAccessException | NoSuchFieldException | SecurityException e) {
142 return null;
143 }
144 }
145
146 @Override
147 protected FieldVector3D<Gradient> accelerationDerivativesGradient(final ForceModel forceModel,
148 final FieldSpacecraftState<Gradient> state) {
149 try {
150 final AbsoluteDate date = state.getDate().toAbsoluteDate();
151 final FieldVector3D<Gradient> position = state.getPVCoordinates().getPosition();
152 java.lang.reflect.Field bodyFrameField = HolmesFeatherstoneAttractionModel.class.getDeclaredField("bodyFrame");
153 bodyFrameField.setAccessible(true);
154 Frame bodyFrame = (Frame) bodyFrameField.get(forceModel);
155
156
157 final Transform fromBodyFrame = bodyFrame.getTransformTo(state.getFrame(), date);
158 final Transform toBodyFrame = fromBodyFrame.getInverse();
159 final Vector3D positionBody = toBodyFrame.transformPosition(position.toVector3D());
160
161
162 final GradientHessian gh = gradientHessian((HolmesFeatherstoneAttractionModel) forceModel,
163 date, positionBody);
164
165
166 final double[] gInertial = fromBodyFrame.transformVector(new Vector3D(gh.getGradient())).toArray();
167
168
169 final RealMatrix hBody = new Array2DRowRealMatrix(gh.getHessian(), false);
170 final RealMatrix rot = new Array2DRowRealMatrix(toBodyFrame.getRotation().getMatrix());
171 final RealMatrix hInertial = rot.transpose().multiply(hBody).multiply(rot);
172
173
174 final double[] derivatives = new double[state.getMass().getFreeParameters()];
175 final Gradient[] accDer = new Gradient[3];
176 for (int i = 0; i < 3; ++i) {
177
178
179 derivatives[0] = hInertial.getEntry(i, 0);
180 derivatives[1] = hInertial.getEntry(i, 1);
181 derivatives[2] = hInertial.getEntry(i, 2);
182
183
184
185 accDer[i] = new Gradient(gInertial[i], derivatives);
186
187 }
188
189 return new FieldVector3D<>(accDer);
190
191 } catch (IllegalArgumentException | IllegalAccessException | NoSuchFieldException | SecurityException e) {
192 return null;
193 }
194 }
195
196 private GradientHessian gradientHessian(final HolmesFeatherstoneAttractionModel hfModel,
197 final AbsoluteDate date, final Vector3D position)
198 {
199 try {
200
201 java.lang.reflect.Field providerField = HolmesFeatherstoneAttractionModel.class.getDeclaredField("provider");
202 providerField.setAccessible(true);
203 NormalizedSphericalHarmonicsProvider provider = (NormalizedSphericalHarmonicsProvider) providerField.get(hfModel);
204 java.lang.reflect.Method createDistancePowersArrayMethod =
205 HolmesFeatherstoneAttractionModel.class.getDeclaredMethod("createDistancePowersArray", Double.TYPE);
206 createDistancePowersArrayMethod.setAccessible(true);
207 java.lang.reflect.Method createCosSinArraysMethod =
208 HolmesFeatherstoneAttractionModel.class.getDeclaredMethod("createCosSinArrays", Double.TYPE, Double.TYPE);
209 createCosSinArraysMethod.setAccessible(true);
210 java.lang.reflect.Method computeTesseralMethod =
211 HolmesFeatherstoneAttractionModel.class.getDeclaredMethod("computeTesseral",
212 Integer.TYPE, Integer.TYPE, Integer.TYPE,
213 Double.TYPE, Double.TYPE, Double.TYPE,
214 double[].class, double[].class, double[].class,
215 double[].class, double[].class, double[].class);
216 computeTesseralMethod.setAccessible(true);
217
218 final int degree = provider.getMaxDegree();
219 final int order = provider.getMaxOrder();
220 final NormalizedSphericalHarmonics harmonics = provider.onDate(date);
221
222
223 double[] pnm0Plus2 = new double[degree + 1];
224 double[] pnm0Plus1 = new double[degree + 1];
225 double[] pnm0 = new double[degree + 1];
226 double[] pnm1Plus1 = new double[degree + 1];
227 double[] pnm1 = new double[degree + 1];
228 final double[] pnm2 = new double[degree + 1];
229
230
231 final double x = position.getX();
232 final double y = position.getY();
233 final double z = position.getZ();
234 final double x2 = x * x;
235 final double y2 = y * y;
236 final double z2 = z * z;
237 final double r2 = x2 + y2 + z2;
238 final double r = FastMath.sqrt (r2);
239 final double rho2 = x2 + y2;
240 final double rho = FastMath.sqrt(rho2);
241 final double t = z / r;
242 final double u = rho / r;
243 final double tOu = z / rho;
244
245
246 final double[] aOrN = (double[]) createDistancePowersArrayMethod.invoke(hfModel, provider.getAe() / r);
247
248
249 final double[][] cosSinLambda = (double[][]) createCosSinArraysMethod.invoke(hfModel, position.getX() / rho, position.getY() / rho);
250
251
252 int index = 0;
253 double value = 0;
254 final double[] gradient = new double[3];
255 final double[][] hessian = new double[3][3];
256 for (int m = degree; m >= 0; --m) {
257
258
259 index = (Integer) computeTesseralMethod.invoke(hfModel, m, degree, index, t, u, tOu,
260 pnm0Plus2, pnm0Plus1, pnm1Plus1, pnm0, pnm1, pnm2);
261
262 if (m <= order) {
263
264
265
266 double sumDegreeS = 0;
267 double sumDegreeC = 0;
268 double dSumDegreeSdR = 0;
269 double dSumDegreeCdR = 0;
270 double dSumDegreeSdTheta = 0;
271 double dSumDegreeCdTheta = 0;
272 double d2SumDegreeSdRdR = 0;
273 double d2SumDegreeSdRdTheta = 0;
274 double d2SumDegreeSdThetadTheta = 0;
275 double d2SumDegreeCdRdR = 0;
276 double d2SumDegreeCdRdTheta = 0;
277 double d2SumDegreeCdThetadTheta = 0;
278 for (int n = FastMath.max(2, m); n <= degree; ++n) {
279 final double qSnm = aOrN[n] * harmonics.getNormalizedSnm(n, m);
280 final double qCnm = aOrN[n] * harmonics.getNormalizedCnm(n, m);
281 final double nOr = n / r;
282 final double nnP1Or2 = nOr * (n + 1) / r;
283 final double s0 = pnm0[n] * qSnm;
284 final double c0 = pnm0[n] * qCnm;
285 final double s1 = pnm1[n] * qSnm;
286 final double c1 = pnm1[n] * qCnm;
287 final double s2 = pnm2[n] * qSnm;
288 final double c2 = pnm2[n] * qCnm;
289 sumDegreeS += s0;
290 sumDegreeC += c0;
291 dSumDegreeSdR -= nOr * s0;
292 dSumDegreeCdR -= nOr * c0;
293 dSumDegreeSdTheta += s1;
294 dSumDegreeCdTheta += c1;
295 d2SumDegreeSdRdR += nnP1Or2 * s0;
296 d2SumDegreeSdRdTheta -= nOr * s1;
297 d2SumDegreeSdThetadTheta += s2;
298 d2SumDegreeCdRdR += nnP1Or2 * c0;
299 d2SumDegreeCdRdTheta -= nOr * c1;
300 d2SumDegreeCdThetadTheta += c2;
301 }
302
303
304 final double sML = cosSinLambda[1][m];
305 final double cML = cosSinLambda[0][m];
306 value = value * u + sML * sumDegreeS + cML * sumDegreeC;
307 gradient[0] = gradient[0] * u + sML * dSumDegreeSdR + cML * dSumDegreeCdR;
308 gradient[1] = gradient[1] * u + m * (cML * sumDegreeS - sML * sumDegreeC);
309 gradient[2] = gradient[2] * u + sML * dSumDegreeSdTheta + cML * dSumDegreeCdTheta;
310 hessian[0][0] = hessian[0][0] * u + sML * d2SumDegreeSdRdR + cML * d2SumDegreeCdRdR;
311 hessian[1][0] = hessian[1][0] * u + m * (cML * dSumDegreeSdR - sML * dSumDegreeCdR);
312 hessian[2][0] = hessian[2][0] * u + sML * d2SumDegreeSdRdTheta + cML * d2SumDegreeCdRdTheta;
313 hessian[1][1] = hessian[1][1] * u - m * m * (sML * sumDegreeS + cML * sumDegreeC);
314 hessian[2][1] = hessian[2][1] * u + m * (cML * dSumDegreeSdTheta - sML * dSumDegreeCdTheta);
315 hessian[2][2] = hessian[2][2] * u + sML * d2SumDegreeSdThetadTheta + cML * d2SumDegreeCdThetadTheta;
316
317 }
318
319
320 final double[] tmp0 = pnm0Plus2;
321 pnm0Plus2 = pnm0Plus1;
322 pnm0Plus1 = pnm0;
323 pnm0 = tmp0;
324 final double[] tmp1 = pnm1Plus1;
325 pnm1Plus1 = pnm1;
326 pnm1 = tmp1;
327
328 }
329
330
331 value = FastMath.scalb(value, SCALING);
332 for (int i = 0; i < 3; ++i) {
333 gradient[i] = FastMath.scalb(gradient[i], SCALING);
334 for (int j = 0; j <= i; ++j) {
335 hessian[i][j] = FastMath.scalb(hessian[i][j], SCALING);
336 }
337 }
338
339
340
341 final double muOr = provider.getMu() / r;
342 value *= muOr;
343 gradient[0] = muOr * gradient[0] - value / r;
344 gradient[1] *= muOr;
345 gradient[2] *= muOr;
346 hessian[0][0] = muOr * hessian[0][0] - 2 * gradient[0] / r;
347 hessian[1][0] = muOr * hessian[1][0] - gradient[1] / r;
348 hessian[2][0] = muOr * hessian[2][0] - gradient[2] / r;
349 hessian[1][1] *= muOr;
350 hessian[2][1] *= muOr;
351 hessian[2][2] *= muOr;
352
353
354 final SphericalCoordinates sc = new SphericalCoordinates(position);
355 return new GradientHessian(sc.toCartesianGradient(gradient),
356 sc.toCartesianHessian(hessian, gradient));
357
358
359 } catch (IllegalArgumentException | IllegalAccessException | NoSuchFieldException |
360 SecurityException | InvocationTargetException | NoSuchMethodException e) {
361 return null;
362 }
363 }
364
365
366 private static class GradientHessian implements Serializable {
367
368
369 private static final long serialVersionUID = 20130219L;
370
371
372 private final double[] gradient;
373
374
375 private final double[][] hessian;
376
377
378
379
380
381
382
383
384 public GradientHessian(final double[] gradient, final double[][] hessian) {
385 this.gradient = gradient;
386 this.hessian = hessian;
387 }
388
389
390
391
392 public double[] getGradient() {
393 return gradient;
394 }
395
396
397
398
399 public double[][] getHessian() {
400 return hessian;
401 }
402
403 }
404
405 @Test
406 void testRelativeNumericPrecision() {
407
408
409
410
411 int max = 50;
412 NormalizedSphericalHarmonicsProvider provider = new GleasonProvider(max, max);
413 HolmesFeatherstoneAttractionModel model =
414 new HolmesFeatherstoneAttractionModel(itrf, provider);
415 Assertions.assertTrue(model.dependsOnPositionOnly());
416
417
418
419
420
421
422
423
424
425
426 int digits = 30;
427 ReferenceFieldModel refModel = new ReferenceFieldModel(provider, digits);
428
429 double r = 1.25;
430 for (double theta = 0.01; theta < 3.14; theta += 0.1) {
431 Vector3D position = new Vector3D(r * FastMath.sin(theta), 0.0, r * FastMath.cos(theta));
432 Dfp refValue = refModel.nonCentralPart(null, position);
433 double value = model.nonCentralPart(null, position, model.getMu());
434 double relativeError = error(refValue, value).divide(refValue).toDouble();
435 Assertions.assertEquals(0, relativeError, 7.0e-15);
436 }
437
438 }
439
440 @Test
441 void testValue() {
442
443 int max = 50;
444 NormalizedSphericalHarmonicsProvider provider = new GleasonProvider(max, max);
445 HolmesFeatherstoneAttractionModel model =
446 new HolmesFeatherstoneAttractionModel(itrf, provider);
447
448 double r = 1.25;
449 for (double lambda = 0; lambda < 2 * FastMath.PI; lambda += 0.5) {
450 for (double theta = 0.05; theta < 3.11; theta += 0.03) {
451 Vector3D position = new Vector3D(r * FastMath.sin(theta) * FastMath.cos(lambda),
452 r * FastMath.sin(theta) * FastMath.sin(lambda),
453 r * FastMath.cos(theta));
454 double refValue = provider.getMu() / position.getNorm() +
455 model.nonCentralPart(AbsoluteDate.GPS_EPOCH, position, model.getMu());
456 double value = model.value(AbsoluteDate.GPS_EPOCH, position, model.getMu());
457 Assertions.assertEquals(refValue, value, 1.0e-15 * FastMath.abs(refValue));
458 }
459 }
460
461 }
462
463
464
465
466
467
468 @Test
469 public void RealFieldTest() {
470 DSFactory factory = new DSFactory(6, 4);
471 DerivativeStructure a_0 = factory.variable(0, 7201009.7124401);
472 DerivativeStructure e_0 = factory.variable(1, 5e-3);
473 DerivativeStructure i_0 = factory.variable(2, 98.7 * FastMath.PI / 180);
474 DerivativeStructure R_0 = factory.variable(3, 15.0 * 22.5 * FastMath.PI / 180);
475 DerivativeStructure O_0 = factory.variable(4, 93.0 * FastMath.PI / 180);
476 DerivativeStructure n_0 = factory.variable(5, 0.1);
477
478 Field<DerivativeStructure> field = a_0.getField();
479 DerivativeStructure zero = field.getZero();
480
481 FieldAbsoluteDate<DerivativeStructure> J2000 = new FieldAbsoluteDate<>(field);
482
483 Frame EME = FramesFactory.getEME2000();
484
485 FieldKeplerianOrbit<DerivativeStructure> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
486 PositionAngleType.MEAN,
487 EME,
488 J2000,
489 zero.add(Constants.EIGEN5C_EARTH_MU));
490
491 FieldSpacecraftState<DerivativeStructure> initialState = new FieldSpacecraftState<>(FKO);
492
493 SpacecraftState iSR = initialState.toSpacecraftState();
494 OrbitType type = OrbitType.EQUINOCTIAL;
495 double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(10.).getTolerances(FKO.toOrbit(), type);
496
497
498 AdaptiveStepsizeFieldIntegrator<DerivativeStructure> integrator =
499 new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
500 integrator.setInitialStepSize(60);
501 AdaptiveStepsizeIntegrator RIntegrator =
502 new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
503 RIntegrator.setInitialStepSize(60);
504
505 FieldNumericalPropagator<DerivativeStructure> FNP = new FieldNumericalPropagator<>(field, integrator);
506 FNP.setOrbitType(type);
507 FNP.setPositionAngleType(PositionAngleType.TRUE);
508 FNP.setInitialState(initialState);
509
510 NumericalPropagator NP = new NumericalPropagator(RIntegrator);
511 NP.setOrbitType(FNP.getOrbitType());
512 NP.setPositionAngleType(FNP.getPositionAngleType());
513 NP.setInitialState(iSR);
514
515 double[][] c = new double[3][1];
516 c[0][0] = 0.0;
517 c[2][0] = normalizedC20;
518 double[][] s = new double[3][1];
519 NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(6378136.460, mu,
520 TideSystem.UNKNOWN,
521 c, s);
522 HolmesFeatherstoneAttractionModel forceModel =
523 new HolmesFeatherstoneAttractionModel(itrf, provider);
524
525 FNP.addForceModel(forceModel);
526 NP.addForceModel(forceModel);
527
528
529 checkRealFieldPropagation(FKO, PositionAngleType.MEAN, 1005., NP, FNP,
530 1.0e-14, 6.0e-8, 1.8e-11, 1.8e-10,
531 1, false);
532 }
533
534
535
536
537
538 @Test
539 public void RealFieldExpectErrorTest() {
540 DSFactory factory = new DSFactory(6, 0);
541 DerivativeStructure a_0 = factory.variable(0, 7201009.7124401);
542 DerivativeStructure e_0 = factory.variable(1, 1e-3);
543 DerivativeStructure i_0 = factory.variable(2, 98.7 * FastMath.PI / 180);
544 DerivativeStructure R_0 = factory.variable(3, 15.0 * 22.5 * FastMath.PI / 180);
545 DerivativeStructure O_0 = factory.variable(4, 93.0 * FastMath.PI / 180);
546 DerivativeStructure n_0 = factory.variable(5, 0.1);
547
548 Field<DerivativeStructure> field = a_0.getField();
549 DerivativeStructure zero = field.getZero();
550
551 FieldAbsoluteDate<DerivativeStructure> J2000 = new FieldAbsoluteDate<>(field);
552
553 Frame EME = FramesFactory.getEME2000();
554 FieldKeplerianOrbit<DerivativeStructure> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
555 PositionAngleType.MEAN,
556 EME,
557 J2000,
558 zero.add(Constants.EIGEN5C_EARTH_MU));
559
560 FieldSpacecraftState<DerivativeStructure> initialState = new FieldSpacecraftState<>(FKO);
561
562 SpacecraftState iSR = initialState.toSpacecraftState();
563 OrbitType type = OrbitType.EQUINOCTIAL;
564 double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(10.).getTolerances(FKO.toOrbit(), type);
565
566
567 AdaptiveStepsizeFieldIntegrator<DerivativeStructure> integrator =
568 new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
569 integrator.setInitialStepSize(60);
570 AdaptiveStepsizeIntegrator RIntegrator =
571 new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
572 RIntegrator.setInitialStepSize(60);
573
574 FieldNumericalPropagator<DerivativeStructure> FNP = new FieldNumericalPropagator<>(field, integrator);
575 FNP.setOrbitType(type);
576 FNP.setInitialState(initialState);
577
578 NumericalPropagator NP = new NumericalPropagator(RIntegrator);
579 NP.setOrbitType(type);
580 NP.setInitialState(iSR);
581
582 double[][] c = new double[3][1];
583 c[0][0] = 0.0;
584 c[2][0] = normalizedC20;
585 double[][] s = new double[3][1];
586 NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(6378136.460, mu,
587 TideSystem.UNKNOWN,
588 c, s);
589 HolmesFeatherstoneAttractionModel forceModel =
590 new HolmesFeatherstoneAttractionModel(itrf, provider);
591
592
593 NP.addForceModel(forceModel);
594
595 FieldAbsoluteDate<DerivativeStructure> target = J2000.shiftedBy(100.);
596 FieldSpacecraftState<DerivativeStructure> finalState_DS = FNP.propagate(target);
597 SpacecraftState finalState_R = NP.propagate(target.toAbsoluteDate());
598 FieldPVCoordinates<DerivativeStructure> finPVC_DS = finalState_DS.getPVCoordinates();
599 PVCoordinates finPVC_R = finalState_R.getPVCoordinates();
600 Assertions.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getX() - finPVC_R.getPosition().getX()) < FastMath.abs(finPVC_R.getPosition().getX()) * 1e-11);
601 Assertions.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getY() - finPVC_R.getPosition().getY()) < FastMath.abs(finPVC_R.getPosition().getY()) * 1e-11);
602 Assertions.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getZ() - finPVC_R.getPosition().getZ()) < FastMath.abs(finPVC_R.getPosition().getZ()) * 1e-11);
603 }
604 @Test
605 void testGradient() {
606
607 int max = 50;
608 NormalizedSphericalHarmonicsProvider provider = new GleasonProvider(max, max);
609 HolmesFeatherstoneAttractionModel model =
610 new HolmesFeatherstoneAttractionModel(itrf, provider);
611
612 double r = 1.25;
613 for (double lambda = 0; lambda < 2 * FastMath.PI; lambda += 0.5) {
614 for (double theta = 0.05; theta < 3.11; theta += 0.03) {
615 Vector3D position = new Vector3D(r * FastMath.sin(theta) * FastMath.cos(lambda),
616 r * FastMath.sin(theta) * FastMath.sin(lambda),
617 r * FastMath.cos(theta));
618 double[] refGradient = gradient(model, null, position, 1.0e-3);
619 double norm = FastMath.sqrt(refGradient[0] * refGradient[0] +
620 refGradient[1] * refGradient[1] +
621 refGradient[2] * refGradient[2]);
622 double[] gradient = model.gradient(null, position, model.getMu());
623 double errorX = refGradient[0] - gradient[0];
624 double errorY = refGradient[1] - gradient[1];
625 double errorZ = refGradient[2] - gradient[2];
626 double relativeError = FastMath.sqrt(errorX * errorX + errorY * errorY + errorZ * errorZ) /
627 norm;
628 Assertions.assertEquals(0, relativeError, 3.0e-12);
629 }
630 }
631
632 }
633
634 @Test
635 void testHessian() {
636
637 int max = 50;
638 NormalizedSphericalHarmonicsProvider provider = new GleasonProvider(max, max);
639 HolmesFeatherstoneAttractionModel model =
640 new HolmesFeatherstoneAttractionModel(itrf, provider);
641
642 double r = 1.25;
643 for (double lambda = 0; lambda < 2 * FastMath.PI; lambda += 0.5) {
644 for (double theta = 0.05; theta < 3.11; theta += 0.03) {
645 Vector3D position = new Vector3D(r * FastMath.sin(theta) * FastMath.cos(lambda),
646 r * FastMath.sin(theta) * FastMath.sin(lambda),
647 r * FastMath.cos(theta));
648 double[][] refHessian = hessian(model, null, position, 1.0e-3);
649 double[][] hessian = gradientHessian(model, null, position).getHessian();
650 double normH2 = 0;
651 double normE2 = 0;
652 for (int i = 0; i < 3; ++i) {
653 for (int j = 0; j < 3; ++j) {
654 double error = refHessian[i][j] - hessian[i][j];
655 normH2 += refHessian[i][j] * refHessian[i][j];
656 normE2 += error * error;
657 }
658 }
659 Assertions.assertEquals(0, FastMath.sqrt(normE2 / normH2), 5.0e-12);
660 }
661 }
662
663 }
664
665 private Dfp error(Dfp refValue, double value) {
666 return refValue.getField().newDfp(value).subtract(refValue);
667 }
668
669 private double[] gradient(final HolmesFeatherstoneAttractionModel model,
670 final AbsoluteDate date, final Vector3D position, final double h)
671 {
672 return new double[] {
673 differential8(model.nonCentralPart(date, position.add(new Vector3D(-4 * h, Vector3D.PLUS_I)), model.getMu()),
674 model.nonCentralPart(date, position.add(new Vector3D(-3 * h, Vector3D.PLUS_I)), model.getMu()),
675 model.nonCentralPart(date, position.add(new Vector3D(-2 * h, Vector3D.PLUS_I)), model.getMu()),
676 model.nonCentralPart(date, position.add(new Vector3D(-1 * h, Vector3D.PLUS_I)), model.getMu()),
677 model.nonCentralPart(date, position.add(new Vector3D(+1 * h, Vector3D.PLUS_I)), model.getMu()),
678 model.nonCentralPart(date, position.add(new Vector3D(+2 * h, Vector3D.PLUS_I)), model.getMu()),
679 model.nonCentralPart(date, position.add(new Vector3D(+3 * h, Vector3D.PLUS_I)), model.getMu()),
680 model.nonCentralPart(date, position.add(new Vector3D(+4 * h, Vector3D.PLUS_I)), model.getMu()),
681 h),
682 differential8(model.nonCentralPart(date, position.add(new Vector3D(-4 * h, Vector3D.PLUS_J)), model.getMu()),
683 model.nonCentralPart(date, position.add(new Vector3D(-3 * h, Vector3D.PLUS_J)), model.getMu()),
684 model.nonCentralPart(date, position.add(new Vector3D(-2 * h, Vector3D.PLUS_J)), model.getMu()),
685 model.nonCentralPart(date, position.add(new Vector3D(-1 * h, Vector3D.PLUS_J)), model.getMu()),
686 model.nonCentralPart(date, position.add(new Vector3D(+1 * h, Vector3D.PLUS_J)), model.getMu()),
687 model.nonCentralPart(date, position.add(new Vector3D(+2 * h, Vector3D.PLUS_J)), model.getMu()),
688 model.nonCentralPart(date, position.add(new Vector3D(+3 * h, Vector3D.PLUS_J)), model.getMu()),
689 model.nonCentralPart(date, position.add(new Vector3D(+4 * h, Vector3D.PLUS_J)), model.getMu()),
690 h),
691 differential8(model.nonCentralPart(date, position.add(new Vector3D(-4 * h, Vector3D.PLUS_K)), model.getMu()),
692 model.nonCentralPart(date, position.add(new Vector3D(-3 * h, Vector3D.PLUS_K)), model.getMu()),
693 model.nonCentralPart(date, position.add(new Vector3D(-2 * h, Vector3D.PLUS_K)), model.getMu()),
694 model.nonCentralPart(date, position.add(new Vector3D(-1 * h, Vector3D.PLUS_K)), model.getMu()),
695 model.nonCentralPart(date, position.add(new Vector3D(+1 * h, Vector3D.PLUS_K)), model.getMu()),
696 model.nonCentralPart(date, position.add(new Vector3D(+2 * h, Vector3D.PLUS_K)), model.getMu()),
697 model.nonCentralPart(date, position.add(new Vector3D(+3 * h, Vector3D.PLUS_K)), model.getMu()),
698 model.nonCentralPart(date, position.add(new Vector3D(+4 * h, Vector3D.PLUS_K)), model.getMu()),
699 h)
700 };
701 }
702
703 private double[][] hessian(final HolmesFeatherstoneAttractionModel model,
704 final AbsoluteDate date, final Vector3D position, final double h)
705 {
706 return new double[][] {
707 differential8(model.gradient(date, position.add(new Vector3D(-4 * h, Vector3D.PLUS_I)), model.getMu()),
708 model.gradient(date, position.add(new Vector3D(-3 * h, Vector3D.PLUS_I)), model.getMu()),
709 model.gradient(date, position.add(new Vector3D(-2 * h, Vector3D.PLUS_I)), model.getMu()),
710 model.gradient(date, position.add(new Vector3D(-1 * h, Vector3D.PLUS_I)), model.getMu()),
711 model.gradient(date, position.add(new Vector3D(+1 * h, Vector3D.PLUS_I)), model.getMu()),
712 model.gradient(date, position.add(new Vector3D(+2 * h, Vector3D.PLUS_I)), model.getMu()),
713 model.gradient(date, position.add(new Vector3D(+3 * h, Vector3D.PLUS_I)), model.getMu()),
714 model.gradient(date, position.add(new Vector3D(+4 * h, Vector3D.PLUS_I)), model.getMu()),
715 h),
716 differential8(model.gradient(date, position.add(new Vector3D(-4 * h, Vector3D.PLUS_J)), model.getMu()),
717 model.gradient(date, position.add(new Vector3D(-3 * h, Vector3D.PLUS_J)), model.getMu()),
718 model.gradient(date, position.add(new Vector3D(-2 * h, Vector3D.PLUS_J)), model.getMu()),
719 model.gradient(date, position.add(new Vector3D(-1 * h, Vector3D.PLUS_J)), model.getMu()),
720 model.gradient(date, position.add(new Vector3D(+1 * h, Vector3D.PLUS_J)), model.getMu()),
721 model.gradient(date, position.add(new Vector3D(+2 * h, Vector3D.PLUS_J)), model.getMu()),
722 model.gradient(date, position.add(new Vector3D(+3 * h, Vector3D.PLUS_J)), model.getMu()),
723 model.gradient(date, position.add(new Vector3D(+4 * h, Vector3D.PLUS_J)), model.getMu()),
724 h),
725 differential8(model.gradient(date, position.add(new Vector3D(-4 * h, Vector3D.PLUS_K)), model.getMu()),
726 model.gradient(date, position.add(new Vector3D(-3 * h, Vector3D.PLUS_K)), model.getMu()),
727 model.gradient(date, position.add(new Vector3D(-2 * h, Vector3D.PLUS_K)), model.getMu()),
728 model.gradient(date, position.add(new Vector3D(-1 * h, Vector3D.PLUS_K)), model.getMu()),
729 model.gradient(date, position.add(new Vector3D(+1 * h, Vector3D.PLUS_K)), model.getMu()),
730 model.gradient(date, position.add(new Vector3D(+2 * h, Vector3D.PLUS_K)), model.getMu()),
731 model.gradient(date, position.add(new Vector3D(+3 * h, Vector3D.PLUS_K)), model.getMu()),
732 model.gradient(date, position.add(new Vector3D(+4 * h, Vector3D.PLUS_K)), model.getMu()),
733 h)
734 };
735 }
736
737
738
739
740
741
742
743 private static class GleasonProvider implements NormalizedSphericalHarmonicsProvider {
744
745 private final int degree;
746 private final int order;
747
748 public GleasonProvider(int degree, int order) {
749 this.degree = degree;
750 this.order = order;
751 }
752
753 public int getMaxDegree() {
754 return degree;
755 }
756
757 public int getMaxOrder() {
758 return order;
759 }
760
761 public double getMu() {
762 return 1;
763 }
764
765 public double getAe() {
766 return 1;
767 }
768
769 public AbsoluteDate getReferenceDate() {
770 return null;
771 }
772
773 public TideSystem getTideSystem() {
774 return TideSystem.UNKNOWN;
775 }
776
777 @Override
778 public NormalizedSphericalHarmonics onDate(final AbsoluteDate date) {
779 return new NormalizedSphericalHarmonics() {
780 @Override
781 public double getNormalizedCnm(int n, int m) {
782 return 1;
783 }
784
785 @Override
786 public double getNormalizedSnm(int n, int m) {
787 return 1;
788 }
789
790 @Override
791 public AbsoluteDate getDate() {
792 return date;
793 }
794 };
795 }
796
797 }
798
799
800 @Test
801 void testHelioSynchronous()
802 {
803
804
805 AbsoluteDate date = new AbsoluteDate(new DateComponents(1970, 7, 1),
806 new TimeComponents(13, 59, 27.816),
807 TimeScalesFactory.getUTC());
808 Transform itrfToEME2000 = itrf.getTransformTo(FramesFactory.getEME2000(), date);
809 Vector3D pole = itrfToEME2000.transformVector(Vector3D.PLUS_K);
810 Frame poleAligned = new Frame(FramesFactory.getEME2000(),
811 new Transform(date, new Rotation(pole, Vector3D.PLUS_K)),
812 "pole aligned", true);
813
814 double i = FastMath.toRadians(98.7);
815 double omega = FastMath.toRadians(93.0);
816 double OMEGA = FastMath.toRadians(15.0 * 22.5);
817 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
818 0, PositionAngleType.MEAN, poleAligned, date, mu);
819 double[][] c = new double[3][1];
820 c[0][0] = 0.0;
821 c[2][0] = normalizedC20;
822 double[][] s = new double[3][1];
823 propagator.addForceModel(new HolmesFeatherstoneAttractionModel(itrf,
824 GravityFieldFactory.getNormalizedProvider(6378136.460, mu,
825 TideSystem.UNKNOWN,
826 c, s)));
827
828
829 propagator.setStepHandler(Constants.JULIAN_DAY, new SpotStepHandler());
830 propagator.setInitialState(new SpacecraftState(orbit));
831 propagator.propagate(date.shiftedBy(7 * Constants.JULIAN_DAY));
832 Assertions.assertTrue(propagator.getCalls() < 9200);
833
834 }
835
836 private static class SpotStepHandler implements OrekitFixedStepHandler {
837
838 private final PVCoordinatesProvider sun;
839 private double previous;
840
841 public SpotStepHandler() {
842 sun = CelestialBodyFactory.getSun();
843 previous = Double.NaN;
844 }
845
846 public void handleStep(SpacecraftState currentState) {
847
848
849 AbsoluteDate current = currentState.getDate();
850 Vector3D sunPos = sun.getPosition(current , FramesFactory.getEME2000());
851 Vector3D normal = currentState.getPVCoordinates().getMomentum();
852 double angle = Vector3D.angle(sunPos , normal);
853 if (! Double.isNaN(previous)) {
854 Assertions.assertEquals(previous, angle, 0.0013);
855 }
856 previous = angle;
857 }
858
859 }
860
861
862 @Test
863 void testEcksteinHechlerReference() {
864
865
866 AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
867 Vector3D position = new Vector3D(3220103., 69623., 6449822.);
868 Vector3D velocity = new Vector3D(6414.7, -2006., -3180.);
869
870 Transform itrfToEME2000 = itrf.getTransformTo(FramesFactory.getEME2000(), date);
871 Vector3D pole = itrfToEME2000.transformVector(Vector3D.PLUS_K);
872 Frame poleAligned = new Frame(FramesFactory.getEME2000(),
873 new Transform(date, new Rotation(pole, Vector3D.PLUS_K)),
874 "pole aligned", true);
875
876 Orbit initialOrbit = new EquinoctialOrbit(new PVCoordinates(position, velocity),
877 poleAligned, date, mu);
878
879 propagator.addForceModel(new HolmesFeatherstoneAttractionModel(itrf,
880 GravityFieldFactory.getNormalizedProvider(ae, mu,
881 TideSystem.UNKNOWN,
882 new double[][] {
883 { 0.0 }, { 0.0 }, { normalizedC20 }, { normalizedC30 },
884 { normalizedC40 }, { normalizedC50 }, { normalizedC60 },
885 },
886 new double[][] {
887 { 0.0 }, { 0.0 }, { 0.0 }, { 0.0 },
888 { 0.0 }, { 0.0 }, { 0.0 },
889 })));
890
891
892 propagator.setInitialState(new SpacecraftState(initialOrbit));
893 propagator.setOrbitType(OrbitType.EQUINOCTIAL);
894 propagator.setPositionAngleType(PositionAngleType.TRUE);
895 propagator.setStepHandler(20, new EckStepHandler(initialOrbit, ae,
896 unnormalizedC20, unnormalizedC30, unnormalizedC40,
897 unnormalizedC50, unnormalizedC60));
898 propagator.propagate(date.shiftedBy(50000));
899 Assertions.assertTrue(propagator.getCalls() < 1100);
900
901 }
902
903 private static class EckStepHandler implements OrekitFixedStepHandler {
904
905
906 private static final double mu = 3.986004415e+14;
907
908 private final EcksteinHechlerPropagator referencePropagator;
909
910 private EckStepHandler(Orbit initialOrbit, double ae,
911 double c20, double c30, double c40, double c50, double c60) {
912 referencePropagator =
913 new EcksteinHechlerPropagator(initialOrbit,
914 ae, mu, c20, c30, c40, c50, c60);
915 }
916
917 public void handleStep(SpacecraftState currentState) {
918
919 SpacecraftState EHPOrbit = referencePropagator.propagate(currentState.getDate());
920 Vector3D posEHP = EHPOrbit.getPosition();
921 Vector3D posDROZ = currentState.getPosition();
922 Vector3D velEHP = EHPOrbit.getPVCoordinates().getVelocity();
923 Vector3D dif = posEHP.subtract(posDROZ);
924
925 Vector3D T = new Vector3D(1 / velEHP.getNorm(), velEHP);
926 Vector3D W = EHPOrbit.getPVCoordinates().getMomentum().normalize();
927 Vector3D N = Vector3D.crossProduct(W, T);
928
929 Assertions.assertTrue(dif.getNorm() < 111);
930 Assertions.assertTrue(FastMath.abs(Vector3D.dotProduct(dif, T)) < 111);
931 Assertions.assertTrue(FastMath.abs(Vector3D.dotProduct(dif, N)) < 54);
932 Assertions.assertTrue(FastMath.abs(Vector3D.dotProduct(dif, W)) < 12);
933
934 }
935
936 }
937
938 @Test
939 void testParameterDerivative() {
940
941 Utils.setDataRoot("regular-data:potential/grgs-format");
942 GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
943
944
945 final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
946 final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
947 final Frame frame = FramesFactory.getGCRF();
948 final AbsoluteDate date = new AbsoluteDate(2005, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI());
949 final CartesianOrbit orbit = new CartesianOrbit(new PVCoordinates(pos, vel), frame,
950 date, Constants.EIGEN5C_EARTH_MU);
951 final LofOffset lofOffset = new LofOffset(frame, LOFType.LVLH_CCSDS);
952 final Attitude attitude = lofOffset.getAttitude(orbit, date, frame);
953 final SpacecraftState state = new SpacecraftState(orbit, attitude);
954
955 final HolmesFeatherstoneAttractionModel holmesFeatherstoneModel =
956 new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true),
957 GravityFieldFactory.getNormalizedProvider(20, 20));
958
959 final String name = NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT;
960 checkParameterDerivative(state, holmesFeatherstoneModel, name, 1.0e-5, 5.0e-12);
961
962 }
963
964 @Test
965 void testParameterDerivativeGradient() {
966
967 Utils.setDataRoot("regular-data:potential/grgs-format");
968 GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
969
970
971 final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
972 final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
973 final SpacecraftState state =
974 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
975 FramesFactory.getGCRF(),
976 new AbsoluteDate(2005, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
977 GravityFieldFactory.getUnnormalizedProvider(1, 1).getMu()));
978
979 final HolmesFeatherstoneAttractionModel holmesFeatherstoneModel =
980 new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true),
981 GravityFieldFactory.getNormalizedProvider(20, 20));
982
983 final String name = NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT;
984 checkParameterDerivativeGradient(state, holmesFeatherstoneModel, name, 1.0e-5, 5.0e-12);
985
986 }
987
988 @Test
989 void testTimeDependentField() {
990
991 Utils.setDataRoot("regular-data:potential/icgem-format");
992 GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("eigen-6s-truncated", true));
993
994 final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
995 final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
996 final SpacecraftState spacecraftState =
997 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
998 FramesFactory.getGCRF(),
999 new AbsoluteDate(2005, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
1000 GravityFieldFactory.getUnnormalizedProvider(1, 1).getMu()));
1001
1002 double dP = 0.1;
1003 double duration = 3 * Constants.JULIAN_DAY;
1004 BoundedPropagator fixedFieldEphemeris = createEphemeris(dP, spacecraftState, duration,
1005 GravityFieldFactory.getConstantNormalizedProvider(8, 8, new AbsoluteDate("2005-01-01T00:00:00.000", TimeScalesFactory.getTAI())));
1006 BoundedPropagator varyingFieldEphemeris = createEphemeris(dP, spacecraftState, duration,
1007 GravityFieldFactory.getNormalizedProvider(8, 8));
1008
1009 double step = 60.0;
1010 double maxDeltaT = 0;
1011 double maxDeltaN = 0;
1012 double maxDeltaW = 0;
1013 for (AbsoluteDate date = fixedFieldEphemeris.getMinDate();
1014 date.compareTo(fixedFieldEphemeris.getMaxDate()) < 0;
1015 date = date.shiftedBy(step)) {
1016 PVCoordinates pvFixedField = fixedFieldEphemeris.getPVCoordinates(date, FramesFactory.getGCRF());
1017 PVCoordinates pvVaryingField = varyingFieldEphemeris.getPVCoordinates(date, FramesFactory.getGCRF());
1018 Vector3D t = pvFixedField.getVelocity().normalize();
1019 Vector3D w = pvFixedField.getMomentum().normalize();
1020 Vector3D n = Vector3D.crossProduct(w, t);
1021 Vector3D delta = pvVaryingField.getPosition().subtract(pvFixedField.getPosition());
1022 maxDeltaT = FastMath.max(maxDeltaT, FastMath.abs(Vector3D.dotProduct(delta, t)));
1023 maxDeltaN = FastMath.max(maxDeltaN, FastMath.abs(Vector3D.dotProduct(delta, n)));
1024 maxDeltaW = FastMath.max(maxDeltaW, FastMath.abs(Vector3D.dotProduct(delta, w)));
1025 }
1026 Assertions.assertTrue(maxDeltaT > 0.65);
1027 Assertions.assertTrue(maxDeltaT < 0.85);
1028 Assertions.assertTrue(maxDeltaN > 0.02);
1029 Assertions.assertTrue(maxDeltaN < 0.03);
1030 Assertions.assertTrue(maxDeltaW > 0.05);
1031 Assertions.assertTrue(maxDeltaW < 0.10);
1032
1033 }
1034
1035 private BoundedPropagator createEphemeris(double dP, SpacecraftState initialState, double duration,
1036 NormalizedSphericalHarmonicsProvider provider)
1037 {
1038 double[][] tol = ToleranceProvider.getDefaultToleranceProvider(dP).getTolerances(initialState.getOrbit(), OrbitType.CARTESIAN);
1039 AbstractIntegrator integrator =
1040 new DormandPrince853Integrator(0.001, 120.0, tol[0], tol[1]);
1041 NumericalPropagator propagator = new NumericalPropagator(integrator);
1042 final EphemerisGenerator generator = propagator.getEphemerisGenerator();
1043 propagator.setOrbitType(OrbitType.CARTESIAN);
1044 propagator.addForceModel(new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider));
1045 propagator.setInitialState(initialState);
1046 propagator.propagate(initialState.getDate().shiftedBy(duration));
1047 return generator.getGeneratedEphemeris();
1048 }
1049
1050 @Test
1051 void testStateJacobian()
1052 {
1053
1054 Utils.setDataRoot("regular-data:potential/grgs-format");
1055 GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
1056
1057
1058 AbsoluteDate date = new AbsoluteDate(new DateComponents(2000, 7, 1),
1059 new TimeComponents(13, 59, 27.816),
1060 TimeScalesFactory.getUTC());
1061 double i = FastMath.toRadians(98.7);
1062 double omega = FastMath.toRadians(93.0);
1063 double OMEGA = FastMath.toRadians(15.0 * 22.5);
1064 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
1065 0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date, mu);
1066 OrbitType integrationType = OrbitType.CARTESIAN;
1067 double[][] tolerances = ToleranceProvider.getDefaultToleranceProvider(0.01).getTolerances(orbit, integrationType);
1068
1069 propagator = new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
1070 tolerances[0], tolerances[1]));
1071 propagator.setOrbitType(integrationType);
1072 HolmesFeatherstoneAttractionModel hfModel =
1073 new HolmesFeatherstoneAttractionModel(itrf, GravityFieldFactory.getNormalizedProvider(50, 50));
1074 Assertions.assertEquals(TideSystem.UNKNOWN, hfModel.getTideSystem());
1075 propagator.addForceModel(hfModel);
1076 SpacecraftState state0 = new SpacecraftState(orbit);
1077 propagator.setInitialState(state0);
1078
1079 checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0),
1080 50000, tolerances[0], 7.8e-6);
1081 }
1082
1083 @Test
1084 void testStateJacobianVs80Implementation()
1085 {
1086
1087 Utils.setDataRoot("regular-data:potential/grgs-format");
1088 GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
1089
1090
1091 AbsoluteDate date = new AbsoluteDate(new DateComponents(2000, 7, 1),
1092 new TimeComponents(13, 59, 27.816),
1093 TimeScalesFactory.getUTC());
1094 double i = FastMath.toRadians(98.7);
1095 double omega = FastMath.toRadians(93.0);
1096 double OMEGA = FastMath.toRadians(15.0 * 22.5);
1097 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
1098 0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date, mu);
1099
1100 HolmesFeatherstoneAttractionModel hfModel =
1101 new HolmesFeatherstoneAttractionModel(itrf, GravityFieldFactory.getNormalizedProvider(50, 50));
1102 Assertions.assertEquals(TideSystem.UNKNOWN, hfModel.getTideSystem());
1103 SpacecraftState state = new SpacecraftState(orbit);
1104
1105 checkStateJacobianVs80Implementation(state, hfModel,
1106 new LofOffset(state.getFrame(), LOFType.LVLH_CCSDS),
1107 2.0e-15, false);
1108
1109 }
1110
1111 @Test
1112 void testStateJacobianVs80ImplementationGradient()
1113 {
1114
1115 Utils.setDataRoot("regular-data:potential/grgs-format");
1116 GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
1117
1118
1119 AbsoluteDate date = new AbsoluteDate(new DateComponents(2000, 7, 1),
1120 new TimeComponents(13, 59, 27.816),
1121 TimeScalesFactory.getUTC());
1122 double i = FastMath.toRadians(98.7);
1123 double omega = FastMath.toRadians(93.0);
1124 double OMEGA = FastMath.toRadians(15.0 * 22.5);
1125 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
1126 0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date, mu);
1127
1128 HolmesFeatherstoneAttractionModel hfModel =
1129 new HolmesFeatherstoneAttractionModel(itrf, GravityFieldFactory.getNormalizedProvider(50, 50));
1130 Assertions.assertEquals(TideSystem.UNKNOWN, hfModel.getTideSystem());
1131 SpacecraftState state = new SpacecraftState(orbit);
1132
1133 checkStateJacobianVs80ImplementationGradient(state, hfModel,
1134 new LofOffset(state.getFrame(), LOFType.LVLH_CCSDS),
1135 2.0e-15, false);
1136
1137 }
1138
1139 @Test
1140 void testStateJacobianVsFiniteDifferences()
1141 {
1142
1143 Utils.setDataRoot("regular-data:potential/grgs-format");
1144 GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
1145
1146
1147 AbsoluteDate date = new AbsoluteDate(new DateComponents(2000, 7, 1),
1148 new TimeComponents(13, 59, 27.816),
1149 TimeScalesFactory.getUTC());
1150 double i = FastMath.toRadians(98.7);
1151 double omega = FastMath.toRadians(93.0);
1152 double OMEGA = FastMath.toRadians(15.0 * 22.5);
1153 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
1154 0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date, mu);
1155
1156 HolmesFeatherstoneAttractionModel hfModel =
1157 new HolmesFeatherstoneAttractionModel(itrf, GravityFieldFactory.getNormalizedProvider(50, 50));
1158 Assertions.assertEquals(TideSystem.UNKNOWN, hfModel.getTideSystem());
1159 SpacecraftState state = new SpacecraftState(orbit);
1160
1161 checkStateJacobianVsFiniteDifferences(state, hfModel, Utils.defaultLaw(),
1162 10.0, 2.0e-10, false);
1163
1164 }
1165
1166 @Test
1167 void testStateJacobianVsFiniteDifferencesGradient()
1168 {
1169
1170 Utils.setDataRoot("regular-data:potential/grgs-format");
1171 GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
1172
1173
1174 AbsoluteDate date = new AbsoluteDate(new DateComponents(2000, 7, 1),
1175 new TimeComponents(13, 59, 27.816),
1176 TimeScalesFactory.getUTC());
1177 double i = FastMath.toRadians(98.7);
1178 double omega = FastMath.toRadians(93.0);
1179 double OMEGA = FastMath.toRadians(15.0 * 22.5);
1180 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
1181 0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date, mu);
1182
1183 HolmesFeatherstoneAttractionModel hfModel =
1184 new HolmesFeatherstoneAttractionModel(itrf, GravityFieldFactory.getNormalizedProvider(50, 50));
1185 Assertions.assertEquals(TideSystem.UNKNOWN, hfModel.getTideSystem());
1186 SpacecraftState state = new SpacecraftState(orbit);
1187
1188 checkStateJacobianVsFiniteDifferencesGradient(state, hfModel, Utils.defaultLaw(),
1189 10.0, 2.0e-10, false);
1190
1191 }
1192
1193 @Test
1194 void testIssue996() {
1195
1196 Utils.setDataRoot("regular-data:potential/grgs-format");
1197 GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
1198
1199
1200 AbsoluteDate date = new AbsoluteDate(new DateComponents(2000, 7, 1),
1201 new TimeComponents(13, 59, 27.816),
1202 TimeScalesFactory.getUTC());
1203 double i = FastMath.toRadians(98.7);
1204 double omega = FastMath.toRadians(93.0);
1205 double OMEGA = FastMath.toRadians(15.0 * 22.5);
1206 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
1207 0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date, mu);
1208
1209 Vector3D pos = orbit.getPosition(itrf);
1210
1211 double[] zeroGradient = new double[] {-0., 0., 0.};
1212
1213 NormalizedSphericalHarmonicsProvider provider00 = GravityFieldFactory.getNormalizedProvider(0, 0);
1214 HolmesFeatherstoneAttractionModel hfModel00 = new HolmesFeatherstoneAttractionModel(itrf, provider00);
1215
1216 Assertions.assertEquals(0., hfModel00.nonCentralPart(date, pos, mu));
1217 Assertions.assertEquals(mu / pos.getNorm(), hfModel00.value(date, pos, mu));
1218 Assertions.assertArrayEquals(zeroGradient, hfModel00.gradient(date, pos, mu));
1219
1220 NormalizedSphericalHarmonicsProvider provider10 = GravityFieldFactory.getNormalizedProvider(1, 0);
1221 HolmesFeatherstoneAttractionModel hfModel10 = new HolmesFeatherstoneAttractionModel(itrf, provider10);
1222
1223 Assertions.assertEquals(0., hfModel10.nonCentralPart(date, pos, mu));
1224 Assertions.assertEquals(mu / pos.getNorm(), hfModel10.value(date, pos, mu));
1225 Assertions.assertArrayEquals(zeroGradient, hfModel10.gradient(date, pos, mu));
1226
1227 NormalizedSphericalHarmonicsProvider provider11 = GravityFieldFactory.getNormalizedProvider(1, 1);
1228 HolmesFeatherstoneAttractionModel hfModel11 = new HolmesFeatherstoneAttractionModel(itrf, provider11);
1229
1230 Assertions.assertEquals(0., hfModel11.nonCentralPart(date, pos, mu));
1231 Assertions.assertEquals(mu / pos.getNorm(), hfModel11.value(date, pos, mu));
1232 Assertions.assertArrayEquals(zeroGradient, hfModel11.gradient(date, pos, mu));
1233
1234 }
1235
1236
1237 @Test
1238 @DisplayName("Test that acceleration derivatives with respect to absolute date are not equal to zero.")
1239 void testIssue1070() {
1240
1241
1242 final int freeParameters = 1;
1243 final GradientField field = GradientField.getField(freeParameters);
1244 final Gradient zero = field.getZero();
1245 final Gradient variable = Gradient.variable(freeParameters, 0, 0.);
1246 final FieldAbsoluteDate<Gradient> fieldAbsoluteDate = new FieldAbsoluteDate<>(field, AbsoluteDate.ARBITRARY_EPOCH).
1247 shiftedBy(variable);
1248
1249
1250 @SuppressWarnings("unchecked")
1251 final FieldSpacecraftState<Gradient> stateMock = Mockito.mock(FieldSpacecraftState.class);
1252 Mockito.when(stateMock.getDate()).thenReturn(fieldAbsoluteDate);
1253 Mockito.when(stateMock.getPosition()).thenReturn(new FieldVector3D<>(zero, zero));
1254 Mockito.when(stateMock.getFrame()).thenReturn(FramesFactory.getGCRF());
1255 Mockito.when(stateMock.getMass()).thenReturn(zero);
1256
1257
1258 final int max = 2;
1259 final NormalizedSphericalHarmonicsProvider provider = new GleasonProvider(max, max);
1260 final HolmesFeatherstoneAttractionModel model = new HolmesFeatherstoneAttractionModel(itrf, provider);
1261
1262
1263 final Gradient dummyGm = zero.add(1.);
1264 final FieldVector3D<Gradient> accelerationVector = model.acceleration(stateMock, new Gradient[] { dummyGm });
1265
1266
1267 final double[] derivatives = accelerationVector.getNormSq().getGradient();
1268 Assertions.assertNotEquals(0., MatrixUtils.createRealVector(derivatives).getNorm());
1269 }
1270
1271 @BeforeEach
1272 public void setUp() {
1273 itrf = null;
1274 propagator = null;
1275 Utils.setDataRoot("regular-data");
1276 try {
1277
1278 mu = 3.986004415e+14;
1279 ae = 6378136.460;
1280 normalizedC20 = -4.84165299820e-04;
1281 normalizedC30 = 9.57211326674e-07;
1282 normalizedC40 = 5.39990167207e-07;
1283 normalizedC50 = 6.86846073356e-08 ;
1284 normalizedC60 = -1.49953256913e-07;
1285 unnormalizedC20 = FastMath.sqrt( 5) * normalizedC20;
1286 unnormalizedC30 = FastMath.sqrt( 7) * normalizedC30;
1287 unnormalizedC40 = FastMath.sqrt( 9) * normalizedC40;
1288 unnormalizedC50 = FastMath.sqrt(11) * normalizedC50;
1289 unnormalizedC60 = FastMath.sqrt(13) * normalizedC60;
1290
1291 itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
1292 double[] absTolerance = {
1293 0.001, 1.0e-9, 1.0e-9, 1.0e-6, 1.0e-6, 1.0e-6, 0.001
1294 };
1295 double[] relTolerance = {
1296 1.0e-7, 1.0e-4, 1.0e-4, 1.0e-7, 1.0e-7, 1.0e-7, 1.0e-7
1297 };
1298 AdaptiveStepsizeIntegrator integrator =
1299 new DormandPrince853Integrator(0.001, 1000, absTolerance, relTolerance);
1300 integrator.setInitialStepSize(60);
1301 propagator = new NumericalPropagator(integrator);
1302 } catch (OrekitException oe) {
1303 Assertions.fail(oe.getMessage());
1304 }
1305 }
1306
1307 @AfterEach
1308 public void tearDown() {
1309 itrf = null;
1310 propagator = null;
1311 }
1312
1313 private double unnormalizedC20;
1314 private double unnormalizedC30;
1315 private double unnormalizedC40;
1316 private double unnormalizedC50;
1317 private double unnormalizedC60;
1318 private double normalizedC20;
1319 private double normalizedC30;
1320 private double normalizedC40;
1321 private double normalizedC50;
1322 private double normalizedC60;
1323 private double mu;
1324 private double ae;
1325
1326 private Frame itrf;
1327 private NumericalPropagator propagator;
1328
1329 }
1330
1331