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.nonstiff.AdaptiveStepsizeFieldIntegrator;
25  import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
26  import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator;
27  import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
28  import org.hipparchus.util.FastMath;
29  import org.junit.jupiter.api.Assertions;
30  import org.junit.jupiter.api.BeforeEach;
31  import org.junit.jupiter.api.Test;
32  import org.orekit.Utils;
33  import org.orekit.attitudes.LofOffset;
34  import org.orekit.forces.AbstractLegacyForceModelTest;
35  import org.orekit.forces.ForceModel;
36  import org.orekit.frames.Frame;
37  import org.orekit.frames.FramesFactory;
38  import org.orekit.frames.LOFType;
39  import org.orekit.frames.StaticTransform;
40  import org.orekit.frames.Transform;
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.PositionAngleType;
47  import org.orekit.propagation.FieldSpacecraftState;
48  import org.orekit.propagation.SpacecraftState;
49  import org.orekit.propagation.ToleranceProvider;
50  import org.orekit.propagation.numerical.FieldNumericalPropagator;
51  import org.orekit.propagation.numerical.NumericalPropagator;
52  import org.orekit.time.AbsoluteDate;
53  import org.orekit.time.DateComponents;
54  import org.orekit.time.FieldAbsoluteDate;
55  import org.orekit.time.TimeComponents;
56  import org.orekit.time.TimeScalesFactory;
57  import org.orekit.utils.Constants;
58  import org.orekit.utils.IERSConventions;
59  import org.orekit.utils.PVCoordinates;
60  
61  public class LenseThirringRelativityTest extends AbstractLegacyForceModelTest {
62  
63      /** speed of light */
64      private static final double c = Constants.SPEED_OF_LIGHT;
65      /** arbitrary date */
66      private static final AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
67      /** inertial frame */
68      private static final Frame frame = FramesFactory.getGCRF();
69  
70      @Override
71      protected FieldVector3D<DerivativeStructure>
72          accelerationDerivatives(ForceModel forceModel, FieldSpacecraftState<DerivativeStructure> state) {
73          try {
74              final AbsoluteDate                       date     = state.getDate().toAbsoluteDate();
75              final FieldVector3D<DerivativeStructure> position = state.getPVCoordinates().getPosition();
76              final FieldVector3D<DerivativeStructure> velocity = state.getPVCoordinates().getVelocity();
77              java.lang.reflect.Field bodyFrameField = LenseThirringRelativity.class.getDeclaredField("bodyFrame");
78              bodyFrameField.setAccessible(true);
79              Frame bodyFrame = (Frame) bodyFrameField.get(forceModel);
80  
81              double gm = forceModel.
82                          getParameterDriver(NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT).
83                          getValue(date);
84              // Radius
85              final DerivativeStructure r  = position.getNorm();
86              final DerivativeStructure r2 = r.multiply(r);
87  
88              // Earth’s angular momentum per unit mass
89              final StaticTransform t = bodyFrame.getStaticTransformTo(frame, date);
90              final Vector3D        j = t.transformVector(Vector3D.PLUS_K).scalarMultiply(9.8e8);
91  
92              return new FieldVector3D<>(position.dotProduct(j).multiply(3.0).divide(r2),
93                                         position.crossProduct(velocity),
94                                         r.getField().getOne(),
95                                         velocity.crossProduct(j))
96                                         .scalarMultiply(r2.multiply(r).multiply(c * c).reciprocal().multiply(2.0 * gm));
97          } catch (IllegalArgumentException | IllegalAccessException | NoSuchFieldException | SecurityException e) {
98              return null;
99          }
100     }
101 
102     @Override
103     protected FieldVector3D<Gradient>
104         accelerationDerivativesGradient(ForceModel forceModel,
105                                         FieldSpacecraftState<Gradient> state) {
106         try {
107             final AbsoluteDate                       date     = state.getDate().toAbsoluteDate();
108             final FieldVector3D<Gradient> position = state.getPVCoordinates().getPosition();
109             final FieldVector3D<Gradient> velocity = state.getPVCoordinates().getVelocity();
110             java.lang.reflect.Field bodyFrameField = LenseThirringRelativity.class.getDeclaredField("bodyFrame");
111             bodyFrameField.setAccessible(true);
112             Frame bodyFrame = (Frame) bodyFrameField.get(forceModel);
113 
114             // Useful constant
115             final double c2 = Constants.SPEED_OF_LIGHT * Constants.SPEED_OF_LIGHT;
116             double gm = forceModel.
117                         getParameterDriver(NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT).
118                         getValue(date);
119             // Radius
120             final Gradient r  = position.getNorm();
121             final Gradient r2 = r.multiply(r);
122 
123             // Earth’s angular momentum per unit mass
124             final Transform t = bodyFrame.getTransformTo(frame, date);
125             final Vector3D  j = t.transformVector(Vector3D.PLUS_K).scalarMultiply(9.8e8);
126 
127             return new FieldVector3D<>(position.dotProduct(j).multiply(3.0).divide(r2),
128                                        position.crossProduct(velocity),
129                                        r.getField().getOne(),
130                                        velocity.crossProduct(j))
131                                        .scalarMultiply(r2.multiply(r).multiply(c2).reciprocal().multiply(2.0 * gm));
132         } catch (IllegalArgumentException | IllegalAccessException | NoSuchFieldException | SecurityException e) {
133             return null;
134         }
135     }
136 
137     @Test
138     void testJacobianVs80Implementation() {
139         double gm = Constants.EIGEN5C_EARTH_MU;
140         LenseThirringRelativity relativity = new LenseThirringRelativity(gm, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
141         final Vector3D p = new Vector3D(3777828.75000531, -5543949.549783845, 2563117.448578311);
142         final Vector3D v = new Vector3D(489.0060271721, -2849.9328929417, -6866.4671013153);
143         SpacecraftState s = new SpacecraftState(new CartesianOrbit(
144                 new PVCoordinates(p, v),
145                 frame,
146                 date,
147                 gm
148         ));
149 
150         checkStateJacobianVs80Implementation(s, relativity,
151                                              new LofOffset(s.getFrame(), LOFType.LVLH_CCSDS),
152                                              1.0e-15, false);
153     }
154 
155     @Test
156     void testJacobianVs80ImplementationGradient() {
157         double gm = Constants.EIGEN5C_EARTH_MU;
158         LenseThirringRelativity relativity = new LenseThirringRelativity(gm, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
159         final Vector3D p = new Vector3D(3777828.75000531, -5543949.549783845, 2563117.448578311);
160         final Vector3D v = new Vector3D(489.0060271721, -2849.9328929417, -6866.4671013153);
161         SpacecraftState s = new SpacecraftState(new CartesianOrbit(
162                 new PVCoordinates(p, v),
163                 frame,
164                 date,
165                 gm
166         ));
167 
168         checkStateJacobianVs80ImplementationGradient(s, relativity,
169                                              new LofOffset(s.getFrame(), LOFType.LVLH_CCSDS),
170                                              1.0e-15, false);
171     }
172 
173     /**Testing if the propagation between the FieldPropagation and the propagation
174      * is equivalent.
175      * Also testing if propagating X+dX with the propagation is equivalent to
176      * propagation X with the FieldPropagation and then applying the taylor
177      * expansion of dX to the result.*/
178     @Test
179     public void RealFieldGradientTest() {
180 
181         final int freeParameters = 6;
182         Gradient a_0 = Gradient.variable(freeParameters, 0, 7e7);
183         Gradient e_0 = Gradient.variable(freeParameters, 1, 0.4);
184         Gradient i_0 = Gradient.variable(freeParameters, 2, 85 * FastMath.PI / 180);
185         Gradient R_0 = Gradient.variable(freeParameters, 3, 0.7);
186         Gradient O_0 = Gradient.variable(freeParameters, 4, 0.5);
187         Gradient n_0 = Gradient.variable(freeParameters, 5, 0.1);
188         Gradient mu  = Gradient.constant(freeParameters, Constants.EIGEN5C_EARTH_MU);
189 
190         Field<Gradient> field = a_0.getField();
191 
192         FieldAbsoluteDate<Gradient> J2000 = new FieldAbsoluteDate<>(field);
193 
194         Frame EME = FramesFactory.getEME2000();
195 
196         FieldKeplerianOrbit<Gradient> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
197                                                                       PositionAngleType.MEAN,
198                                                                       EME,
199                                                                       J2000,
200                                                                       mu);
201 
202         FieldSpacecraftState<Gradient> initialState = new FieldSpacecraftState<>(FKO);
203 
204         SpacecraftState iSR = initialState.toSpacecraftState();
205         OrbitType type = OrbitType.KEPLERIAN;
206         double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(0.001).getTolerances(FKO.toOrbit(), type);
207 
208 
209         AdaptiveStepsizeFieldIntegrator<Gradient> integrator =
210                         new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
211         integrator.setInitialStepSize(60);
212         AdaptiveStepsizeIntegrator RIntegrator =
213                         new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
214         RIntegrator.setInitialStepSize(60);
215 
216         FieldNumericalPropagator<Gradient> FNP = new FieldNumericalPropagator<>(field, integrator);
217         FNP.setOrbitType(type);
218         FNP.setInitialState(initialState);
219 
220         NumericalPropagator NP = new NumericalPropagator(RIntegrator);
221         NP.setOrbitType(type);
222         NP.setInitialState(iSR);
223 
224         LenseThirringRelativity relativity = new LenseThirringRelativity(Constants.EIGEN5C_EARTH_MU, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
225 
226         FNP.addForceModel(relativity);
227         NP.addForceModel(relativity);
228 
229         // Do the test
230         checkRealFieldPropagationGradient(FKO, PositionAngleType.MEAN, 1005., NP, FNP,
231                                   1.0e-15, 1.3e-2, 2.9e-4, 1.4e-3,
232                                   1, false);
233     }
234 
235     @Test
236     void testParameterDerivativeGradient() {
237 
238         final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
239         final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
240         final SpacecraftState state =
241                 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
242                                                        FramesFactory.getGCRF(),
243                                                        new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
244                                                        Constants.EIGEN5C_EARTH_MU));
245 
246         LenseThirringRelativity relativity = new LenseThirringRelativity(Constants.EIGEN5C_EARTH_MU, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
247         Assertions.assertFalse(relativity.dependsOnPositionOnly());
248         final String name = NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT;
249         checkParameterDerivativeGradient(state, relativity, name, 1.0, 1.0e-15);
250 
251     }
252 
253     @Test
254     void testGlobalStateJacobian()
255         {
256 
257         // initialization
258         AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
259                                              new TimeComponents(13, 59, 27.816),
260                                              TimeScalesFactory.getUTC());
261         double i     = FastMath.toRadians(98.7);
262         double omega = FastMath.toRadians(93.0);
263         double OMEGA = FastMath.toRadians(15.0 * 22.5);
264         Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
265                                          0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date,
266                                          Constants.EIGEN5C_EARTH_MU);
267         OrbitType integrationType = OrbitType.CARTESIAN;
268         double[][] tolerances = ToleranceProvider.getDefaultToleranceProvider(0.01).getTolerances(orbit, integrationType);
269 
270         NumericalPropagator propagator =
271                 new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
272                                                                        tolerances[0], tolerances[1]));
273         propagator.setOrbitType(integrationType);
274         LenseThirringRelativity relativity = new LenseThirringRelativity(Constants.EIGEN5C_EARTH_MU, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
275         propagator.addForceModel(relativity);
276         SpacecraftState state0 = new SpacecraftState(orbit);
277 
278         checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0),
279                            1e4, tolerances[0], 2.5e-9);
280 
281     }
282 
283     @BeforeEach
284     public void setUp() {
285         Utils.setDataRoot("regular-data");
286     }
287 
288 }