L2TransformProvider.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;

  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. /** L2 Transform provider for a frame on the L2 Lagrange point of two celestial bodies.
  39.  *
  40.  * @author Luc Maisonobe
  41.  * @author Julio Hernanz
  42.  */
  43. class L2TransformProvider 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.     /** Serializable UID.*/
  61.     private static final long serialVersionUID = 20170725L;

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

  66.     /** Celestial body with bigger mass, m1.*/
  67.     private final CelestialBody primaryBody;

  69.     /** Celestial body with smaller mass, m2.*/
  70.     private final CelestialBody secondaryBody;

  72.     /** Simple constructor.
  73.      * @param primaryBody Primary body.
  74.      * @param secondaryBody Secondary body.
  75.      */
  76.     L2TransformProvider(final CelestialBody primaryBody, final CelestialBody secondaryBody) {
  77.         this.primaryBody = primaryBody;
  78.         this.secondaryBody = secondaryBody;
  79.         this.frame = primaryBody.getInertiallyOrientedFrame();
  80.     }

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

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

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

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

  132.     /** Compute the coordinates of the L2 point.
  133.      * @param primaryToSecondary relative position of secondary body with respect to primary body
  134.      * @return coordinates of the L2 point given in frame: primaryBody.getInertiallyOrientedFrame()
  135.      */
  136.     private Vector3D getL2(final Vector3D primaryToSecondary) {

  138.         // mass ratio
  139.         final double massRatio = secondaryBody.getGM() / primaryBody.getGM();

  141.         // Approximate position of L2 point, valid when m2 << m1
  142.         final double bigR  = primaryToSecondary.getNorm();
  143.         final double baseR = bigR * (FastMath.cbrt(massRatio / 3) + 1);

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

  166.         // L2 point is built
  167.         return new Vector3D(r / bigR, primaryToSecondary);

  169.     }

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

  179.         // mass ratio
  180.         final double massRatio = secondaryBody.getGM() / primaryBody.getGM();

  182.         // Approximate position of L2 point, valid when m2 << m1
  183.         final T bigR  = primaryToSecondary.getNorm();
  184.         final T baseR = bigR.multiply(FastMath.cbrt(massRatio / 3) + 1);

  186.         // Accurate position of L2 point, by solving the L2 equilibrium equation
  187.         final CalculusFieldUnivariateFunction<T> l2Equation = r -> {
  188.             final T rminusbigR = r.subtract(bigR);
  189.             final T lhs1       = r.multiply(r).reciprocal();
  190.             final T lhs2       = rminusbigR.multiply(rminusbigR).reciprocal().multiply(massRatio);
  191.             final T rhs1       = bigR.multiply(bigR).reciprocal();
  192.             final T rhs2       = rminusbigR.multiply(rhs1).multiply(1 + massRatio).divide(bigR);
  193.             return lhs1.add(lhs2).subtract(rhs1.add(rhs2));
  194.         };
  195.         final T zero             = primaryToSecondary.getX().getField().getZero();
  196.         final T[] searchInterval = UnivariateSolverUtils.bracket(l2Equation,
  197.                                                                  baseR, zero, bigR.multiply(2),
  198.                                                                  bigR.multiply(0.01), zero,
  199.                                                                  MAX_EVALUATIONS);


  202.         final T relativeAccuracy = zero.newInstance(RELATIVE_ACCURACY);
  203.         final T absoluteAccuracy = zero.newInstance(ABSOLUTE_ACCURACY);
  204.         final T functionAccuracy = zero.newInstance(FUNCTION_ACCURACY);

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

  215.         // L2 point is built
  216.         return new FieldVector3D<>(r.divide(bigR), primaryToSecondary);

  218.     }

  220. }
Created with JaCoCo 0.8.12.202403310830