CR3BPForceModel.java

  1. /* Copyright 2002-2020 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.propagation.numerical.cr3bp;

  18. import java.util.stream.Stream;

  19. import org.hipparchus.Field;
  20. import org.hipparchus.RealFieldElement;
  21. import org.hipparchus.analysis.differentiation.DSFactory;
  22. import org.hipparchus.analysis.differentiation.DerivativeStructure;
  23. import org.hipparchus.analysis.differentiation.FDSFactory;
  24. import org.hipparchus.analysis.differentiation.FieldDerivativeStructure;
  25. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  26. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  27. import org.hipparchus.util.FastMath;
  28. import org.orekit.bodies.CR3BPSystem;
  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.ParameterDriver;

  35. /** Class calculating the acceleration induced by CR3BP model.
  36.  * @see "Dynamical systems, the three-body problem, and space mission design, Koon, Lo, Marsden, Ross"
  37.  * @author Vincent Mouraux
  38.  * @since 10.2
  39.  */
  40. public class CR3BPForceModel extends AbstractForceModel {

  41.     /** Suffix for parameter name for Mass Ratio enabling Jacobian processing. */
  42.     public static final String MASS_RATIO_SUFFIX =  "CR3BP System Mass Ratio";

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

  51.     /** Driver for gravitational parameter. */
  52.     private final ParameterDriver muParameterDriver;

  53.     /** Simple constructor.
  54.      * @param cr3bp Name of the CR3BP System
  55.      */
  56.     public CR3BPForceModel(final CR3BPSystem cr3bp) {
  57.         muParameterDriver = new ParameterDriver(cr3bp.getName() + MASS_RATIO_SUFFIX,
  58.                                                 cr3bp.getMassRatio(), MU_SCALE, 0.0,
  59.                                                 Double.POSITIVE_INFINITY);
  60.     }

  61.     /** {@inheritDoc} */
  62.     public Vector3D acceleration(final SpacecraftState s,
  63.                                  final double[] parameters) {

  64.         // Spacecraft Velocity
  65.         final double vx = s.getPVCoordinates().getVelocity().getX();
  66.         final double vy = s.getPVCoordinates().getVelocity().getY();

  67.         // Spacecraft Potential
  68.         final DerivativeStructure potential = getPotential(s);

  69.         // Potential derivatives
  70.         final double[] dU = potential.getAllDerivatives();

  71.         // first order derivatives index
  72.         final int idX = potential.getFactory().getCompiler().getPartialDerivativeIndex(1, 0, 0);
  73.         final int idY = potential.getFactory().getCompiler().getPartialDerivativeIndex(0, 1, 0);
  74.         final int idZ = potential.getFactory().getCompiler().getPartialDerivativeIndex(0, 0, 1);

  75.         // Acceleration calculation according to CR3BP Analytical Model
  76.         final double accx = dU[idX] + 2.0 * vy;
  77.         final double accy = dU[idY] - 2.0 * vx;
  78.         final double accz = dU[idZ];

  79.         // compute absolute acceleration
  80.         return new Vector3D(accx, accy, accz);

  81.     }

  82.     /** {@inheritDoc} */
  83.     public <T extends RealFieldElement<T>> FieldVector3D<T> acceleration(final FieldSpacecraftState<T> s,
  84.                                                                          final T[] parameters) {

  85.         // Spacecraft Velocity
  86.         final T vx = s.getPVCoordinates().getVelocity().getX();
  87.         final T vy = s.getPVCoordinates().getVelocity().getY();

  88.         // Spacecraft Potential
  89.         final FieldDerivativeStructure<T> fieldPotential = getPotential(s);
  90.         // Potential derivatives
  91.         final T[] dU = fieldPotential.getAllDerivatives();

  92.         // first order derivatives index
  93.         final int idX = fieldPotential.getFactory().getCompiler().getPartialDerivativeIndex(1, 0, 0);
  94.         final int idY = fieldPotential.getFactory().getCompiler().getPartialDerivativeIndex(0, 1, 0);
  95.         final int idZ = fieldPotential.getFactory().getCompiler().getPartialDerivativeIndex(0, 0, 1);

  96.         // Acceleration calculation according to CR3BP Analytical Model
  97.         final T accx = dU[idX].add(vy.multiply(2.0));
  98.         final T accy = dU[idY].subtract(vx.multiply(2.0));
  99.         final T accz = dU[idZ];

  100.         // compute absolute acceleration
  101.         return new FieldVector3D<>(accx, accy, accz);

  102.     }

  103.     /**
  104.      * Calculate spacecraft potential.
  105.      * @param s SpacecraftState
  106.      * @return Spacecraft Potential
  107.      */
  108.     public DerivativeStructure getPotential(final SpacecraftState s) {

  109.         // Spacecraft Position
  110.         final double x = s.getPVCoordinates().getPosition().getX();
  111.         final double y = s.getPVCoordinates().getPosition().getY();
  112.         final double z = s.getPVCoordinates().getPosition().getZ();

  113.         final DSFactory factoryP = new DSFactory(3, 2);
  114.         final DerivativeStructure fpx = factoryP.variable(0, x);
  115.         final DerivativeStructure fpy = factoryP.variable(1, y);
  116.         final DerivativeStructure fpz = factoryP.variable(2, z);

  117.         final DerivativeStructure zero = fpx.getField().getZero();

  118.         // Get CR3BP System mass ratio
  119.         final DerivativeStructure mu = zero.add(muParameterDriver.getValue());

  120.         // Normalized distances between primaries and barycenter in CR3BP
  121.         final DerivativeStructure d1 = mu;
  122.         final DerivativeStructure d2 = mu.negate().add(1.0);

  123.         // Norm of the Spacecraft position relative to the primary body
  124.         final DerivativeStructure r1 =
  125.             FastMath.sqrt((fpx.add(d1)).multiply(fpx.add(d1)).add(fpy.multiply(fpy))
  126.                 .add(fpz.multiply(fpz)));

  127.         // Norm of the Spacecraft position relative to the secondary body
  128.         final DerivativeStructure r2 =
  129.             FastMath.sqrt((fpx.subtract(d2)).multiply(fpx.subtract(d2))
  130.                 .add(fpy.multiply(fpy)).add(fpz.multiply(fpz)));

  131.         // Potential of the Spacecraft
  132.         return (mu.negate().add(1.0).divide(r1)).add(mu.divide(r2))
  133.                 .add(fpx.multiply(fpx).add(fpy.multiply(fpy)).multiply(0.5)).add(d1.multiply(d2).multiply(0.5));
  134.     }

  135.     /**
  136.      * Calculate spacecraft potential.
  137.      * @param <T> Field element
  138.      * @param s SpacecraftState
  139.      * @return Spacecraft Potential
  140.      */
  141.     public <T extends RealFieldElement<T>> FieldDerivativeStructure<T> getPotential(final FieldSpacecraftState<T> s) {

  142.         // Spacecraft Position
  143.         final T x = s.getPVCoordinates().getPosition().getX();
  144.         final T y = s.getPVCoordinates().getPosition().getY();
  145.         final T z = s.getPVCoordinates().getPosition().getZ();

  146.         final FDSFactory<T> factoryP = new FDSFactory<>(s.getDate().getField(), 3, 2);
  147.         final FieldDerivativeStructure<T> fpx = factoryP.variable(0, x);
  148.         final FieldDerivativeStructure<T> fpy = factoryP.variable(1, y);
  149.         final FieldDerivativeStructure<T> fpz = factoryP.variable(2, z);
  150.         final FieldDerivativeStructure<T> zero = fpx.getField().getZero();

  151.         // Get CR3BP System mass ratio
  152.         final FieldDerivativeStructure<T> mu = zero.add(muParameterDriver.getValue());

  153.         // Normalized distances between primaries and barycenter in CR3BP
  154.         final FieldDerivativeStructure<T> d1 = mu;
  155.         final FieldDerivativeStructure<T> d2 = mu.negate().add(1.0);

  156.         // Norm of the Spacecraft position relative to the primary body
  157.         final FieldDerivativeStructure<T> r1 =
  158.             FastMath.sqrt((fpx.add(d1)).multiply(fpx.add(d1)).add(fpy.multiply(fpy))
  159.                 .add(fpz.multiply(fpz)));

  160.         // Norm of the Spacecraft position relative to the secondary body
  161.         final FieldDerivativeStructure<T> r2 =
  162.             FastMath.sqrt((fpx.subtract(d2)).multiply(fpx.subtract(d2))
  163.                 .add(fpy.multiply(fpy)).add(fpz.multiply(fpz)));

  164.         // Potential of the Spacecraft
  165.         return (mu.negate().add(1.0).divide(r1)).add(mu.divide(r2))
  166.                 .add(fpx.multiply(fpx).add(fpy.multiply(fpy)).multiply(0.5)).add(d1.multiply(d2).multiply(0.5));
  167.     }

  168.     /** {@inheritDoc} */
  169.     public Stream<EventDetector> getEventsDetectors() {
  170.         return Stream.empty();
  171.     }

  172.     @Override
  173.     /** {@inheritDoc} */
  174.     public <T extends RealFieldElement<T>> Stream<FieldEventDetector<T>> getFieldEventsDetectors(final Field<T> field) {
  175.         return Stream.empty();
  176.     }

  177.     /** {@inheritDoc} */
  178.     public ParameterDriver[] getParametersDrivers() {
  179.         return new ParameterDriver[] {
  180.             muParameterDriver
  181.         };
  182.     }

  183.     /** {@inheritDoc} */
  184.     public boolean dependsOnPositionOnly() {
  185.         return true;
  186.     }
  187. }