L2TransformProvider.java

  1. /* Copyright 2002-2022 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 <T extends CalculusFieldElement<T>> FieldTransform<T> getTransform(final FieldAbsoluteDate<T> date) {
  95.         final FieldPVCoordinates<T> pv21        = secondaryBody.getPVCoordinates(date, frame);
  96.         final FieldVector3D<T>      translation = getL2(pv21.getPosition()).negate();
  97.         final Field<T>              field       = pv21.getPosition().getX().getField();
  98.         final FieldRotation<T>      rotation    = new FieldRotation<>(pv21.getPosition(), pv21.getVelocity(),
  99.                                                                       FieldVector3D.getPlusI(field),
  100.                                                                       FieldVector3D.getPlusJ(field));
  101.         return new FieldTransform<T>(date,
  102.                                      new FieldTransform<>(date, translation),
  103.                                      new FieldTransform<>(date, rotation));
  104.     }

  106.     /** Compute the coordinates of the L2 point.
  107.      * @param primaryToSecondary relative position of secondary body with respect to primary body
  108.      * @return coordinates of the L2 point given in frame: primaryBody.getInertiallyOrientedFrame()
  109.      */
  110.     private Vector3D getL2(final Vector3D primaryToSecondary) {

  112.         // mass ratio
  113.         final double massRatio = secondaryBody.getGM() / primaryBody.getGM();

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

  119.         // Accurate position of L2 point, by solving the L2 equilibrium equation
  120.         final UnivariateFunction l2Equation = r -> {
  121.             final double rminusbigR  = r - bigR;
  122.             final double lhs1        = 1.0 / (r * r);
  123.             final double lhs2        = massRatio / (rminusbigR * rminusbigR);
  124.             final double rhs1        = 1.0 / (bigR * bigR);
  125.             final double rhs2        = (1 + massRatio) * rminusbigR * rhs1 / bigR;
  126.             return (lhs1 + lhs2) - (rhs1 + rhs2);
  127.         };
  128.         final double[] searchInterval = UnivariateSolverUtils.bracket(l2Equation,
  129.                                                                       baseR, 0, 2 * bigR,
  130.                                                                       0.01 * bigR, 1, MAX_EVALUATIONS);
  131.         final BracketingNthOrderBrentSolver solver =
  132.                         new BracketingNthOrderBrentSolver(RELATIVE_ACCURACY,
  133.                                                           ABSOLUTE_ACCURACY,
  134.                                                           FUNCTION_ACCURACY,
  135.                                                           MAX_ORDER);
  136.         final double r = solver.solve(MAX_EVALUATIONS, l2Equation,
  137.                                       searchInterval[0], searchInterval[1],
  138.                                       AllowedSolution.ANY_SIDE);

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

  143.     }

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

  153.         // mass ratio
  154.         final double massRatio = secondaryBody.getGM() / primaryBody.getGM();

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

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


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

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

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

  192.     }

  194. }
Created with JaCoCo 0.8.7.202105040129