1   /* Copyright 2002-2022 CS GROUP
2    * Licensed to CS GROUP (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *   http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.orekit.forces.gravity;
18  
19  import org.hipparchus.Field;
20  import org.hipparchus.analysis.differentiation.DerivativeStructure;
21  import org.hipparchus.analysis.differentiation.Gradient;
22  import org.hipparchus.geometry.euclidean.threed.FieldRotation;
23  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
24  import org.hipparchus.geometry.euclidean.threed.Vector3D;
25  import org.hipparchus.ode.AbstractIntegrator;
26  import org.hipparchus.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
27  import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
28  import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator;
29  import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
30  import org.hipparchus.util.FastMath;
31  import org.junit.Assert;
32  import org.junit.Before;
33  import org.junit.Test;
34  import org.orekit.Utils;
35  import org.orekit.attitudes.LofOffset;
36  import org.orekit.forces.AbstractLegacyForceModelTest;
37  import org.orekit.forces.ForceModel;
38  import org.orekit.frames.Frame;
39  import org.orekit.frames.FramesFactory;
40  import org.orekit.frames.LOFType;
41  import org.orekit.orbits.CartesianOrbit;
42  import org.orekit.orbits.FieldKeplerianOrbit;
43  import org.orekit.orbits.KeplerianOrbit;
44  import org.orekit.orbits.Orbit;
45  import org.orekit.orbits.OrbitType;
46  import org.orekit.orbits.PositionAngle;
47  import org.orekit.propagation.FieldSpacecraftState;
48  import org.orekit.propagation.SpacecraftState;
49  import org.orekit.propagation.numerical.FieldNumericalPropagator;
50  import org.orekit.propagation.numerical.NumericalPropagator;
51  import org.orekit.time.AbsoluteDate;
52  import org.orekit.time.DateComponents;
53  import org.orekit.time.FieldAbsoluteDate;
54  import org.orekit.time.TimeComponents;
55  import org.orekit.time.TimeScalesFactory;
56  import org.orekit.utils.Constants;
57  import org.orekit.utils.FieldPVCoordinates;
58  import org.orekit.utils.PVCoordinates;
59  
60  public class DeSitterRelativityTest extends AbstractLegacyForceModelTest {
61  
62      /** speed of light */
63      private static final double c = Constants.SPEED_OF_LIGHT;
64      /** arbitrary date */
65      private static final AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
66      /** inertial frame */
67      private static final Frame frame = FramesFactory.getGCRF();
68  
69      @Override
70      protected FieldVector3D<DerivativeStructure>
71          accelerationDerivatives(ForceModel forceModel, AbsoluteDate date,
72                                  Frame frame,
73                                  FieldVector3D<DerivativeStructure> position,
74                                  FieldVector3D<DerivativeStructure> velocity,
75                                  FieldRotation<DerivativeStructure> rotation,
76                                  DerivativeStructure mass) {
77  
78          final DeSitterRelativity model = (DeSitterRelativity) forceModel;
79  
80          // Useful constant
81          final double c2 = c * c;
82  
83          // Sun's gravitational parameter
84          final double gm = model.getParametersDrivers().get(0).getValue();
85  
86          // Coordinates of the Earth with respect to the Sun
87          final FieldPVCoordinates<DerivativeStructure> pvEarth = model.getEarth().getPVCoordinates(new FieldAbsoluteDate<>(mass.getField(), date), model.getSun().getInertiallyOrientedFrame());
88          final FieldVector3D<DerivativeStructure> pEarth = pvEarth.getPosition();
89          final FieldVector3D<DerivativeStructure> vEarth = pvEarth .getVelocity();
90  
91          // Radius
92          final DerivativeStructure r  = pEarth.getNorm();
93          final DerivativeStructure r3 = r.multiply(r).multiply(r);
94  
95          // Eq. 10.12
96          return new FieldVector3D<>(r3.multiply(c2).reciprocal().multiply(-3.0 * gm), vEarth.crossProduct(pEarth).crossProduct(velocity));
97      }
98  
99      @Override
100     protected FieldVector3D<Gradient>
101         accelerationDerivativesGradient(ForceModel forceModel,
102                                         AbsoluteDate date, Frame frame,
103                                         FieldVector3D<Gradient> position,
104                                         FieldVector3D<Gradient> velocity,
105                                         FieldRotation<Gradient> rotation,
106                                         Gradient mass) {
107 
108         final DeSitterRelativity model = (DeSitterRelativity) forceModel;
109 
110         // Useful constant
111         final double c2 = Constants.SPEED_OF_LIGHT * Constants.SPEED_OF_LIGHT;
112 
113         // Sun's gravitational parameter
114         final double gm = model.getParametersDrivers().get(0).getValue();
115 
116         // Coordinates of the Earth with respect to the Sun
117         final FieldPVCoordinates<Gradient> pvEarth = model.getEarth().getPVCoordinates(new FieldAbsoluteDate<>(mass.getField(), date), model.getSun().getInertiallyOrientedFrame());
118         final FieldVector3D<Gradient> pEarth = pvEarth.getPosition();
119         final FieldVector3D<Gradient> vEarth = pvEarth .getVelocity();
120 
121         // Radius
122         final Gradient r  = pEarth.getNorm();
123         final Gradient r3 = r.multiply(r).multiply(r);
124 
125         // Eq. 10.12
126         return new FieldVector3D<>(r3.multiply(c2).reciprocal().multiply(-3.0 * gm), vEarth.crossProduct(pEarth).crossProduct(velocity));
127     }
128 
129     /**
130      * Check against prediction in
131      * 
132      * "C IUFOLINI, Ignazio, MATZNER, Richard, GURZADYAN, Vahe, et al.
133      * A new laser-ranged satellite for General Relativityand space geodesy:
134      * III. De Sitter effect and the LARES 2 space experiment.
135      * The European Physical Journal C, 2017, vol. 77, no 12, p. 819"
136      * 
137      * They predict a precession of the orbital plane at a rate of the order
138      * of -19.2 milliarcsecs per year.
139      *
140      * As this effect is very small, a propagation of several years is needed to
141      * find the same value. However, it takes a lot of time to be performed.
142      * As a result, we propose a propagation of 60 days for which the impact is
143      * equal to -20.2 milliarcsecs per year.
144      */
145     @Test
146     public void testSmallEffectOnOrbit() {
147         // Setup
148         final double gm = Constants.EIGEN5C_EARTH_MU;
149         Orbit orbit =
150                 new KeplerianOrbit(7000000.0, 0.01, FastMath.toRadians(80.), FastMath.toRadians(80.), FastMath.toRadians(20.),
151                                    FastMath.toRadians(40.), PositionAngle.MEAN,
152                                    frame, date, gm
153                 );
154         double[][] tol = NumericalPropagator.tolerances(0.1, orbit, OrbitType.KEPLERIAN);
155         AbstractIntegrator integrator = new DormandPrince853Integrator(1, 3600, tol[0], tol[1]);
156         NumericalPropagator propagator = new NumericalPropagator(integrator);
157         propagator.setOrbitType(OrbitType.CARTESIAN);
158         propagator.addForceModel(new DeSitterRelativity());
159         propagator.setInitialState(new SpacecraftState(orbit));
160 
161         // Action: propagate a period
162         AbsoluteDate end = orbit.getDate().shiftedBy(60 * Constants.JULIAN_DAY);
163         PVCoordinates actual = propagator.getPVCoordinates(end, frame);
164 
165         // Verify
166         KeplerianOrbit endOrbit = new KeplerianOrbit(actual, frame, end, gm);
167         KeplerianOrbit startOrbit = new KeplerianOrbit(orbit);
168         double dp = startOrbit.getRightAscensionOfAscendingNode() - endOrbit.getRightAscensionOfAscendingNode();
169         double dtYears = end.durationFrom(orbit.getDate()) / Constants.JULIAN_YEAR;
170         double dpDeg = FastMath.toDegrees(dp);
171         // change in right ascension of the ascending node in milliarcseconds per year
172         double milliArcsecPerYear = 1.0e3 * dpDeg * 3600 / dtYears;
173         Assert.assertEquals(-19.2, milliArcsecPerYear, 1.0);
174     }
175 
176     @Test
177     public void testGlobalStateJacobian()
178         {
179 
180         // initialization
181         AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
182                                              new TimeComponents(13, 59, 27.816),
183                                              TimeScalesFactory.getUTC());
184         double i     = FastMath.toRadians(98.7);
185         double omega = FastMath.toRadians(93.0);
186         double OMEGA = FastMath.toRadians(15.0 * 22.5);
187         Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
188                                          0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
189                                          Constants.EIGEN5C_EARTH_MU);
190         OrbitType integrationType = OrbitType.CARTESIAN;
191         double[][] tolerances = NumericalPropagator.tolerances(0.01, orbit, integrationType);
192 
193         NumericalPropagator propagator =
194                 new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
195                                                                        tolerances[0], tolerances[1]));
196         propagator.setOrbitType(integrationType);
197         DeSitterRelativity relativity = new DeSitterRelativity();
198         propagator.addForceModel(relativity);
199         SpacecraftState state0 = new SpacecraftState(orbit);
200 
201         checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0),
202                            1e4, tolerances[0], 2.5e-8);
203 
204     }
205 
206     @Test
207     public void testParameterDerivativeGradient() {
208 
209         final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
210         final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
211         final SpacecraftState state =
212                 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
213                                                        FramesFactory.getGCRF(),
214                                                        new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
215                                                        Constants.EIGEN5C_EARTH_MU));
216 
217         DeSitterRelativity relativity = new DeSitterRelativity();
218         Assert.assertFalse(relativity.dependsOnPositionOnly());
219         final String name = relativity.getSun().getName() + ThirdBodyAttraction.ATTRACTION_COEFFICIENT_SUFFIX;
220         checkParameterDerivativeGradient(state, relativity, name, 1.0, 1.0e-15);
221 
222     }
223 
224     @Test
225     public void RealFieldGradientTest() {
226 
227         final int freeParameters = 6;
228         Gradient a_0 = Gradient.variable(freeParameters, 0, 7e7);
229         Gradient e_0 = Gradient.variable(freeParameters, 1, 0.4);
230         Gradient i_0 = Gradient.variable(freeParameters, 2, 85 * FastMath.PI / 180);
231         Gradient R_0 = Gradient.variable(freeParameters, 3, 0.7);
232         Gradient O_0 = Gradient.variable(freeParameters, 4, 0.5);
233         Gradient n_0 = Gradient.variable(freeParameters, 5, 0.1);
234         Gradient mu  = Gradient.constant(freeParameters, Constants.EIGEN5C_EARTH_MU);
235 
236         Field<Gradient> field = a_0.getField();
237 
238         FieldAbsoluteDate<Gradient> J2000 = new FieldAbsoluteDate<>(field);
239 
240         Frame EME = FramesFactory.getEME2000();
241 
242         FieldKeplerianOrbit<Gradient> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
243                                                                       PositionAngle.MEAN,
244                                                                       EME,
245                                                                       J2000,
246                                                                       mu);
247 
248         FieldSpacecraftState<Gradient> initialState = new FieldSpacecraftState<>(FKO);
249 
250         SpacecraftState iSR = initialState.toSpacecraftState();
251         OrbitType type = OrbitType.KEPLERIAN;
252         double[][] tolerance = NumericalPropagator.tolerances(0.001, FKO.toOrbit(), type);
253 
254 
255         AdaptiveStepsizeFieldIntegrator<Gradient> integrator =
256                         new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
257         integrator.setInitialStepSize(60);
258         AdaptiveStepsizeIntegrator RIntegrator =
259                         new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
260         RIntegrator.setInitialStepSize(60);
261 
262         FieldNumericalPropagator<Gradient> FNP = new FieldNumericalPropagator<>(field, integrator);
263         FNP.setOrbitType(type);
264         FNP.setInitialState(initialState);
265 
266         NumericalPropagator NP = new NumericalPropagator(RIntegrator);
267         NP.setOrbitType(type);
268         NP.setInitialState(iSR);
269 
270         DeSitterRelativity relativity = new DeSitterRelativity();
271 
272         FNP.addForceModel(relativity);
273         NP.addForceModel(relativity);
274         
275         // Do the test
276         checkRealFieldPropagationGradient(FKO, PositionAngle.MEAN, 1005., NP, FNP,
277                                   1.0e-15, 1.3e-2, 2.9e-4, 1.4e-3,
278                                   1, false);
279     }
280 
281     @Test
282     public void testJacobianVs80ImplementationGradient() {
283         double gm = Constants.EIGEN5C_EARTH_MU;
284         DeSitterRelativity relativity = new DeSitterRelativity();
285         final Vector3D p = new Vector3D(3777828.75000531, -5543949.549783845, 2563117.448578311);
286         final Vector3D v = new Vector3D(489.0060271721, -2849.9328929417, -6866.4671013153);
287         SpacecraftState s = new SpacecraftState(new CartesianOrbit(
288                 new PVCoordinates(p, v),
289                 frame,
290                 date,
291                 gm
292         ));
293 
294         checkStateJacobianVs80ImplementationGradient(s, relativity,
295                                              new LofOffset(s.getFrame(), LOFType.LVLH_CCSDS),
296                                              1.0e-15, false);
297     }
298 
299     @Test
300     public void testJacobianVs80Implementation() {
301         double gm = Constants.EIGEN5C_EARTH_MU;
302         DeSitterRelativity relativity = new DeSitterRelativity();
303         final Vector3D p = new Vector3D(3777828.75000531, -5543949.549783845, 2563117.448578311);
304         final Vector3D v = new Vector3D(489.0060271721, -2849.9328929417, -6866.4671013153);
305         SpacecraftState s = new SpacecraftState(new CartesianOrbit(
306                 new PVCoordinates(p, v),
307                 frame,
308                 date,
309                 gm
310         ));
311 
312         checkStateJacobianVs80Implementation(s, relativity,
313                                              new LofOffset(s.getFrame(), LOFType.LVLH_CCSDS),
314                                              1.0e-50, false);
315     }
316 
317     @Before
318     public void setUp() {
319         Utils.setDataRoot("regular-data");
320     }
321 
322 }