L1TransformProvider.java
- /* Copyright 2002-2024 CS GROUP
- * Licensed to CS GROUP (CS) under one or more
- * contributor license agreements. See the NOTICE file distributed with
- * this work for additional information regarding copyright ownership.
- * CS licenses this file to You under the Apache License, Version 2.0
- * (the "License"); you may not use this file except in compliance with
- * the License. You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
- package org.orekit.frames;
- import org.hipparchus.Field;
- import org.hipparchus.CalculusFieldElement;
- import org.hipparchus.analysis.CalculusFieldUnivariateFunction;
- import org.hipparchus.analysis.UnivariateFunction;
- import org.hipparchus.analysis.solvers.AllowedSolution;
- import org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver;
- import org.hipparchus.analysis.solvers.FieldBracketingNthOrderBrentSolver;
- import org.hipparchus.analysis.solvers.UnivariateSolverUtils;
- import org.hipparchus.geometry.euclidean.threed.FieldRotation;
- import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
- import org.hipparchus.geometry.euclidean.threed.Rotation;
- import org.hipparchus.geometry.euclidean.threed.Vector3D;
- import org.hipparchus.util.FastMath;
- import org.orekit.bodies.CelestialBody;
- import org.orekit.time.AbsoluteDate;
- import org.orekit.time.FieldAbsoluteDate;
- import org.orekit.utils.FieldPVCoordinates;
- import org.orekit.utils.PVCoordinates;
- /** L1 Transform provider for a frame on the L1 Lagrange point of two celestial bodies.
- *
- * @author Luc Maisonobe
- * @author Julio Hernanz
- */
- public class L1TransformProvider implements TransformProvider {
- /** Relative accuracy on position for solver. */
- private static final double RELATIVE_ACCURACY = 1e-14;
- /** Absolute accuracy on position for solver (1mm). */
- private static final double ABSOLUTE_ACCURACY = 1e-3;
- /** Function value ccuracy for solver (set to 0 so we rely only on position for convergence). */
- private static final double FUNCTION_ACCURACY = 0;
- /** Maximal order for solver. */
- private static final int MAX_ORDER = 5;
- /** Maximal number of evaluations for solver. */
- private static final int MAX_EVALUATIONS = 1000;
- /** Serializable UID.*/
- private static final long serialVersionUID = 20170824L;
- /** Frame for results. Always defined as primaryBody's inertially oriented frame.*/
- private final Frame frame;
- /** Celestial body with bigger mass, m1.*/
- private final CelestialBody primaryBody;
- /** Celestial body with smaller mass, m2.*/
- private final CelestialBody secondaryBody;
- /** Simple constructor.
- * @param primaryBody Primary body.
- * @param secondaryBody Secondary body.
- */
- public L1TransformProvider(final CelestialBody primaryBody, final CelestialBody secondaryBody) {
- this.primaryBody = primaryBody;
- this.secondaryBody = secondaryBody;
- this.frame = primaryBody.getInertiallyOrientedFrame();
- }
- /** {@inheritDoc} */
- @Override
- public Transform getTransform(final AbsoluteDate date) {
- final PVCoordinates pv21 = secondaryBody.getPVCoordinates(date, frame);
- final Vector3D translation = getL1(pv21.getPosition()).negate();
- final Rotation rotation = new Rotation(pv21.getPosition(), pv21.getVelocity(),
- Vector3D.PLUS_I, Vector3D.PLUS_J);
- return new Transform(date, new Transform(date, translation), new Transform(date, rotation));
- }
- /** {@inheritDoc} */
- @Override
- public StaticTransform getStaticTransform(final AbsoluteDate date) {
- final PVCoordinates pv21 = secondaryBody.getPVCoordinates(date, frame);
- final Vector3D translation = getL1(pv21.getPosition()).negate();
- final Rotation rotation = new Rotation(pv21.getPosition(), pv21.getVelocity(),
- Vector3D.PLUS_I, Vector3D.PLUS_J);
- return StaticTransform.compose(
- date,
- StaticTransform.of(date, translation),
- StaticTransform.of(date, rotation));
- }
- /** {@inheritDoc} */
- @Override
- public <T extends CalculusFieldElement<T>> FieldTransform<T> getTransform(final FieldAbsoluteDate<T> date) {
- final FieldPVCoordinates<T> pv21 = secondaryBody.getPVCoordinates(date, frame);
- final FieldVector3D<T> translation = getL1(pv21.getPosition()).negate();
- final Field<T> field = pv21.getPosition().getX().getField();
- final FieldRotation<T> rotation = new FieldRotation<>(pv21.getPosition(), pv21.getVelocity(),
- FieldVector3D.getPlusI(field),
- FieldVector3D.getPlusJ(field));
- return new FieldTransform<>(date,
- new FieldTransform<>(date, translation),
- new FieldTransform<>(date, rotation));
- }
- /** {@inheritDoc} */
- @Override
- public <T extends CalculusFieldElement<T>> FieldStaticTransform<T> getStaticTransform(final FieldAbsoluteDate<T> date) {
- final FieldPVCoordinates<T> pv21 = secondaryBody.getPVCoordinates(date, frame);
- final FieldVector3D<T> translation = getL1(pv21.getPosition()).negate();
- final FieldRotation<T> rotation = new FieldRotation<>(pv21.getPosition(), pv21.getVelocity(),
- FieldVector3D.getPlusI(date.getField()), FieldVector3D.getPlusJ(date.getField()));
- return FieldStaticTransform.compose(
- date,
- FieldStaticTransform.of(date, translation),
- FieldStaticTransform.of(date, rotation));
- }
- /** Compute the coordinates of the L1 point.
- * @param primaryToSecondary relative position of secondary body with respect to primary body
- * @return coordinates of the L1 point given in frame: primaryBody.getInertiallyOrientedFrame()
- */
- private Vector3D getL1(final Vector3D primaryToSecondary) {
- // mass ratio
- final double massRatio = secondaryBody.getGM() / primaryBody.getGM();
- // Approximate position of L1 point, valid when m2 << m1
- final double bigR = primaryToSecondary.getNorm();
- final double baseR = bigR * (1 - FastMath.cbrt(massRatio / 3));
- // Accurate position of L1 point, by solving the L1 equilibrium equation
- final UnivariateFunction l1Equation = r -> {
- final double bigrminusR = bigR - r;
- final double lhs = 1.0 / ( r * r );
- final double rhs1 = massRatio / (bigrminusR * bigrminusR);
- final double rhs2 = 1.0 / (bigR * bigR);
- final double rhs3 = (1 + massRatio) * bigrminusR * rhs2 / bigR;
- return lhs - (rhs1 + rhs2 - rhs3);
- };
- final double[] searchInterval = UnivariateSolverUtils.bracket(l1Equation,
- baseR, 0, bigR,
- 0.01 * bigR, 1, MAX_EVALUATIONS);
- final BracketingNthOrderBrentSolver solver =
- new BracketingNthOrderBrentSolver(RELATIVE_ACCURACY,
- ABSOLUTE_ACCURACY,
- FUNCTION_ACCURACY,
- MAX_ORDER);
- final double r = solver.solve(MAX_EVALUATIONS, l1Equation,
- searchInterval[0], searchInterval[1],
- AllowedSolution.ANY_SIDE);
- // L1 point is built
- return new Vector3D(r / bigR, primaryToSecondary);
- }
- /** Compute the coordinates of the L1 point.
- * @param <T> type of the field elements
- * @param primaryToSecondary relative position of secondary body with respect to primary body
- * @return coordinates of the L1 point given in frame: primaryBody.getInertiallyOrientedFrame()
- */
- private <T extends CalculusFieldElement<T>> FieldVector3D<T>
- getL1(final FieldVector3D<T> primaryToSecondary) {
- // mass ratio
- final double massRatio = secondaryBody.getGM() / primaryBody.getGM();
- // Approximate position of L1 point, valid when m2 << m1
- final T bigR = primaryToSecondary.getNorm();
- final T baseR = bigR.multiply(1 - FastMath.cbrt(massRatio / 3));
- // Accurate position of L1 point, by solving the L1 equilibrium equation
- final CalculusFieldUnivariateFunction<T> l1Equation = r -> {
- final T bigrminusR = bigR.subtract(r);
- final T lhs = r.multiply(r).reciprocal();
- final T rhs1 = bigrminusR.multiply(bigrminusR).reciprocal().multiply(massRatio);
- final T rhs2 = bigR.multiply(bigR).reciprocal();
- final T rhs3 = bigrminusR.multiply(rhs2).multiply(1 + massRatio).divide(bigR);
- return lhs.subtract(rhs1.add(rhs2).add(rhs3));
- };
- final T zero = primaryToSecondary.getX().getField().getZero();
- final T[] searchInterval = UnivariateSolverUtils.bracket(l1Equation,
- baseR, zero, bigR.multiply(2),
- bigR.multiply(0.01), zero.add(1),
- MAX_EVALUATIONS);
- final T relativeAccuracy = zero.add(RELATIVE_ACCURACY);
- final T absoluteAccuracy = zero.add(ABSOLUTE_ACCURACY);
- final T functionAccuracy = zero.add(FUNCTION_ACCURACY);
- final FieldBracketingNthOrderBrentSolver<T> solver =
- new FieldBracketingNthOrderBrentSolver<>(relativeAccuracy,
- absoluteAccuracy,
- functionAccuracy,
- MAX_ORDER);
- final T r = solver.solve(MAX_EVALUATIONS, l1Equation,
- searchInterval[0], searchInterval[1],
- AllowedSolution.ANY_SIDE);
- // L1 point is built
- return new FieldVector3D<>(r.divide(bigR), primaryToSecondary);
- }
- }