1 /* Copyright 2002-2025 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.orbits;
18
19 import org.hipparchus.complex.Complex;
20 import org.hipparchus.geometry.euclidean.threed.Vector3D;
21 import org.hipparchus.linear.EigenDecompositionNonSymmetric;
22 import org.hipparchus.linear.FieldVector;
23 import org.hipparchus.linear.MatrixUtils;
24 import org.hipparchus.linear.RealMatrix;
25 import org.hipparchus.linear.RealVector;
26 import org.orekit.bodies.CR3BPSystem;
27 import org.orekit.propagation.SpacecraftState;
28 import org.orekit.propagation.numerical.cr3bp.STMEquations;
29 import org.orekit.utils.PVCoordinates;
30
31 /**
32 * Base class for libration orbits.
33 * @see HaloOrbit
34 * @see LyapunovOrbit
35 * @author Vincent Mouraux
36 * @author Bryan Cazabonne
37 * @since 10.2
38 */
39 public abstract class LibrationOrbit {
40
41 /** CR3BP System of the libration Orbit. */
42 private final CR3BPSystem syst;
43
44 /** Position-Velocity initial position on a libration Orbit. */
45 private PVCoordinates initialPV;
46
47 /** Orbital Period of the libration Orbit. */
48 private double orbitalPeriod;
49
50 /**
51 * Constructor.
52 * @param system CR3BP System considered
53 * @param initialPV initial position on a libration Orbit
54 * @param orbitalPeriod initial orbital period of the libration Orbit
55 */
56 protected LibrationOrbit(final CR3BPSystem system,
57 final PVCoordinates initialPV,
58 final double orbitalPeriod) {
59 this.syst = system;
60 this.initialPV = initialPV;
61 this.orbitalPeriod = orbitalPeriod;
62 }
63
64 /** Return the orbital period of the libration orbit.
65 * @return orbitalPeriod orbital period of the libration orbit
66 */
67 public double getOrbitalPeriod() {
68 return orbitalPeriod;
69 }
70
71 /** Return the initialPV on the libration orbit.
72 * <p>
73 * This will return the exact initialPV only if you applied a prior
74 * differential correction. If you did not, you can use the method
75 * {@link #applyCorrectionOnPV(CR3BPDifferentialCorrection)}
76 * </p>
77 * @return initialPV initialPV on the libration orbit
78 */
79 public PVCoordinates getInitialPV() {
80 return initialPV;
81 }
82
83 /** Apply differential correction.
84 * <p>
85 * This will update {@link #initialPV} and
86 * {@link #orbitalPeriod} parameters.
87 * </p>
88 */
89 public void applyDifferentialCorrection() {
90 final CR3BPDifferentialCorrection diff = new CR3BPDifferentialCorrection(initialPV, syst, orbitalPeriod);
91 initialPV = applyCorrectionOnPV(diff);
92 orbitalPeriod = diff.getOrbitalPeriod();
93 }
94
95 /** Return a manifold direction from one position on a libration Orbit.
96 * @param s SpacecraftState with additional equations
97 * @param isStable true if the manifold is stable
98 * @return manifold first guess Position-Velocity of a point on the libration Orbit
99 */
100 public PVCoordinates getManifolds(final SpacecraftState s, final boolean isStable) {
101
102 // Get the index of the eigen vector of the state transition matrix,
103 // depending on the stability or unstability of the manifold
104 final int eigenVectorIndex = isStable ? 1 : 0;
105
106 // Small delta, linked to the characteristic velocity of the CR3BP system
107 final double epsilon = syst.getVdim() * 1E2 / syst.getDdim();
108
109 // Get monodromy (i.e. state transition) matrix and its eigen decomposition
110 final RealMatrix phi = new STMEquations(syst).getStateTransitionMatrix(s);
111 final EigenDecompositionNonSymmetric eigen = new EigenDecompositionNonSymmetric(phi);
112
113 // Get normalized eigen vector linked to the stability of the manifold
114 final FieldVector<Complex> cv = eigen.getEigenvector(eigenVectorIndex);
115
116 // Get real vector value and normalize
117 final RealVector rv = MatrixUtils.createRealVector(cv.getDimension());
118 for (int i = 0; i < cv.getDimension(); ++i) {
119 rv.setEntry(i, cv.getEntry(i).getRealPart());
120 }
121 final RealVector eigenVector = rv.unitVector();
122
123 // New PVCoordinates following the manifold
124 return new PVCoordinates(s.getPosition().add(new Vector3D(eigenVector.getEntry(0),
125 eigenVector.getEntry(1),
126 eigenVector.getEntry(2)).scalarMultiply(epsilon)),
127 s.getPVCoordinates().getVelocity().add(new Vector3D(eigenVector.getEntry(3),
128 eigenVector.getEntry(4),
129 eigenVector.getEntry(5)).scalarMultiply(epsilon)));
130 }
131
132 /**
133 * Apply the differential correction to compute more accurate initial PV.
134 * @param diff cr3bp differential correction
135 * @return corrected PV coordinates
136 */
137 protected abstract PVCoordinates applyCorrectionOnPV(CR3BPDifferentialCorrection diff);
138
139 }