3.3. Simulation for Sun-Earth CRTBP 47
Start
Set an initial guess and Transform the variable
Integrate until
Calculate variations &
Set new initial conditions No
Yes
Finish Retransform the variable
Center Manifold Design StepDifferential Correction Step
Propagate the sequences Calculate the first term of the sequences
Yes
No No
ε ε
ε' ε'
u s
Figure 3.2: Flowchart of center manifold design method for periodic orbits.
3.3. Simulation for Sun-Earth CRTBP 49
-2 -1 0
5 105
z [km]
1 2
-5
y [km] 10
0 5
x [km]
105
0 5 -5
Figure 3.3: Nominal halo orbit.
To verify the proposed method, the first to fourth elements of zn(t0) are set as the initial conditionξ:
ξ= [
0 −1.8058 −11.000 0 ]T
×10−3 (3.34)
The histories ofzc,k(t;ξ), zs,k(t;ξ), and zu,k(t;ξ) for k=0,1,5,10 by propagating the initial condition (3.32) are plotted in Fig. 3.4. It can be seen that zc,k(t;ξ), zs,k(t;ξ), and zu,k(t;ξ) converge to the nominal halo orbit as k increases. Consequently, the approximate initial conditionz(t0) =
[
ξT zs(t0;ξ) zu(t0;ξ) ]T
converges to the nominal initial conditionzn(t0) as k increases. Therefore, the center manifold design method can compute a libration point orbit by giving an initial parameter vectorξ.
3.3.2 Physical meaning of the parameter vector
All periodic or quasi-periodic orbits lying on the center manifold of the system are parameterized by a single four-dimensional vector ξ. From Eq. (3.11), only zc1 and zc2
affect the out-of-plane motion. Therefore, planar libration point orbits can be obtained by
0 1 2 3 4 5 6 t
-3 -2 -1 0 1 2 3
z c1 10-3
Nominal k=0 k=1 k=5 k=10
(a) zc1,k(t;ξ)
0 1 2 3 4 5 6
t -3
-2 -1 0 1 2 3
z c2 10-3
(b) zc2,k(t;ξ)
0 1 2 3 4 5 6
t -3
-2 -1 0 1 2 3
z c3 10-2
(c) zc3,k(t;ξ)
0 1 2 3 4 5 6
t -3
-2 -1 0 1 2 3
z c4 10-2
(d) zc4,k(t;ξ)
0 1 2 3 4 5 6
t -2
0 2 4 6 8 10
z s 10-4
(e) zs,k(t;ξ)
0 1 2 3 4 5 6
t -2
0 2 4 6 8 10
z u 10-4
(f) zu,k(t;ξ)
Figure 3.4: Time histories of the sequences (k=0,1,5,10) and the nominal halo orbit.
3.3. Simulation for Sun-Earth CRTBP 51
1 -5
0
z [km]
105
5
0
y [km] 10
5
5 0
x [km]
104
10 15 -1
(a) 3-D view
0 5 10 15
x [km] 104
-8 -6 -4 -2 0 2 4 6 8
z [km]
105
(b) projection onto x-z plane
Figure 3.5: Quasi-periodic orbit (quasi-vertical Lyapunov orbit): ξ= [1,0,0,0]T×10−2. The initial positions are shown by asterisks.
setting zc1(t0) =zc2(t0) =0. Otherwise, the orbits converge to three-dimensional periodic or quasi-periodic orbits.
As examples, the libration point orbits are computed by setting one element of the parameter vector to 0.01 and the remaining to 0. The Lagrangian point is fixed at L1. The orbits obtained by each parameter vector are shown in Figs. 3.5–3.7. When either zc1(t0) or zc2(t0) takes nonzero value, the obtained orbits are quasi-vertical Lyapunov orbits whose amplitudes in the z-direction are much higher than those in the x and y directions (Figs. 3.5 and 3.6).
By contrast, the orbits using only zc3(t0)or zc4(t0)are planar periodic orbits called Lyapunov orbits (Fig. 3.7).
3.3.3 Family of Lyapunov orbits
From the findings stated in Section 3.3.2, a Lyapunov orbit is found by a parameter vector ξLyapunov=
[
0 0 zc3(t0) zc4(t0) ]T
(3.35)
-5 0
5 105
z [km]
5
0
104 y [km]
5 0
x [km]
104
10 15 -5
(a) 3-D view
0 5 10 15
x [km] 104
-8 -6 -4 -2 0 2 4 6 8
z [km]
105
(b) projection onto x-z plane
Figure 3.6: Quasi-periodic orbit (quasi-vertical Lyapunov orbit): ξ= [0,1,0,0]T×10−2. The initial positions are shown by asterisks.
where at least either zc3(t0) or zc4(t0) is nonzero. For example, by setting zc4(t0) =0 and zc3(t0) = 2n×10−3 (n = 1,2,···,6), the L1 Lyapunov orbits and L2 Lyapunov orbits are obtained by the center manifold design step as shown in Fig. 3.8.
3.3.4 Finding vertical Lyapunov orbit from quasi-periodic orbit around L
1For a vertical Lyapunov orbit, the initial parameter vectorξis randomly chosen as ξ=
[
0 10 0 0 ]T
×10−3 (3.36)
This parameter vector generates a quasi-periodic orbit, which has already been shown in Fig. 3.6. By setting a Poincaré section Σ to the x-y plane and the one-sided Poincaré map as P+ including only the intersections of the quasi-periodic orbit with Σ that have positive values of ˙z, the discrete points onΣfall on a closed curve. Therefore, the quasi-periodic orbit can be characterized by two periods: the returning time Trand the winding time Tw. Define the
3.3. Simulation for Sun-Earth CRTBP 53
-5 0 5
x [km] 105
-8 -6 -4 -2 0 2 4 6 8
y [km]
105
(a) zc3(t0) =0.01 and zc4(t0) =0
-5 0 5
x [km] 105
-8 -6 -4 -2 0 2 4 6 8
y [km]
105
(b) zc3(t0) =0 and zc4(t0) =0.01
Figure 3.7: Periodic orbits (Lyapunov orbits): ξ= [0,0,zc3(t0),zc4(t0)]T. The initial positions are shown by asterisks.
-5 0 5
x [km] 105
-8 -6 -4 -2 0 2 4 6 8
y [km]
105
n=1 n=2 n=3 n=4 n=5 n=6
(a) L1Lyapunov family
-5 0 5
x [km] 105
-8 -6 -4 -2 0 2 4 6 8
y [km]
105
n=1 n=2 n=3 n=4 n=5 n=6
(b) L2Lyapunov family
Figure 3.8: Lyapunov families obtained by the center manifold design method.
returning time Trto be the average time taken for points starting fromΣto return back to it, i.e.
Tr= lim
n→∞
1 n
∑
n i=1(Ti+1−Ti) (3.37)
where Tiis the discrete time obtained by the ith iterate ofP+. On the other hand, according to the literature [105], the winding time represents the average number of iterates ofP+required to return back tox0. DefinePi to be the discrete point obtained by the ith iterate of P+ and assume thatPik passes overP1 after the discrete points move k times around the closed loop.
Then, the winding time is defined by
Tw= lim
k→∞
ik
k (3.38)
From Eqs. (3.37) and (3.38), the returning time and winding time of the obtained quasi-periodic orbit are computed as Tr=3.1714 and Tw=19.515, respectively. Thus, the initial guess of the half-period of the periodic orbit for the differential correction step is set as half of the returning time: tc=1.5857.
When the parameter vector ξ= [
0 zc2(t0) 0 0 ]T
, a purely periodic orbit never exists except for zc2(t0) =0; this is the Lagrangian point. Therefore, zc2 is fixed and zc3 is modified in the differential correction step, as is tc. The initial set(ξ,tc)is corrected by Eq. (3.30), and it finally converges to
ξ= [
0 10 0.61160 0 ]T
×10−3, tc=1.5855 (3.39) This corrected set (ξ,tc)forms the exact periodic orbit, which is the vertical Lyapunov orbit shown by the red line in Fig. 3.9. Note that the z-component after the differential correction step is the same as that of the quasi-vertical Lyapunov orbit. This is because fixing zc2 is equivalent to fixing z from Eq. (3.11).
3.3.5 Finding halo orbit from quasi-periodic orbit around L
2Because the proposed method is applicable to any collinear libration point, another application for a halo orbit around the L2point is demonstrated. The initial parameter vectorξ
3.3. Simulation for Sun-Earth CRTBP 55
-5 0
0 5
105
z [km]
5
5
104 y [km]
0 104
x [km]
10 15 -5
Vertical Lyapunov orbit Quasi vertical Lyapunov orbit
(a) 3-D view
0 5 10 15
x [km] 104
-8 -6 -4 -2 0 2 4 6 8
y [km]
104
(b) projection onto x-y plane
0 5 10 15
x [km] 104
-8 -6 -4 -2 0 2 4 6 8
z [km]
105
(c) projection onto x-z plane
-1 -0.5 0 0.5 1
y [km] 105
-8 -6 -4 -2 0 2 4 6 8
z [km]
105
(d) projection onto y-z plane
Figure 3.9: Periodic orbit obtained by differential correction step (zc2-fixed case). The red and blue orbits are obtained by the center manifold step and differential correction step, respectively. The initial positions are shown by asterisks.
is randomly chosen as
ξ= [
0 −5 −11 0 ]T
×10−3 (3.40)
This parameter vector generates a quasi-halo orbit, which is shown by blue lines in Figs. 3.10 and 3.11. By setting the Poincaré section to the x-y plane, the obtained quasi-periodic orbit is characterized by the returning time Tr=3.1036 and the winding time Tw=34.378. Therefore, the initial guess of the half-period of the periodic orbit is set as half of the returning time:
tc=1.5518.
First, zc2 is fixed and zc3 is modified in the differential correction step, as is tc. The initial set(ξ,tc)is corrected by Eq. (3.30), and it finally converges to
ξ= [
0 −5 −12.688 0 ]T
×10−3, tc=1.5470 (3.41) This corrected set(ξ,tc)forms the exact periodic orbit, which is the halo orbit shown by the red line in Fig. 3.10. Note that the z-component after the differential correction step is the same as that of the quasi-halo orbit.
Similarly, the set(ξ,tc)can be corrected by Eq. (3.31) after the iteration of convergence to ξ=
[
0 −1.8058 −11 0 ]T
×10−3, tc=1.5511 (3.42) This corrected set(ξ,tc)generates the halo orbit shown by the red line in Fig. 3.11. Note that in contrast the zc2-fixed case, all the components of the initial position are slightly changed by the differential correction step. The reason is that fixing zc3 does not imply fixing any component of the position vector but fixing the ratio of x and ˙y from Eq. (3.11).
Whereas the classical differential correction method cannot provide a quasi-periodic orbit, quasi-periodic orbits can be found by a simple procedure, which can be modified to obtain periodic orbits using the proposed method.