CR3BPSystem.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.bodies;

  18. import org.hipparchus.analysis.UnivariateFunction;
  19. import org.hipparchus.analysis.solvers.AllowedSolution;
  20. import org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver;
  21. import org.hipparchus.analysis.solvers.UnivariateSolverUtils;
  22. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  23. import org.hipparchus.util.FastMath;
  24. import org.orekit.frames.CR3BPRotatingFrame;
  25. import org.orekit.frames.Frame;
  26. import org.orekit.frames.Transform;
  27. import org.orekit.time.AbsoluteDate;
  28. import org.orekit.utils.AbsolutePVCoordinates;
  29. import org.orekit.utils.LagrangianPoints;
  30. import org.orekit.utils.PVCoordinates;

  31. /**
  32.  * Class creating, from two different celestial bodies, the corresponding system
  33.  * with respect to the Circular Restricted Three Body problem hypotheses.
  34.  * @see "Dynamical systems, the three-body problem, and space mission design, Koon, Lo, Marsden, Ross"
  35.  * @author Vincent Mouraux
  36.  * @author William Desprats
  37.  * @since 10.2
  38.  */
  39. public class CR3BPSystem {

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

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

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

  46.     /** Maximal order for solver. */
  47.     private static final int MAX_ORDER = 5;

  48.     /** Maximal number of evaluations for solver. */
  49.     private static final int MAX_EVALUATIONS = 10000;

  50.     /** Mass ratio. */
  51.     private final double mu;

  52.     /** Distance between the two primaries, meters. */
  53.     private final double dDim;

  54.     /** Orbital Velocity of m1, m/s. */
  55.     private final double vDim;

  56.     /** Orbital Period of m1 and m2, seconds. */
  57.     private final double tDim;

  58.     /** CR3BP System name. */
  59.     private final String name;

  60.     /** Rotating Frame for the system. */
  61.     private final Frame rotatingFrame;

  62.     /** Primary body. */
  63.     private final CelestialBody primaryBody;

  64.     /** Secondary body. */
  65.     private final CelestialBody secondaryBody;

  66.     /** L1 Point position in the rotating frame. */
  67.     private Vector3D l1Position;

  68.     /** L2 Point position in the rotating frame. */
  69.     private Vector3D l2Position;

  70.     /** L3 Point position in the rotating frame. */
  71.     private Vector3D l3Position;

  72.     /** L4 Point position in the rotating frame. */
  73.     private Vector3D l4Position;

  74.     /** L5 Point position in the rotating frame. */
  75.     private Vector3D l5Position;

  76.     /** Distance between a L1 and the secondaryBody. */
  77.     private final double gamma1;

  78.     /** Distance between a L2 and the secondaryBody. */
  79.     private final double gamma2;

  80.     /** Distance between a L3 and the primaryBody. */
  81.     private final double gamma3;

  82.     /**
  83.      * Simple constructor.
  84.      * <p>
  85.      * Standard constructor to use to define a CR3BP System. Mass ratio is
  86.      * calculated from JPL data.
  87.      * </p>
  88.      * @param primaryBody Primary Body in the CR3BP System
  89.      * @param secondaryBody Secondary Body in the CR3BP System
  90.      * @param a Semi-Major Axis of the secondary body
  91.      */
  92.     public CR3BPSystem(final CelestialBody primaryBody, final CelestialBody secondaryBody, final double a) {
  93.         this(primaryBody, secondaryBody, a, secondaryBody.getGM() / (secondaryBody.getGM() + primaryBody.getGM()));
  94.     }

  95.     /**
  96.      * Simple constructor.
  97.      * <p>
  98.      * This constructor can be used to define a CR3BP System if the user wants
  99.      * to use a specified mass ratio.
  100.      * </p>
  101.      * @param primaryBody Primary Body in the CR3BP System
  102.      * @param secondaryBody Secondary Body in the CR3BP System
  103.      * @param a Semi-Major Axis of the secondary body
  104.      * @param mu Mass ratio of the CR3BPSystem
  105.      */
  106.     public CR3BPSystem(final CelestialBody primaryBody, final CelestialBody secondaryBody, final double a, final double mu) {
  107.         this.primaryBody = primaryBody;
  108.         this.secondaryBody = secondaryBody;

  109.         this.name = primaryBody.getName() + "_" + secondaryBody.getName();

  110.         final double mu1 = primaryBody.getGM();

  111.         this.mu = mu;
  112.         this.dDim = a;
  113.         this.vDim = FastMath.sqrt(mu1 / dDim);
  114.         this.tDim = 2 * FastMath.PI * dDim / vDim;
  115.         this.rotatingFrame = new CR3BPRotatingFrame(mu, primaryBody, secondaryBody);

  116.         computeLagrangianPointsPosition();

  117.         // Calculation of collinear points gamma

  118.         // L1 Gamma
  119.         final double x1 = l1Position.getX();
  120.         final double dCP1 = 1 - mu;
  121.         this.gamma1 = dCP1 - x1;

  122.         // L2 Gamma
  123.         final double x2 = l2Position.getX();
  124.         final double dCP2 = 1 - mu;
  125.         this.gamma2 = x2 - dCP2;

  126.         // L3 Gamma
  127.         final double x3 = l3Position.getX();
  128.         final double dCP3 = -mu;
  129.         this.gamma3 = dCP3 - x3;

  130.     }

  131.     /** Calculation of Lagrangian Points position using CR3BP equations.
  132.      */
  133.     private void computeLagrangianPointsPosition() {
  134.         // Calculation of Lagrangian Points position using CR3BP equations

  135.         // L1
  136.         final BracketingNthOrderBrentSolver solver =
  137.                         new BracketingNthOrderBrentSolver(RELATIVE_ACCURACY,
  138.                                                           ABSOLUTE_ACCURACY,
  139.                                                           FUNCTION_ACCURACY, MAX_ORDER);

  140.         final double baseR1 = 1 - FastMath.cbrt(mu / 3);
  141.         final UnivariateFunction l1Equation = x -> {
  142.             final double leq1 =
  143.                             x * (x + mu) * (x + mu) * (x + mu - 1) * (x + mu - 1);
  144.             final double leq2 = (1 - mu) * (x + mu - 1) * (x + mu - 1);
  145.             final double leq3 = mu * (x + mu) * (x + mu);
  146.             return leq1 - leq2 + leq3;
  147.         };
  148.         final double[] searchInterval1 =
  149.                         UnivariateSolverUtils.bracket(l1Equation, baseR1, -mu,
  150.                                                       1 - mu, 1E-6, 1,
  151.                                                       MAX_EVALUATIONS);

  152.         final double r1 =
  153.                         solver.solve(MAX_EVALUATIONS, l1Equation, searchInterval1[0],
  154.                                      searchInterval1[1], AllowedSolution.ANY_SIDE);

  155.         this.l1Position = new Vector3D(r1, 0.0, 0.0);

  156.         // L2
  157.         final double baseR2 = 1 + FastMath.cbrt(mu / 3);
  158.         final UnivariateFunction l2Equation = x -> {
  159.             final double leq21 =
  160.                             x * (x + mu) * (x + mu) * (x + mu - 1) * (x + mu - 1);
  161.             final double leq22 = (1 - mu) * (x + mu - 1) * (x + mu - 1);
  162.             final double leq23 = mu * (x + mu) * (x + mu);
  163.             return leq21 - leq22 - leq23;
  164.         };
  165.         final double[] searchInterval2 =
  166.                         UnivariateSolverUtils.bracket(l2Equation, baseR2, 1 - mu, 2, 1E-6,
  167.                                                       1, MAX_EVALUATIONS);

  168.         final double r2 =
  169.                         solver.solve(MAX_EVALUATIONS, l2Equation, searchInterval2[0],
  170.                                      searchInterval2[1], AllowedSolution.ANY_SIDE);

  171.         this.l2Position = new Vector3D(r2, 0.0, 0.0);

  172.         // L3
  173.         final double baseR3 = -(1 + 5 * mu / 12);
  174.         final UnivariateFunction l3Equation = x -> {
  175.             final double leq31 =
  176.                             x * (x + mu) * (x + mu) * (x + mu - 1) * (x + mu - 1);
  177.             final double leq32 = (1 - mu) * (x + mu - 1) * (x + mu - 1);
  178.             final double leq33 = mu * (x + mu) * (x + mu);
  179.             return leq31 + leq32 + leq33;
  180.         };
  181.         final double[] searchInterval3 =
  182.                         UnivariateSolverUtils.bracket(l3Equation, baseR3, -2, -mu, 1E-6, 1,
  183.                                                       MAX_EVALUATIONS);

  184.         final double r3 =
  185.                         solver.solve(MAX_EVALUATIONS, l3Equation, searchInterval3[0],
  186.                                      searchInterval3[1], AllowedSolution.ANY_SIDE);

  187.         this.l3Position = new Vector3D(r3, 0.0, 0.0);

  188.         // L4
  189.         this.l4Position = new Vector3D(0.5 - mu, FastMath.sqrt(3) / 2, 0);

  190.         // L5
  191.         this.l5Position = new Vector3D(0.5 - mu, -FastMath.sqrt(3) / 2, 0);
  192.     }

  193.     /** Get the CR3BP mass ratio of the system mu2/(mu1+mu2).
  194.      * @return CR3BP mass ratio of the system mu2/(mu1+mu2)
  195.      */
  196.     public double getMassRatio() {
  197.         return mu;
  198.     }

  199.     /** Get the CR3BP distance between the two bodies.
  200.      * @return CR3BP distance between the two bodies(m)
  201.      */
  202.     public double getDdim() {
  203.         return dDim;
  204.     }

  205.     /** Get the CR3BP orbital velocity of m2.
  206.      * @return CR3BP orbital velocity of m2(m/s)
  207.      */
  208.     public double getVdim() {
  209.         return vDim;
  210.     }

  211.     /** Get the CR3BP orbital period of m2 around m1.
  212.      * @return CR3BP orbital period of m2 around m1(s)
  213.      */
  214.     public double getTdim() {
  215.         return tDim;
  216.     }

  217.     /** Get the name of the CR3BP system.
  218.      * @return name of the CR3BP system
  219.      */
  220.     public String getName() {
  221.         return name;
  222.     }

  223.     /** Get the primary CelestialBody.
  224.      * @return primary CelestialBody
  225.      */
  226.     public CelestialBody getPrimary() {
  227.         return primaryBody;
  228.     }

  229.     /** Get the secondary CelestialBody.
  230.      * @return secondary CelestialBody
  231.      */
  232.     public CelestialBody getSecondary() {
  233.         return secondaryBody;
  234.     }

  235.     /** Get the CR3BP Rotating Frame.
  236.      * @return CR3BP Rotating Frame
  237.      */
  238.     public Frame getRotatingFrame() {
  239.         return rotatingFrame;
  240.     }

  241.     /** Get the position of the Lagrangian point in the CR3BP Rotating frame.
  242.      * @param lagrangianPoint Lagrangian Point to consider
  243.      * @return position of the Lagrangian point in the CR3BP Rotating frame (-)
  244.      */
  245.     public Vector3D getLPosition(final LagrangianPoints lagrangianPoint) {

  246.         final Vector3D lPosition;

  247.         switch (lagrangianPoint) {

  248.             case L1:
  249.                 lPosition = l1Position;
  250.                 break;

  251.             case L3:
  252.                 lPosition = l3Position;
  253.                 break;

  254.             case L4:
  255.                 lPosition = l4Position;
  256.                 break;

  257.             case L5:
  258.                 lPosition = l5Position;
  259.                 break;

  260.             default:
  261.                 lPosition = l2Position;
  262.                 break;

  263.         }
  264.         return lPosition;
  265.     }

  266.     /**
  267.      * Get the position of the Lagrangian point in the CR3BP Rotating frame.
  268.      * @param lagrangianPoint Lagrangian Point to consider
  269.      * @return Distance between a Lagrangian Point and its closest primary.
  270.      */
  271.     public double getGamma(final LagrangianPoints lagrangianPoint) {

  272.         final double gamma;

  273.         switch (lagrangianPoint) {

  274.             case L1:
  275.                 gamma = gamma1;
  276.                 break;

  277.             case L2:
  278.                 gamma = gamma2;
  279.                 break;

  280.             case L3:
  281.                 gamma = gamma3;
  282.                 break;

  283.             default:
  284.                 gamma = 0;
  285.         }
  286.         return gamma;
  287.     }

  288.     /** Get the PVCoordinates from normalized units to standard units in an output frame.
  289.      * @param pv0 Normalized PVCoordinates in the rotating frame
  290.      * @param date Date of the transform
  291.      * @param outputFrame Frame in which the output PVCoordinates will be
  292.      * @return PVCoordinates in the output frame [m,m/s]
  293.      */
  294.     private PVCoordinates getRealPV(final PVCoordinates pv0, final AbsoluteDate date, final Frame outputFrame) {
  295.         // 1.   Dimensionalize  the  primary-centered  rotating  state  using  the  instantaneously
  296.         //      defined characteristic quantities
  297.         // 2.   Apply the transformation to primary inertial frame
  298.         // 3.   Apply the transformation to output frame

  299.         final Frame primaryInertialFrame = primaryBody.getInertiallyOrientedFrame();
  300.         final Vector3D pv21 = secondaryBody.getPosition(date, primaryInertialFrame);

  301.         // Distance and Velocity to dimensionalize the state vector
  302.         final double dist12 = pv21.getNorm();
  303.         final double vCircular  = FastMath.sqrt(primaryBody.getGM() / dist12);

  304.         // Dimensionalized state vector centered on primary body
  305.         final PVCoordinates pvDim = new PVCoordinates(pv0.getPosition().scalarMultiply(dist12),
  306.                                                       pv0.getVelocity().scalarMultiply(vCircular));

  307.         // Transformation between rotating frame and the primary inertial
  308.         final Transform rotatingToPrimaryInertial = rotatingFrame.getTransformTo(primaryInertialFrame, date);

  309.         // State vector in the primary inertial frame
  310.         final PVCoordinates pv2 = rotatingToPrimaryInertial.transformPVCoordinates(pvDim);


  311.         // Transformation between primary inertial frame and the output frame
  312.         final Transform primaryInertialToOutputFrame = primaryInertialFrame.getTransformTo(outputFrame, date);

  313.         return primaryInertialToOutputFrame.transformPVCoordinates(pv2);
  314.     }

  315.     /** Get the AbsolutePVCoordinates from normalized units to standard units in an output frame.
  316.      * This method ensure the constituency of the date of returned AbsolutePVCoordinate, especially
  317.      * when apv0 is the result of a propagation in CR3BP normalized model.
  318.      * @param apv0 Normalized AbsolutePVCoordinates in the rotating frame
  319.      * @param initialDate Date of the at the beginning of the propagation
  320.      * @param outputFrame Frame in which the output AbsolutePVCoordinates will be
  321.      * @return AbsolutePVCoordinates in the output frame [m,m/s]
  322.      */
  323.     public AbsolutePVCoordinates getRealAPV(final AbsolutePVCoordinates apv0, final AbsoluteDate initialDate, final Frame outputFrame) {

  324.         final double duration = apv0.getDate().durationFrom(initialDate) * tDim / (2 * FastMath.PI);
  325.         final AbsoluteDate date = initialDate.shiftedBy(duration);

  326.         // PVCoordinate in the output frame
  327.         final PVCoordinates pv3 = getRealPV(apv0.getPVCoordinates(), date, outputFrame);

  328.         return new AbsolutePVCoordinates(outputFrame, date, pv3);
  329.     }

  330. }