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