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.analysis.differentiation.Gradient;
20  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
21  import org.hipparchus.geometry.euclidean.threed.Vector3D;
22  import org.orekit.bodies.CelestialBody;
23  import org.orekit.propagation.SpacecraftState;
24  
25  /** Third body attraction force model.
26   * This class is a copy of {@link ThirdBodyAttraction} class.
27   * The computation of derivatives of the acceleration
28   * w.r.t. the Epoch has been added.
29   *
30   * @author Fabien Maussion
31   * @author Véronique Pommier-Maurussane
32   * @since 10.2
33   */
34  public class ThirdBodyAttractionEpoch extends ThirdBodyAttraction {
35  
36      /** Simple constructor.
37       * @param body the third body to consider
38       * (ex: {@link org.orekit.bodies.CelestialBodyFactory#getSun()} or
39       * {@link org.orekit.bodies.CelestialBodyFactory#getMoon()})
40       */
41      public ThirdBodyAttractionEpoch(final CelestialBody body) {
42          super(body);
43      }
44  
45      /** Compute acceleration.
46       * @param s current state information: date, kinematics, attitude
47       * @param parameters values of the force model parameters
48       * @return acceleration in same frame as state
49       */
50      private FieldVector3D<Gradient> accelerationToEpoch(final SpacecraftState s, final double[] parameters) {
51  
52          final double gm = parameters[0];
53  
54          // compute bodies separation vectors and squared norm
55          final Vector3D centralToBody = getBodyPosition(s.getDate(), s.getFrame());
56  
57          // Spacecraft Position
58          final double rx = centralToBody.getX();
59          final double ry = centralToBody.getY();
60          final double rz = centralToBody.getZ();
61  
62          final int freeParameters = 3;
63          final Gradient fpx = Gradient.variable(freeParameters, 0, rx);
64          final Gradient fpy = Gradient.variable(freeParameters, 1, ry);
65          final Gradient fpz = Gradient.variable(freeParameters, 2, rz);
66  
67          final FieldVector3D<Gradient> centralToBodyFV = new FieldVector3D<>(new Gradient[] {fpx, fpy, fpz});
68  
69  
70          final Gradient                r2Central = centralToBodyFV.getNormSq();
71          final FieldVector3D<Gradient> satToBody = centralToBodyFV.subtract(s.getPosition());
72          final Gradient                r2Sat     = satToBody.getNormSq();
73  
74          return new FieldVector3D<>(gm, satToBody.scalarMultiply(r2Sat.multiply(r2Sat.sqrt()).reciprocal()),
75                                    -gm, centralToBodyFV.scalarMultiply(r2Central.multiply(r2Central.sqrt()).reciprocal()));
76      }
77  
78      /** Compute derivatives of the state w.r.t epoch.
79       * @param s current state information: date, kinematics, attitude
80       * @param parameters values of the force model parameters
81       * @return derivatives
82       */
83      public double[] getDerivativesToEpoch(final SpacecraftState s, final double[] parameters) {
84  
85          final FieldVector3D<Gradient> acc = accelerationToEpoch(s, parameters);
86          final Vector3D centralToBodyVelocity = getBodyPVCoordinates(s.getDate(), s.getFrame()).getVelocity();
87  
88          final double[] dAccxdR1i = acc.getX().getGradient();
89          final double[] dAccydR1i = acc.getY().getGradient();
90          final double[] dAcczdR1i = acc.getZ().getGradient();
91          final double[] v = centralToBodyVelocity.toArray();
92  
93          return new double[] {
94              dAccxdR1i[0] * v[0] + dAccxdR1i[1] * v[1] + dAccxdR1i[2] * v[2],
95              dAccydR1i[0] * v[0] + dAccydR1i[1] * v[1] + dAccydR1i[2] * v[2],
96              dAcczdR1i[0] * v[0] + dAcczdR1i[1] * v[1] + dAcczdR1i[2] * v[2]
97          };
98  
99      }
100 
101 }