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

  21. import org.hipparchus.Field;
  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.annotation.DefaultDataContext;
  27. import org.orekit.bodies.CelestialBody;
  28. import org.orekit.data.DataContext;
  29. import org.orekit.forces.AbstractForceModel;
  30. import org.orekit.propagation.FieldSpacecraftState;
  31. import org.orekit.propagation.SpacecraftState;
  32. import org.orekit.propagation.events.EventDetector;
  33. import org.orekit.propagation.events.FieldEventDetector;
  34. import org.orekit.utils.Constants;
  35. import org.orekit.utils.FieldPVCoordinates;
  36. import org.orekit.utils.PVCoordinates;
  37. import org.orekit.utils.ParameterDriver;

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

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

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

  59.     /** The Sun. */
  60.     private final CelestialBody sun;

  61.     /** The Earth. */
  62.     private final CelestialBody earth;

  63.     /** Driver for gravitational parameter. */
  64.     private final ParameterDriver gmParameterDriver;

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

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

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

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

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

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

  108.         // Useful constant
  109.         final double c2 = Constants.SPEED_OF_LIGHT * Constants.SPEED_OF_LIGHT;

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

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

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

  119.         // Radius
  120.         final double r  = pEarth.getNorm();
  121.         final double r3 = r * r * r;

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

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

  129.         // Useful constant
  130.         final double c2 = Constants.SPEED_OF_LIGHT * Constants.SPEED_OF_LIGHT;

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

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

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

  140.         // Radius
  141.         final T r  = pEarth.getNorm();
  142.         final T r3 = r.multiply(r).multiply(r);

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

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

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

  156.     /** {@inheritDoc} */
  157.     @Override
  158.     public List<ParameterDriver> getParametersDrivers() {
  159.         return Collections.singletonList(gmParameterDriver);
  160.     }

  161. }