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.stream.Stream;

  19. import org.hipparchus.Field;
  20. import org.hipparchus.RealFieldElement;
  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.AbstractForceModel;
  28. import org.orekit.propagation.FieldSpacecraftState;
  29. import org.orekit.propagation.SpacecraftState;
  30. import org.orekit.propagation.events.EventDetector;
  31. import org.orekit.propagation.events.FieldEventDetector;
  32. import org.orekit.utils.Constants;
  33. import org.orekit.utils.FieldPVCoordinates;
  34. import org.orekit.utils.PVCoordinates;
  35. import org.orekit.utils.ParameterDriver;

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

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

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

  57.     /** The Sun. */
  58.     private final CelestialBody sun;

  59.     /** The Earth. */
  60.     private final CelestialBody earth;

  61.     /** Driver for gravitational parameter. */
  62.     private final ParameterDriver gmParameterDriver;

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

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

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

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

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

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

  106.         // Useful constant
  107.         final double c2 = Constants.SPEED_OF_LIGHT * Constants.SPEED_OF_LIGHT;

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

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

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

  117.         // Radius
  118.         final double r  = pEarth.getNorm();
  119.         final double r3 = r * r * r;

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

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

  127.         // Useful constant
  128.         final double c2 = Constants.SPEED_OF_LIGHT * Constants.SPEED_OF_LIGHT;

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

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

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

  138.         // Radius
  139.         final T r  = pEarth.getNorm();
  140.         final T r3 = r.multiply(r).multiply(r);

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

  144.     /** {@inheritDoc} */
  145.     @Override
  146.     public Stream<EventDetector> getEventsDetectors() {
  147.         return Stream.empty();
  148.     }

  149.     /** {@inheritDoc} */
  150.     @Override
  151.     public <T extends RealFieldElement<T>> Stream<FieldEventDetector<T>> getFieldEventsDetectors(final Field<T> field) {
  152.         return Stream.empty();
  153.     }

  154.     /** {@inheritDoc} */
  155.     @Override
  156.     public ParameterDriver[] getParametersDrivers() {
  157.         return new ParameterDriver[] {
  158.             gmParameterDriver
  159.         };
  160.     }

  161. }