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.orbits;
18
19 import org.hipparchus.geometry.euclidean.threed.Vector3D;
20 import org.hipparchus.linear.EigenDecomposition;
21 import org.hipparchus.linear.RealMatrix;
22 import org.hipparchus.linear.RealVector;
23 import org.orekit.attitudes.AttitudeProvider;
24 import org.orekit.bodies.CR3BPSystem;
25 import org.orekit.propagation.SpacecraftState;
26 import org.orekit.propagation.numerical.cr3bp.STMEquations;
27 import org.orekit.time.TimeScale;
28 import org.orekit.utils.PVCoordinates;
29
30 /**
31 * Base class for libration orbits.
32 * @see HaloOrbit
33 * @see LyapunovOrbit
34 * @author Vincent Mouraux
35 * @author Bryan Cazabonne
36 * @since 10.2
37 */
38 public abstract class LibrationOrbit {
39
40 /** CR3BP System of the libration Orbit. */
41 private final CR3BPSystem syst;
42
43 /** Position-Velocity initial position on a libration Orbit. */
44 private PVCoordinates initialPV;
45
46 /** Orbital Period of the libration Orbit. */
47 private double orbitalPeriod;
48
49 /**
50 * Constructor.
51 * @param system CR3BP System considered
52 * @param initialPV initial position on a libration Orbit
53 * @param orbitalPeriod initial orbital period of the libration Orbit
54 */
55 protected LibrationOrbit(final CR3BPSystem system,
56 final PVCoordinates initialPV,
57 final double orbitalPeriod) {
58 this.syst = system;
59 this.initialPV = initialPV;
60 this.orbitalPeriod = orbitalPeriod;
61 }
62
63 /** Return the orbital period of the libration orbit.
64 * @return orbitalPeriod orbital period of the libration orbit
65 */
66 public double getOrbitalPeriod() {
67 return orbitalPeriod;
68 }
69
70 /** Return the initialPV on the libration orbit.
71 * <p>
72 * This will return the exact initialPV only if you applied a prior
73 * differential correction. If you did not, you can use the method
74 * {@link #applyCorrectionOnPV(CR3BPDifferentialCorrection)}
75 * </p>
76 * @return initialPV initialPV on the libration orbit
77 */
78 public PVCoordinates getInitialPV() {
79 return initialPV;
80 }
81
82 /** Apply differential correction.
83 * <p>
84 * This will update {@link #initialPV} and
85 * {@link #orbitalPeriod} parameters.
86 * </p>
87 */
88 public void applyDifferentialCorrection() {
89 final CR3BPDifferentialCorrection diff = new CR3BPDifferentialCorrection(initialPV, syst, orbitalPeriod);
90 initialPV = applyCorrectionOnPV(diff);
91 orbitalPeriod = diff.getOrbitalPeriod();
92 }
93
94 /** Apply differential correction.
95 * <p>
96 * This will update {@link #initialPV} and
97 * {@link #orbitalPeriod} parameters.
98 * </p>
99 * @param attitudeProvider the attitude law for the numerocal propagator
100 * @param utc UTC time scale
101 * @deprecated as of 11.1, replaced by {@link #applyDifferentialCorrection()}
102 */
103 @Deprecated
104 public void applyDifferentialCorrection(final AttitudeProvider attitudeProvider, final TimeScale utc) {
105 applyDifferentialCorrection();
106 }
107
108 /** Return a manifold direction from one position on a libration Orbit.
109 * @param s SpacecraftState with additional equations
110 * @param isStable true if the manifold is stable
111 * @return manifold first guess Position-Velocity of a point on the libration Orbit
112 */
113 public PVCoordinates getManifolds(final SpacecraftState s,
114 final boolean isStable) {
115 return isStable ? getStableManifolds(s) : getUnstableManifolds(s);
116 }
117
118 /** Return the stable manifold direction for several positions on a libration Orbit.
119 * @param s SpacecraftStates (with STM equations) to compute from
120 * @return Stable manifold first direction from a point on the libration Orbit
121 */
122 private PVCoordinates getStableManifolds(final SpacecraftState s) {
123
124 // Small delta, linked to the characteristic velocity of the CR3BP system
125 final double epsilon = syst.getVdim() * 1E2 / syst.getDdim();
126
127 // Get Normalize eigen vector linked to the stability of the manifold
128 final RealMatrix phi = new STMEquations(syst).getStateTransitionMatrix(s);
129 final RealVector eigenVector = new EigenDecomposition(phi).getEigenvector(1).unitVector();
130
131 // New PVCoordinates following the manifold
132 return new PVCoordinates(s.getPVCoordinates().getPosition()
133 .add(new Vector3D(eigenVector.getEntry(0), eigenVector
134 .getEntry(1), eigenVector.getEntry(2))
135 .scalarMultiply(epsilon)), s.getPVCoordinates()
136 .getVelocity()
137 .add(new Vector3D(eigenVector.getEntry(3),
138 eigenVector.getEntry(4),
139 eigenVector.getEntry(5))
140 .scalarMultiply(epsilon)));
141 }
142
143 /** Get the Unstable manifold direction for several positions on a libration Orbit.
144 * @param s spacecraft state (with STM equations) to compute from
145 * @return pv coordinates representing the unstable manifold first direction
146 * from a point on the libration Orbit
147 */
148 private PVCoordinates getUnstableManifolds(final SpacecraftState s) {
149
150 // Small delta, linked to the characteristic velocity of the CR3BP system
151 final double epsilon =
152 syst.getVdim() * 1E2 / syst.getDdim();
153
154 // Get Normalize eigen vector linked to the stability of the manifold
155 final RealMatrix phi = new STMEquations(syst).getStateTransitionMatrix(s);
156 final RealVector eigenVector = new EigenDecomposition(phi).getEigenvector(0).unitVector();
157
158 // New PVCoordinates following the manifold
159 return new PVCoordinates(s.getPVCoordinates().getPosition()
160 .add(new Vector3D(eigenVector.getEntry(0), eigenVector
161 .getEntry(1), eigenVector.getEntry(2))
162 .scalarMultiply(epsilon)), s.getPVCoordinates()
163 .getVelocity()
164 .add(new Vector3D(eigenVector.getEntry(3),
165 eigenVector.getEntry(4),
166 eigenVector.getEntry(5))
167 .scalarMultiply(epsilon)));
168 }
169
170 /**
171 * Apply the differential correction to compute more accurate initial PV.
172 * @param diff cr3bp differential correction
173 * @return corrected PV coordinates
174 */
175 protected abstract PVCoordinates applyCorrectionOnPV(CR3BPDifferentialCorrection diff);
176
177 }