CR3BPSystem.java
/* Copyright 2002-2023 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.bodies;
import org.hipparchus.analysis.UnivariateFunction;
import org.hipparchus.analysis.solvers.AllowedSolution;
import org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver;
import org.hipparchus.analysis.solvers.UnivariateSolverUtils;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.orekit.frames.CR3BPRotatingFrame;
import org.orekit.frames.Frame;
import org.orekit.frames.Transform;
import org.orekit.time.AbsoluteDate;
import org.orekit.utils.AbsolutePVCoordinates;
import org.orekit.utils.LagrangianPoints;
import org.orekit.utils.PVCoordinates;
/**
* Class creating, from two different celestial bodies, the corresponding system
* with respect to the Circular Restricted Three Body problem hypotheses.
* @see "Dynamical systems, the three-body problem, and space mission design, Koon, Lo, Marsden, Ross"
* @author Vincent Mouraux
* @author William Desprats
* @since 10.2
*/
public class CR3BPSystem {
/** 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 accuracy 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 = 10000;
/** Mass ratio. */
private final double mu;
/** Distance between the two primaries, meters. */
private final double dDim;
/** Orbital Velocity of m1, m/s. */
private final double vDim;
/** Orbital Period of m1 and m2, seconds. */
private final double tDim;
/** CR3BP System name. */
private final String name;
/** Rotating Frame for the system. */
private final Frame rotatingFrame;
/** Primary body. */
private final CelestialBody primaryBody;
/** Secondary body. */
private final CelestialBody secondaryBody;
/** L1 Point position in the rotating frame. */
private Vector3D l1Position;
/** L2 Point position in the rotating frame. */
private Vector3D l2Position;
/** L3 Point position in the rotating frame. */
private Vector3D l3Position;
/** L4 Point position in the rotating frame. */
private Vector3D l4Position;
/** L5 Point position in the rotating frame. */
private Vector3D l5Position;
/** Distance between a L1 and the secondaryBody. */
private final double gamma1;
/** Distance between a L2 and the secondaryBody. */
private final double gamma2;
/** Distance between a L3 and the primaryBody. */
private final double gamma3;
/**
* Simple constructor.
* <p>
* Standard constructor to use to define a CR3BP System. Mass ratio is
* calculated from JPL data.
* </p>
* @param primaryBody Primary Body in the CR3BP System
* @param secondaryBody Secondary Body in the CR3BP System
* @param a Semi-Major Axis of the secondary body
*/
public CR3BPSystem(final CelestialBody primaryBody, final CelestialBody secondaryBody, final double a) {
this(primaryBody, secondaryBody, a, secondaryBody.getGM() / (secondaryBody.getGM() + primaryBody.getGM()));
}
/**
* Simple constructor.
* <p>
* This constructor can be used to define a CR3BP System if the user wants
* to use a specified mass ratio.
* </p>
* @param primaryBody Primary Body in the CR3BP System
* @param secondaryBody Secondary Body in the CR3BP System
* @param a Semi-Major Axis of the secondary body
* @param mu Mass ratio of the CR3BPSystem
*/
public CR3BPSystem(final CelestialBody primaryBody, final CelestialBody secondaryBody, final double a, final double mu) {
this.primaryBody = primaryBody;
this.secondaryBody = secondaryBody;
this.name = primaryBody.getName() + "_" + secondaryBody.getName();
final double mu1 = primaryBody.getGM();
this.mu = mu;
this.dDim = a;
this.vDim = FastMath.sqrt(mu1 / dDim);
this.tDim = 2 * FastMath.PI * dDim / vDim;
this.rotatingFrame = new CR3BPRotatingFrame(mu, primaryBody, secondaryBody);
computeLagrangianPointsPosition();
// Calculation of collinear points gamma
// L1 Gamma
final double x1 = l1Position.getX();
final double dCP1 = 1 - mu;
this.gamma1 = dCP1 - x1;
// L2 Gamma
final double x2 = l2Position.getX();
final double dCP2 = 1 - mu;
this.gamma2 = x2 - dCP2;
// L3 Gamma
final double x3 = l3Position.getX();
final double dCP3 = -mu;
this.gamma3 = dCP3 - x3;
}
/** Calculation of Lagrangian Points position using CR3BP equations.
*/
private void computeLagrangianPointsPosition() {
// Calculation of Lagrangian Points position using CR3BP equations
// L1
final BracketingNthOrderBrentSolver solver =
new BracketingNthOrderBrentSolver(RELATIVE_ACCURACY,
ABSOLUTE_ACCURACY,
FUNCTION_ACCURACY, MAX_ORDER);
final double baseR1 = 1 - FastMath.cbrt(mu / 3);
final UnivariateFunction l1Equation = x -> {
final double leq1 =
x * (x + mu) * (x + mu) * (x + mu - 1) * (x + mu - 1);
final double leq2 = (1 - mu) * (x + mu - 1) * (x + mu - 1);
final double leq3 = mu * (x + mu) * (x + mu);
return leq1 - leq2 + leq3;
};
final double[] searchInterval1 =
UnivariateSolverUtils.bracket(l1Equation, baseR1, -mu,
1 - mu, 1E-6, 1,
MAX_EVALUATIONS);
final double r1 =
solver.solve(MAX_EVALUATIONS, l1Equation, searchInterval1[0],
searchInterval1[1], AllowedSolution.ANY_SIDE);
this.l1Position = new Vector3D(r1, 0.0, 0.0);
// L2
final double baseR2 = 1 + FastMath.cbrt(mu / 3);
final UnivariateFunction l2Equation = x -> {
final double leq21 =
x * (x + mu) * (x + mu) * (x + mu - 1) * (x + mu - 1);
final double leq22 = (1 - mu) * (x + mu - 1) * (x + mu - 1);
final double leq23 = mu * (x + mu) * (x + mu);
return leq21 - leq22 - leq23;
};
final double[] searchInterval2 =
UnivariateSolverUtils.bracket(l2Equation, baseR2, 1 - mu, 2, 1E-6,
1, MAX_EVALUATIONS);
final double r2 =
solver.solve(MAX_EVALUATIONS, l2Equation, searchInterval2[0],
searchInterval2[1], AllowedSolution.ANY_SIDE);
this.l2Position = new Vector3D(r2, 0.0, 0.0);
// L3
final double baseR3 = -(1 + 5 * mu / 12);
final UnivariateFunction l3Equation = x -> {
final double leq31 =
x * (x + mu) * (x + mu) * (x + mu - 1) * (x + mu - 1);
final double leq32 = (1 - mu) * (x + mu - 1) * (x + mu - 1);
final double leq33 = mu * (x + mu) * (x + mu);
return leq31 + leq32 + leq33;
};
final double[] searchInterval3 =
UnivariateSolverUtils.bracket(l3Equation, baseR3, -2, -mu, 1E-6, 1,
MAX_EVALUATIONS);
final double r3 =
solver.solve(MAX_EVALUATIONS, l3Equation, searchInterval3[0],
searchInterval3[1], AllowedSolution.ANY_SIDE);
this.l3Position = new Vector3D(r3, 0.0, 0.0);
// L4
this.l4Position = new Vector3D(0.5 - mu, FastMath.sqrt(3) / 2, 0);
// L5
this.l5Position = new Vector3D(0.5 - mu, -FastMath.sqrt(3) / 2, 0);
}
/** Get the CR3BP mass ratio of the system mu2/(mu1+mu2).
* @return CR3BP mass ratio of the system mu2/(mu1+mu2)
*/
public double getMassRatio() {
return mu;
}
/** Get the CR3BP distance between the two bodies.
* @return CR3BP distance between the two bodies(m)
*/
public double getDdim() {
return dDim;
}
/** Get the CR3BP orbital velocity of m2.
* @return CR3BP orbital velocity of m2(m/s)
*/
public double getVdim() {
return vDim;
}
/** Get the CR3BP orbital period of m2 around m1.
* @return CR3BP orbital period of m2 around m1(s)
*/
public double getTdim() {
return tDim;
}
/** Get the name of the CR3BP system.
* @return name of the CR3BP system
*/
public String getName() {
return name;
}
/** Get the primary CelestialBody.
* @return primary CelestialBody
*/
public CelestialBody getPrimary() {
return primaryBody;
}
/** Get the secondary CelestialBody.
* @return secondary CelestialBody
*/
public CelestialBody getSecondary() {
return secondaryBody;
}
/** Get the CR3BP Rotating Frame.
* @return CR3BP Rotating Frame
*/
public Frame getRotatingFrame() {
return rotatingFrame;
}
/** Get the position of the Lagrangian point in the CR3BP Rotating frame.
* @param lagrangianPoint Lagrangian Point to consider
* @return position of the Lagrangian point in the CR3BP Rotating frame (-)
*/
public Vector3D getLPosition(final LagrangianPoints lagrangianPoint) {
final Vector3D lPosition;
switch (lagrangianPoint) {
case L1:
lPosition = l1Position;
break;
case L3:
lPosition = l3Position;
break;
case L4:
lPosition = l4Position;
break;
case L5:
lPosition = l5Position;
break;
default:
lPosition = l2Position;
break;
}
return lPosition;
}
/**
* Get the position of the Lagrangian point in the CR3BP Rotating frame.
* @param lagrangianPoint Lagrangian Point to consider
* @return Distance between a Lagrangian Point and its closest primary.
*/
public double getGamma(final LagrangianPoints lagrangianPoint) {
final double gamma;
switch (lagrangianPoint) {
case L1:
gamma = gamma1;
break;
case L2:
gamma = gamma2;
break;
case L3:
gamma = gamma3;
break;
default:
gamma = 0;
}
return gamma;
}
/** Get the PVCoordinates from normalized units to standard units in an output frame.
* @param pv0 Normalized PVCoordinates in the rotating frame
* @param date Date of the transform
* @param outputFrame Frame in which the output PVCoordinates will be
* @return PVCoordinates in the output frame [m,m/s]
*/
private PVCoordinates getRealPV(final PVCoordinates pv0, final AbsoluteDate date, final Frame outputFrame) {
// 1. Dimensionalize the primary-centered rotating state using the instantaneously
// defined characteristic quantities
// 2. Apply the transformation to primary inertial frame
// 3. Apply the transformation to output frame
final Frame primaryInertialFrame = primaryBody.getInertiallyOrientedFrame();
final Vector3D pv21 = secondaryBody.getPosition(date, primaryInertialFrame);
// Distance and Velocity to dimensionalize the state vector
final double dist12 = pv21.getNorm();
final double vCircular = FastMath.sqrt(primaryBody.getGM() / dist12);
// Dimensionalized state vector centered on primary body
final PVCoordinates pvDim = new PVCoordinates(pv0.getPosition().scalarMultiply(dist12),
pv0.getVelocity().scalarMultiply(vCircular));
// Transformation between rotating frame and the primary inertial
final Transform rotatingToPrimaryInertial = rotatingFrame.getTransformTo(primaryInertialFrame, date);
// State vector in the primary inertial frame
final PVCoordinates pv2 = rotatingToPrimaryInertial.transformPVCoordinates(pvDim);
// Transformation between primary inertial frame and the output frame
final Transform primaryInertialToOutputFrame = primaryInertialFrame.getTransformTo(outputFrame, date);
return primaryInertialToOutputFrame.transformPVCoordinates(pv2);
}
/** Get the AbsolutePVCoordinates from normalized units to standard units in an output frame.
* This method ensure the constituency of the date of returned AbsolutePVCoordinate, especially
* when apv0 is the result of a propagation in CR3BP normalized model.
* @param apv0 Normalized AbsolutePVCoordinates in the rotating frame
* @param initialDate Date of the at the beginning of the propagation
* @param outputFrame Frame in which the output AbsolutePVCoordinates will be
* @return AbsolutePVCoordinates in the output frame [m,m/s]
*/
public AbsolutePVCoordinates getRealAPV(final AbsolutePVCoordinates apv0, final AbsoluteDate initialDate, final Frame outputFrame) {
final double duration = apv0.getDate().durationFrom(initialDate) * tDim / (2 * FastMath.PI);
final AbsoluteDate date = initialDate.shiftedBy(duration);
// PVCoordinate in the output frame
final PVCoordinates pv3 = getRealPV(apv0.getPVCoordinates(), date, outputFrame);
return new AbsolutePVCoordinates(outputFrame, date, pv3);
}
}