DeSitterRelativity.java

  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. import java.util.Collections;
  19. import java.util.List;

  20. import org.hipparchus.CalculusFieldElement;
  21. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  22. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  23. import org.hipparchus.util.FastMath;
  24. import org.orekit.annotation.DefaultDataContext;
  25. import org.orekit.bodies.CelestialBody;
  26. import org.orekit.data.DataContext;
  27. import org.orekit.forces.ForceModel;
  28. import org.orekit.propagation.FieldSpacecraftState;
  29. import org.orekit.propagation.SpacecraftState;
  30. import org.orekit.utils.Constants;
  31. import org.orekit.utils.FieldPVCoordinates;
  32. import org.orekit.utils.PVCoordinates;
  33. import org.orekit.utils.ParameterDriver;

  34. /**
  35.  * De Sitter post-Newtonian correction force due to general relativity.
  36.  * <p>
  37.  * De Sitter term causes a precession of the orbital plane at a rate of 19 mas per year.
  38.  * </p>
  39.  * @see "Petit, G. and Luzum, B. (eds.), IERS Conventions (2010), Chapter 10,
  40.  * General relativistic models for space-time coordinates and equations of motion (2010)"
  41.  *
  42.  * @author Bryan Cazabonne
  43.  * @since 10.3
  44.  */
  45. public class DeSitterRelativity implements ForceModel {

  46.     /** Suffix for parameter name for attraction coefficient enabling Jacobian processing. */
  47.     public static final String ATTRACTION_COEFFICIENT_SUFFIX = " attraction coefficient";

  48.     /** Central attraction scaling factor.
  49.      * <p>
  50.      * We use a power of 2 to avoid numeric noise introduction
  51.      * in the multiplications/divisions sequences.
  52.      * </p>
  53.      */
  54.     private static final double MU_SCALE = FastMath.scalb(1.0, 32);

  55.     /** The Sun. */
  56.     private final CelestialBody sun;

  57.     /** The Earth. */
  58.     private final CelestialBody earth;

  59.     /** Driver for gravitational parameter. */
  60.     private final ParameterDriver gmParameterDriver;

  61.     /**
  62.      * Constructor.
  63.      * <p>It uses the {@link DataContext#getDefault()} to initialize the celestial bodies.</p>
  64.      */
  65.     @DefaultDataContext
  66.     public DeSitterRelativity() {
  67.         this(DataContext.getDefault().getCelestialBodies().getEarth(),
  68.              DataContext.getDefault().getCelestialBodies().getSun());
  69.     }

  70.     /**
  71.      * Simple constructor.
  72.      * @param earth the Earth
  73.      * @param sun the Sun
  74.      */
  75.     public DeSitterRelativity(final CelestialBody earth, final CelestialBody sun) {
  76.         gmParameterDriver = new ParameterDriver(sun.getName() + ThirdBodyAttraction.ATTRACTION_COEFFICIENT_SUFFIX,
  77.                                                 sun.getGM(), MU_SCALE,
  78.                                                 0.0, Double.POSITIVE_INFINITY);
  79.         this.earth = earth;
  80.         this.sun   = sun;
  81.     }

  82.     /**
  83.      * Get the sun model used to compute De Sitter effect.
  84.      * @return the sun model
  85.      */
  86.     public CelestialBody getSun() {
  87.         return sun;
  88.     }

  89.     /**
  90.      * Get the Earth model used to compute De Sitter effect.
  91.      * @return the earth model
  92.      */
  93.     public CelestialBody getEarth() {
  94.         return earth;
  95.     }

  96.     /** {@inheritDoc} */
  97.     @Override
  98.     public boolean dependsOnPositionOnly() {
  99.         return false;
  100.     }

  101.     /** {@inheritDoc} */
  102.     @Override
  103.     public Vector3D acceleration(final SpacecraftState s, final double[] parameters) {

  104.         // Useful constant
  105.         final double c2 = Constants.SPEED_OF_LIGHT * Constants.SPEED_OF_LIGHT;

  106.         // Sun's gravitational parameter
  107.         final double gm = parameters[0];

  108.         // Satellite velocity with respect to the Earth
  109.         final PVCoordinates pvSat = s.getPVCoordinates();
  110.         final Vector3D vSat = pvSat.getVelocity();

  111.         // Coordinates of the Earth with respect to the Sun
  112.         final PVCoordinates pvEarth = earth.getPVCoordinates(s.getDate(), sun.getInertiallyOrientedFrame());
  113.         final Vector3D pEarth = pvEarth.getPosition();
  114.         final Vector3D vEarth = pvEarth.getVelocity();

  115.         // Radius
  116.         final double r  = pEarth.getNorm();
  117.         final double r3 = r * r * r;

  118.         // Eq. 10.12
  119.         return new Vector3D((-3.0 * gm) / (c2 * r3), vEarth.crossProduct(pEarth).crossProduct(vSat));
  120.     }

  121.     /** {@inheritDoc} */
  122.     @Override
  123.     public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(final FieldSpacecraftState<T> s,
  124.                                                                          final T[] parameters) {

  125.         // Useful constant
  126.         final double c2 = Constants.SPEED_OF_LIGHT * Constants.SPEED_OF_LIGHT;

  127.         // Sun's gravitational parameter
  128.         final T gm = parameters[0];

  129.         // Satellite velocity with respect to the Earth
  130.         final FieldPVCoordinates<T> pvSat = s.getPVCoordinates();
  131.         final FieldVector3D<T> vSat = pvSat.getVelocity();

  132.         // Coordinates of the Earth with respect to the Sun
  133.         final FieldPVCoordinates<T> pvEarth = earth.getPVCoordinates(s.getDate(), sun.getInertiallyOrientedFrame());
  134.         final FieldVector3D<T> pEarth = pvEarth.getPosition();
  135.         final FieldVector3D<T> vEarth = pvEarth .getVelocity();

  136.         // Radius
  137.         final T r  = pEarth.getNorm();
  138.         final T r3 = r.multiply(r).multiply(r);

  139.         // Eq. 10.12
  140.         return new FieldVector3D<>(gm.multiply(-3.0).divide(r3.multiply(c2)), vEarth.crossProduct(pEarth).crossProduct(vSat));
  141.     }

  142.     /** {@inheritDoc} */
  143.     @Override
  144.     public List<ParameterDriver> getParametersDrivers() {
  145.         return Collections.singletonList(gmParameterDriver);
  146.     }

  147. }