L1TransformProvider.java

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

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

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

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

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

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

  54.     /** Maximal order for solver. */
  55.     private static final int MAX_ORDER = 5;

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

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

  63.     /** Celestial body with bigger mass, m1.*/
  64.     private final CelestialBody primaryBody;

  66.     /** Celestial body with smaller mass, m2.*/
  67.     private final CelestialBody secondaryBody;

  69.     /** Simple constructor.
  70.      * @param primaryBody Primary body.
  71.      * @param secondaryBody Secondary body.
  72.      */
  73.     public L1TransformProvider(final CelestialBody primaryBody, final CelestialBody secondaryBody) {
  74.         this.primaryBody   = primaryBody;
  75.         this.secondaryBody = secondaryBody;
  76.         this.frame         = primaryBody.getInertiallyOrientedFrame();
  77.     }

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

  89.     /** {@inheritDoc} */
  90.     @Override
  91.     public StaticTransform getStaticTransform(final AbsoluteDate date) {
  92.         final PVCoordinates pv21        = secondaryBody.getPVCoordinates(date, frame);
  93.         final Vector3D      translation = getL1(pv21.getPosition()).negate();
  94.         final Rotation      rotation    = new Rotation(pv21.getPosition(), pv21.getVelocity(),
  95.                 Vector3D.PLUS_I, Vector3D.PLUS_J);
  96.         return StaticTransform.compose(
  97.                 date,
  98.                 StaticTransform.of(date, translation),
  99.                 StaticTransform.of(date, rotation));
  100.     }

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

  116.     /** {@inheritDoc} */
  117.     @Override
  118.     public <T extends CalculusFieldElement<T>> FieldStaticTransform<T> getStaticTransform(final FieldAbsoluteDate<T> date) {
  119.         final FieldPVCoordinates<T> pv21        = secondaryBody.getPVCoordinates(date, frame);
  120.         final FieldVector3D<T>      translation = getL1(pv21.getPosition()).negate();
  121.         final FieldRotation<T>      rotation    = new FieldRotation<>(pv21.getPosition(), pv21.getVelocity(),
  122.                 FieldVector3D.getPlusI(date.getField()), FieldVector3D.getPlusJ(date.getField()));
  123.         return FieldStaticTransform.compose(
  124.                 date,
  125.                 FieldStaticTransform.of(date, translation),
  126.                 FieldStaticTransform.of(date, rotation));
  127.     }

  129.     /** Compute the coordinates of the L1 point.
  130.      * @param primaryToSecondary relative position of secondary body with respect to primary body
  131.      * @return coordinates of the L1 point given in frame: primaryBody.getInertiallyOrientedFrame()
  132.      */
  133.     private Vector3D getL1(final Vector3D primaryToSecondary) {

  135.         // mass ratio
  136.         final double massRatio = secondaryBody.getGM() / primaryBody.getGM();

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

  142.         // Accurate position of L1 point, by solving the L1 equilibrium equation
  143.         final UnivariateFunction l1Equation = r -> {
  144.             final double bigrminusR  = bigR - r;
  145.             final double lhs         = 1.0 / ( r * r );
  146.             final double rhs1        = massRatio / (bigrminusR * bigrminusR);
  147.             final double rhs2        = 1.0 / (bigR * bigR);
  148.             final double rhs3        = (1 + massRatio) * bigrminusR * rhs2 / bigR;
  149.             return lhs - (rhs1 + rhs2 - rhs3);
  150.         };
  151.         final double[] searchInterval = UnivariateSolverUtils.bracket(l1Equation,
  152.                                                                       baseR, 0, bigR,
  153.                                                                       0.01 * bigR, 1, MAX_EVALUATIONS);
  154.         final BracketingNthOrderBrentSolver solver =
  155.                         new BracketingNthOrderBrentSolver(RELATIVE_ACCURACY,
  156.                                                           ABSOLUTE_ACCURACY,
  157.                                                           FUNCTION_ACCURACY,
  158.                                                           MAX_ORDER);
  159.         final double r = solver.solve(MAX_EVALUATIONS, l1Equation,
  160.                                       searchInterval[0], searchInterval[1],
  161.                                       AllowedSolution.ANY_SIDE);

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

  166.     }

  168.     /** Compute the coordinates of the L1 point.
  169.      * @param <T> type of the field elements
  170.      * @param primaryToSecondary relative position of secondary body with respect to primary body
  171.      * @return coordinates of the L1 point given in frame: primaryBody.getInertiallyOrientedFrame()
  172.      */
  173.     private <T extends CalculusFieldElement<T>> FieldVector3D<T>
  174.         getL1(final FieldVector3D<T> primaryToSecondary) {

  176.         // mass ratio
  177.         final double massRatio = secondaryBody.getGM() / primaryBody.getGM();

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

  183.         // Accurate position of L1 point, by solving the L1 equilibrium equation
  184.         final CalculusFieldUnivariateFunction<T> l1Equation = r -> {
  185.             final T bigrminusR = bigR.subtract(r);
  186.             final T lhs        = r.multiply(r).reciprocal();
  187.             final T rhs1       = bigrminusR.multiply(bigrminusR).reciprocal().multiply(massRatio);
  188.             final T rhs2       = bigR.multiply(bigR).reciprocal();
  189.             final T rhs3       = bigrminusR.multiply(rhs2).multiply(1 + massRatio).divide(bigR);
  190.             return lhs.subtract(rhs1.add(rhs2).add(rhs3));
  191.         };
  192.         final T zero             = primaryToSecondary.getX().getField().getZero();
  193.         final T[] searchInterval = UnivariateSolverUtils.bracket(l1Equation,
  194.                                                                  baseR, zero, bigR.multiply(2),
  195.                                                                  bigR.multiply(0.01), zero.newInstance(1),
  196.                                                                  MAX_EVALUATIONS);


  199.         final T relativeAccuracy = zero.newInstance(RELATIVE_ACCURACY);
  200.         final T absoluteAccuracy = zero.newInstance(ABSOLUTE_ACCURACY);
  201.         final T functionAccuracy = zero.newInstance(FUNCTION_ACCURACY);

  203.         final FieldBracketingNthOrderBrentSolver<T> solver =
  204.                         new FieldBracketingNthOrderBrentSolver<>(relativeAccuracy,
  205.                                                                  absoluteAccuracy,
  206.                                                                  functionAccuracy,
  207.                                                                  MAX_ORDER);
  208.         final T r = solver.solve(MAX_EVALUATIONS, l1Equation,
  209.                                  searchInterval[0], searchInterval[1],
  210.                                  AllowedSolution.ANY_SIDE);

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

  215.     }

  217. }
Created with JaCoCo 0.8.12.202403310830