1   /* Contributed to the public domain
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 java.util.Collections;
20  import java.util.List;
21  
22  import org.hipparchus.CalculusFieldElement;
23  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
24  import org.hipparchus.geometry.euclidean.threed.Vector3D;
25  import org.hipparchus.util.FastMath;
26  import org.orekit.forces.ForceModel;
27  import org.orekit.propagation.FieldSpacecraftState;
28  import org.orekit.propagation.SpacecraftState;
29  import org.orekit.utils.Constants;
30  import org.orekit.utils.FieldPVCoordinates;
31  import org.orekit.utils.PVCoordinates;
32  import org.orekit.utils.ParameterDriver;
33  
34  /**
35   * Post-Newtonian correction force due to general relativity. The main effect is the
36   * precession of perigee by a few arcseconds per year.
37   *
38   * <p> Implemented from Montenbruck and Gill equation 3.146.
39   *
40   * @author Evan Ward
41   * @see "Montenbruck, Oliver, and Gill, Eberhard. Satellite orbits : models, methods, and
42   * applications. Berlin New York: Springer, 2000."
43   */
44  public class Relativity implements ForceModel {
45  
46      /** Central attraction scaling factor.
47       * <p>
48       * We use a power of 2 to avoid numeric noise introduction
49       * in the multiplications/divisions sequences.
50       * </p>
51       */
52      private static final double MU_SCALE = FastMath.scalb(1.0, 32);
53  
54      /** Driver for gravitational parameter. */
55      private final ParameterDriver gmParameterDriver;
56  
57      /**
58       * Create a force model to add post-Newtonian acceleration corrections to an Earth
59       * orbit.
60       *
61       * @param gm Earth's gravitational parameter.
62       */
63      public Relativity(final double gm) {
64          gmParameterDriver = new ParameterDriver(NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT,
65                                                  gm, MU_SCALE,
66                                                  0.0, Double.POSITIVE_INFINITY);
67      }
68  
69      /** {@inheritDoc} */
70      @Override
71      public boolean dependsOnPositionOnly() {
72          return false;
73      }
74  
75      /** {@inheritDoc} */
76      @Override
77      public Vector3D acceleration(final SpacecraftState s, final double[] parameters) {
78  
79          final double gm = parameters[0];
80  
81          final PVCoordinates pv = s.getPVCoordinates();
82          final Vector3D p = pv.getPosition();
83          final Vector3D v = pv.getVelocity();
84          //radius
85          final double r2 = p.getNormSq();
86          final double r = FastMath.sqrt(r2);
87          //speed
88          final double s2 = v.getNormSq();
89          final double c2 = Constants.SPEED_OF_LIGHT * Constants.SPEED_OF_LIGHT;
90          //eq. 3.146
91          return new Vector3D(
92                  4 * gm / r - s2,
93                  p,
94                  4 * p.dotProduct(v),
95                  v)
96                  .scalarMultiply(gm / (r2 * r * c2));
97  
98      }
99  
100     /** {@inheritDoc} */
101     @Override
102     public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(final FieldSpacecraftState<T> s,
103                                                                          final T[] parameters) {
104 
105         final T gm = parameters[0];
106 
107         final FieldPVCoordinates<T> pv = s.getPVCoordinates();
108         final FieldVector3D<T> p = pv.getPosition();
109         final FieldVector3D<T> v = pv.getVelocity();
110         //radius
111         final T r2 = p.getNormSq();
112         final T r = r2.sqrt();
113         //speed
114         final T s2 = v.getNormSq();
115         final double c2 = Constants.SPEED_OF_LIGHT * Constants.SPEED_OF_LIGHT;
116         //eq. 3.146
117         return new FieldVector3D<>(r.reciprocal().multiply(4).multiply(gm).subtract(s2),
118                                    p,
119                                    p.dotProduct(v).multiply(4),
120                                    v).scalarMultiply(r2.multiply(r).multiply(c2).reciprocal().multiply(gm));
121 
122     }
123 
124     /** {@inheritDoc} */
125     @Override
126     public List<ParameterDriver> getParametersDrivers() {
127         return Collections.singletonList(gmParameterDriver);
128     }
129 
130 }