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 org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.complex.Complex;
21 import org.hipparchus.complex.ComplexField;
22 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
23 import org.hipparchus.geometry.euclidean.threed.Vector3D;
24 import org.hipparchus.util.MathArrays;
25 import org.junit.jupiter.api.Assertions;
26 import org.junit.jupiter.api.BeforeEach;
27 import org.junit.jupiter.api.Test;
28 import org.orekit.Utils;
29 import org.orekit.forces.gravity.potential.GravityFieldFactory;
30 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
31 import org.orekit.frames.Frame;
32 import org.orekit.frames.FramesFactory;
33 import org.orekit.orbits.CartesianOrbit;
34 import org.orekit.propagation.FieldSpacecraftState;
35 import org.orekit.propagation.SpacecraftState;
36 import org.orekit.time.AbsoluteDate;
37 import org.orekit.time.FieldAbsoluteDate;
38 import org.orekit.time.TimeScalarFunction;
39 import org.orekit.utils.Constants;
40 import org.orekit.utils.PVCoordinates;
41
42 class J2OnlyPerturbationTest {
43
44 @BeforeEach
45 public void setUp() {
46 Utils.setDataRoot("regular-data:potential");
47 }
48
49 @Test
50 void testAcceleration() {
51
52 final UnnormalizedSphericalHarmonicsProvider unnormalizedProvider = GravityFieldFactory.getUnnormalizedProvider(2 ,0);
53 final Frame frame = FramesFactory.getGTOD(false);
54 final J2OnlyPerturbation j2OnlyPerturbation = new J2OnlyPerturbation(unnormalizedProvider, frame);
55 final SpacecraftState spacecraftState = buildState();
56
57 final Vector3D actualAcceleration = j2OnlyPerturbation.acceleration(spacecraftState, new double[0]);
58
59 final HolmesFeatherstoneAttractionModel holmesFeatherstoneAttractionModel = new HolmesFeatherstoneAttractionModel(
60 frame, GravityFieldFactory.getNormalizedProvider(unnormalizedProvider));
61 final Vector3D expectedAcceleration = holmesFeatherstoneAttractionModel.acceleration(spacecraftState,
62 new double[] { unnormalizedProvider.getMu() });
63 final double tolerance = 1e-12;
64 Assertions.assertEquals(expectedAcceleration.getX(), actualAcceleration.getX(), tolerance);
65 Assertions.assertEquals(expectedAcceleration.getY(), actualAcceleration.getY(), tolerance);
66 Assertions.assertEquals(expectedAcceleration.getZ(), actualAcceleration.getZ(), tolerance);
67 }
68
69 @Test
70 void testAccelerationField() {
71
72 final Frame frame = FramesFactory.getGTOD(false);
73 final J2OnlyPerturbation j2OnlyPerturbation = new J2OnlyPerturbation(Constants.EGM96_EARTH_MU,
74 Constants.EGM96_EARTH_EQUATORIAL_RADIUS, -Constants.EGM96_EARTH_C20, frame);
75 final SpacecraftState spacecraftState = buildState();
76 final ComplexField field = ComplexField.getInstance();
77 final FieldSpacecraftState<Complex> fieldState = new FieldSpacecraftState<>(field, spacecraftState);
78
79 final FieldVector3D<Complex> fieldAcceleration = j2OnlyPerturbation.acceleration(fieldState, null);
80 final Vector3D actualAcceleration = fieldAcceleration.toVector3D();
81
82 final Vector3D expectedAcceleration = j2OnlyPerturbation.acceleration(spacecraftState, new double[0]);
83 final double tolerance = 1e-12;
84 Assertions.assertEquals(expectedAcceleration.getX(), actualAcceleration.getX(), tolerance);
85 Assertions.assertEquals(expectedAcceleration.getY(), actualAcceleration.getY(), tolerance);
86 Assertions.assertEquals(expectedAcceleration.getZ(), actualAcceleration.getZ(), tolerance);
87 }
88
89 @Test
90 void testAccelerationFieldAgainstGeopotential() {
91
92 final UnnormalizedSphericalHarmonicsProvider unnormalizedProvider = GravityFieldFactory.getUnnormalizedProvider(2 ,0);
93 final Frame frame = FramesFactory.getTOD(false);
94 final J2OnlyPerturbation j2OnlyPerturbation = new J2OnlyPerturbation(unnormalizedProvider, frame);
95 final SpacecraftState spacecraftState = buildState();
96 final ComplexField field = ComplexField.getInstance();
97 final FieldSpacecraftState<Complex> fieldState = new FieldSpacecraftState<>(field, spacecraftState);
98
99 final FieldVector3D<Complex> actualAcceleration = j2OnlyPerturbation.acceleration(fieldState, null);
100
101 final Frame rotatingFrame = FramesFactory.getGTOD(false);
102 final HolmesFeatherstoneAttractionModel holmesFeatherstoneAttractionModel = new HolmesFeatherstoneAttractionModel(
103 rotatingFrame, GravityFieldFactory.getNormalizedProvider(unnormalizedProvider));
104 final Complex[] mu = MathArrays.buildArray(field, 1);
105 mu[0] = field.getOne().newInstance(unnormalizedProvider.getMu());
106 final FieldVector3D<Complex> expectedAcceleration = holmesFeatherstoneAttractionModel.acceleration(fieldState,
107 mu);
108 final double tolerance = 1e-12;
109 Assertions.assertEquals(expectedAcceleration.getX().getReal(), actualAcceleration.getX().getReal(), tolerance);
110 Assertions.assertEquals(expectedAcceleration.getY().getReal(), actualAcceleration.getY().getReal(), tolerance);
111 Assertions.assertEquals(expectedAcceleration.getZ().getReal(), actualAcceleration.getZ().getReal(), tolerance);
112 Assertions.assertEquals(expectedAcceleration.getX().getImaginary(), actualAcceleration.getX().getImaginary(), tolerance);
113 Assertions.assertEquals(expectedAcceleration.getY().getImaginary(), actualAcceleration.getY().getImaginary(), tolerance);
114 Assertions.assertEquals(expectedAcceleration.getZ().getImaginary(), actualAcceleration.getZ().getImaginary(), tolerance);
115 }
116
117 @Test
118 void testGetters() {
119
120 final Frame expectedFrame = FramesFactory.getGTOD(false);
121 final double expectedMu = Constants.EGM96_EARTH_MU;
122 final double expectedrEq = Constants.EGM96_EARTH_EQUATORIAL_RADIUS;
123 final double expectedJ2 = -Constants.EGM96_EARTH_C20;
124
125 final J2OnlyPerturbation j2OnlyPerturbation = new J2OnlyPerturbation(expectedMu, expectedrEq,
126 getJ2OverTime(expectedJ2), expectedFrame);
127
128 Assertions.assertEquals(expectedFrame, j2OnlyPerturbation.getFrame());
129 Assertions.assertEquals(expectedMu, j2OnlyPerturbation.getMu());
130 Assertions.assertEquals(expectedrEq, j2OnlyPerturbation.getrEq());
131 Assertions.assertTrue(j2OnlyPerturbation.getParametersDrivers().isEmpty());
132 Assertions.assertTrue(j2OnlyPerturbation.dependsOnPositionOnly());
133 }
134
135 @Test
136 void testGetJ2() {
137
138 final Frame expectedFrame = FramesFactory.getGTOD(false);
139 final double expectedMu = Constants.EGM96_EARTH_MU;
140 final double expectedrEq = Constants.EGM96_EARTH_EQUATORIAL_RADIUS;
141 final double expectedJ2 = -Constants.EGM96_EARTH_C20;
142 final AbsoluteDate absoluteDate = AbsoluteDate.ARBITRARY_EPOCH;
143
144 final J2OnlyPerturbation j2OnlyPerturbation = new J2OnlyPerturbation(expectedMu, expectedrEq,
145 getJ2OverTime(expectedJ2), expectedFrame);
146
147 final FieldAbsoluteDate<Complex> fieldAbsoluteDate = new FieldAbsoluteDate<>(ComplexField.getInstance(),
148 absoluteDate);
149 Assertions.assertEquals(j2OnlyPerturbation.getJ2(absoluteDate),
150 j2OnlyPerturbation.getJ2(fieldAbsoluteDate).getReal());
151 }
152
153 private TimeScalarFunction getJ2OverTime(final double constantJ2) {
154 return new TimeScalarFunction() {
155 @Override
156 public double value(AbsoluteDate date) {
157 return constantJ2;
158 }
159
160 @Override
161 public <T extends CalculusFieldElement<T>> T value(FieldAbsoluteDate<T> date) {
162 return date.getField().getZero().newInstance(constantJ2);
163 }
164 };
165 }
166
167 private SpacecraftState buildState() {
168 final Vector3D position = new Vector3D(1e7, 1e3, 1e2);
169 final Vector3D velocity = new Vector3D(0., 4e3, 1e1);
170 final PVCoordinates pvCoordinates = new PVCoordinates(position, velocity);
171 final CartesianOrbit orbit = new CartesianOrbit(pvCoordinates, FramesFactory.getGCRF(),
172 AbsoluteDate.ARBITRARY_EPOCH, Constants.EGM96_EARTH_MU);
173 return new SpacecraftState(orbit);
174 }
175
176 }