L1TransformProvider.java

  1. /* Copyright 2002-2024 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.frames;

  18. import org.hipparchus.Field;
  19. import org.hipparchus.CalculusFieldElement;
  20. import org.hipparchus.analysis.CalculusFieldUnivariateFunction;
  21. import org.hipparchus.analysis.UnivariateFunction;
  22. import org.hipparchus.analysis.solvers.AllowedSolution;
  23. import org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver;
  24. import org.hipparchus.analysis.solvers.FieldBracketingNthOrderBrentSolver;
  25. import org.hipparchus.analysis.solvers.UnivariateSolverUtils;
  26. import org.hipparchus.geometry.euclidean.threed.FieldRotation;
  27. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  28. import org.hipparchus.geometry.euclidean.threed.Rotation;
  29. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  30. import org.hipparchus.util.FastMath;
  31. import org.orekit.bodies.CelestialBody;
  32. import org.orekit.time.AbsoluteDate;
  33. import org.orekit.time.FieldAbsoluteDate;
  34. import org.orekit.utils.FieldPVCoordinates;
  35. import org.orekit.utils.PVCoordinates;

  36. /** L1 Transform provider for a frame on the L1 Lagrange point of two celestial bodies.
  37.  *
  38.  * @author Luc Maisonobe
  39.  * @author Julio Hernanz
  40.  */
  41. public class L1TransformProvider implements TransformProvider {

  42.     /** Relative accuracy on position for solver. */
  43.     private static final double RELATIVE_ACCURACY = 1e-14;

  44.     /** Absolute accuracy on position for solver (1mm). */
  45.     private static final double ABSOLUTE_ACCURACY = 1e-3;

  46.     /** Function value ccuracy for solver (set to 0 so we rely only on position for convergence). */
  47.     private static final double FUNCTION_ACCURACY = 0;

  48.     /** Maximal order for solver. */
  49.     private static final int MAX_ORDER = 5;

  50.     /** Maximal number of evaluations for solver. */
  51.     private static final int MAX_EVALUATIONS = 1000;

  52.     /** Serializable UID.*/
  53.     private static final long serialVersionUID = 20170824L;

  54.     /** Frame for results. Always defined as primaryBody's inertially oriented frame.*/
  55.     private final Frame frame;

  56.     /** Celestial body with bigger mass, m1.*/
  57.     private final CelestialBody primaryBody;

  58.     /** Celestial body with smaller mass, m2.*/
  59.     private final CelestialBody secondaryBody;

  60.     /** Simple constructor.
  61.      * @param primaryBody Primary body.
  62.      * @param secondaryBody Secondary body.
  63.      */
  64.     public L1TransformProvider(final CelestialBody primaryBody, final CelestialBody secondaryBody) {
  65.         this.primaryBody   = primaryBody;
  66.         this.secondaryBody = secondaryBody;
  67.         this.frame         = primaryBody.getInertiallyOrientedFrame();
  68.     }

  69.     /** {@inheritDoc} */
  70.     @Override
  71.     public Transform getTransform(final AbsoluteDate date) {
  72.         final PVCoordinates pv21        = secondaryBody.getPVCoordinates(date, frame);
  73.         final Vector3D      translation = getL1(pv21.getPosition()).negate();
  74.         final Rotation      rotation    = new Rotation(pv21.getPosition(), pv21.getVelocity(),
  75.                                                        Vector3D.PLUS_I, Vector3D.PLUS_J);
  76.         return new Transform(date, new Transform(date, translation), new Transform(date, rotation));
  77.     }

  78.     /** {@inheritDoc} */
  79.     @Override
  80.     public StaticTransform getStaticTransform(final AbsoluteDate date) {
  81.         final PVCoordinates pv21        = secondaryBody.getPVCoordinates(date, frame);
  82.         final Vector3D      translation = getL1(pv21.getPosition()).negate();
  83.         final Rotation      rotation    = new Rotation(pv21.getPosition(), pv21.getVelocity(),
  84.                 Vector3D.PLUS_I, Vector3D.PLUS_J);
  85.         return StaticTransform.compose(
  86.                 date,
  87.                 StaticTransform.of(date, translation),
  88.                 StaticTransform.of(date, rotation));
  89.     }

  90.     /** {@inheritDoc} */
  91.     @Override
  92.     public <T extends CalculusFieldElement<T>> FieldTransform<T> getTransform(final FieldAbsoluteDate<T> date) {
  93.         final FieldPVCoordinates<T> pv21        = secondaryBody.getPVCoordinates(date, frame);
  94.         final FieldVector3D<T>      translation = getL1(pv21.getPosition()).negate();
  95.         final Field<T>              field       = pv21.getPosition().getX().getField();
  96.         final FieldRotation<T>      rotation    = new FieldRotation<>(pv21.getPosition(), pv21.getVelocity(),
  97.                                                                       FieldVector3D.getPlusI(field),
  98.                                                                       FieldVector3D.getPlusJ(field));
  99.         return new FieldTransform<>(date,
  100.                                     new FieldTransform<>(date, translation),
  101.                                     new FieldTransform<>(date, rotation));
  102.     }

  103.     /** {@inheritDoc} */
  104.     @Override
  105.     public <T extends CalculusFieldElement<T>> FieldStaticTransform<T> getStaticTransform(final FieldAbsoluteDate<T> date) {
  106.         final FieldPVCoordinates<T> pv21        = secondaryBody.getPVCoordinates(date, frame);
  107.         final FieldVector3D<T>      translation = getL1(pv21.getPosition()).negate();
  108.         final FieldRotation<T>      rotation    = new FieldRotation<>(pv21.getPosition(), pv21.getVelocity(),
  109.                 FieldVector3D.getPlusI(date.getField()), FieldVector3D.getPlusJ(date.getField()));
  110.         return FieldStaticTransform.compose(
  111.                 date,
  112.                 FieldStaticTransform.of(date, translation),
  113.                 FieldStaticTransform.of(date, rotation));
  114.     }

  115.     /** Compute the coordinates of the L1 point.
  116.      * @param primaryToSecondary relative position of secondary body with respect to primary body
  117.      * @return coordinates of the L1 point given in frame: primaryBody.getInertiallyOrientedFrame()
  118.      */
  119.     private Vector3D getL1(final Vector3D primaryToSecondary) {

  120.         // mass ratio
  121.         final double massRatio = secondaryBody.getGM() / primaryBody.getGM();

  122.         // Approximate position of L1 point, valid when m2 << m1
  123.         final double bigR  = primaryToSecondary.getNorm();
  124.         final double baseR = bigR * (1 - FastMath.cbrt(massRatio / 3));

  125.         // Accurate position of L1 point, by solving the L1 equilibrium equation
  126.         final UnivariateFunction l1Equation = r -> {
  127.             final double bigrminusR  = bigR - r;
  128.             final double lhs         = 1.0 / ( r * r );
  129.             final double rhs1        = massRatio / (bigrminusR * bigrminusR);
  130.             final double rhs2        = 1.0 / (bigR * bigR);
  131.             final double rhs3        = (1 + massRatio) * bigrminusR * rhs2 / bigR;
  132.             return lhs - (rhs1 + rhs2 - rhs3);
  133.         };
  134.         final double[] searchInterval = UnivariateSolverUtils.bracket(l1Equation,
  135.                                                                       baseR, 0, bigR,
  136.                                                                       0.01 * bigR, 1, MAX_EVALUATIONS);
  137.         final BracketingNthOrderBrentSolver solver =
  138.                         new BracketingNthOrderBrentSolver(RELATIVE_ACCURACY,
  139.                                                           ABSOLUTE_ACCURACY,
  140.                                                           FUNCTION_ACCURACY,
  141.                                                           MAX_ORDER);
  142.         final double r = solver.solve(MAX_EVALUATIONS, l1Equation,
  143.                                       searchInterval[0], searchInterval[1],
  144.                                       AllowedSolution.ANY_SIDE);

  145.         // L1 point is built
  146.         return new Vector3D(r / bigR, primaryToSecondary);

  147.     }

  148.     /** Compute the coordinates of the L1 point.
  149.      * @param <T> type of the field elements
  150.      * @param primaryToSecondary relative position of secondary body with respect to primary body
  151.      * @return coordinates of the L1 point given in frame: primaryBody.getInertiallyOrientedFrame()
  152.      */
  153.     private <T extends CalculusFieldElement<T>> FieldVector3D<T>
  154.         getL1(final FieldVector3D<T> primaryToSecondary) {

  155.         // mass ratio
  156.         final double massRatio = secondaryBody.getGM() / primaryBody.getGM();

  157.         // Approximate position of L1 point, valid when m2 << m1
  158.         final T bigR  = primaryToSecondary.getNorm();
  159.         final T baseR = bigR.multiply(1 - FastMath.cbrt(massRatio / 3));

  160.         // Accurate position of L1 point, by solving the L1 equilibrium equation
  161.         final CalculusFieldUnivariateFunction<T> l1Equation = r -> {
  162.             final T bigrminusR = bigR.subtract(r);
  163.             final T lhs        = r.multiply(r).reciprocal();
  164.             final T rhs1       = bigrminusR.multiply(bigrminusR).reciprocal().multiply(massRatio);
  165.             final T rhs2       = bigR.multiply(bigR).reciprocal();
  166.             final T rhs3       = bigrminusR.multiply(rhs2).multiply(1 + massRatio).divide(bigR);
  167.             return lhs.subtract(rhs1.add(rhs2).add(rhs3));
  168.         };
  169.         final T zero             = primaryToSecondary.getX().getField().getZero();
  170.         final T[] searchInterval = UnivariateSolverUtils.bracket(l1Equation,
  171.                                                                  baseR, zero, bigR.multiply(2),
  172.                                                                  bigR.multiply(0.01), zero.add(1),
  173.                                                                  MAX_EVALUATIONS);


  174.         final T relativeAccuracy = zero.add(RELATIVE_ACCURACY);
  175.         final T absoluteAccuracy = zero.add(ABSOLUTE_ACCURACY);
  176.         final T functionAccuracy = zero.add(FUNCTION_ACCURACY);

  177.         final FieldBracketingNthOrderBrentSolver<T> solver =
  178.                         new FieldBracketingNthOrderBrentSolver<>(relativeAccuracy,
  179.                                                                  absoluteAccuracy,
  180.                                                                  functionAccuracy,
  181.                                                                  MAX_ORDER);
  182.         final T r = solver.solve(MAX_EVALUATIONS, l1Equation,
  183.                                  searchInterval[0], searchInterval[1],
  184.                                  AllowedSolution.ANY_SIDE);

  185.         // L1 point is built
  186.         return new FieldVector3D<>(r.divide(bigR), primaryToSecondary);

  187.     }

  188. }