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