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