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