- Email: [email protected]
A new method of patched-conic for interplanetary orbit
Optik 156 (2018) 121–127
Contents lists available at ScienceDirect
Optik journal homepage: www.elsevier.de/ijleo
Original research article
A new method of patched-conic for interplanetary orbit Jin Li ∗ , Jianhui Zhao, Fan Li School of Instrumentation Science and Opto-Electronics Engineering, Beihang University, Beijing 100191, China
a r t i c l e
i n f o
Article history: Received 8 August 2017 Accepted 23 October 2017 Keywords: Patched-conic Sphere of influence Orbit constraints Continuous Iteration
a b s t r a c t In this paper, a new method of patched-conic for designing the interplanetary orbit is proposed. It is based on iterative computation and the equations are derived. The new method avoids the errors caused by assuming the velocity vector at the boundary of planetary sphere of influence parallel to the asymptote of the hyperbolic orbit. Based on the orbit constraint, the spacecraft’s velocity vector and the position vector at the boundary of planetary sphere of influence are calculated during each iteration to meet the continuous condition. This method is suitable for implementation in computer program and has a fast convergence speed. © 2017 Elsevier GmbH. All rights reserved.
1. Introduction The method of patched-conic plays a fundamental role in orbit design of space exploration [1–4]. It provides an appropriately approximation to precision orbit by sequence of two-body orbits. Taking the orbit from earth to mars as an example, the departure orbit is a hyperbolic orbit relative to earth; the transit orbit is an elliptical orbit relative to sun; the arrival orbit is a hyperbolic orbit relative to mars. For each of the three orbits is within the sphere of influence of the focus celestial body, the assumption is that only one gravitational center is considered. Each orbit can be calculated by classical two-body orbit functions. To connect three orbits, the differences between the velocity vectors of the two orbits at boundary should be within the ordered threshold. This is the continuous condition. There are some orbit constraints such as inclination, periapsis radius and so on. In order to meet the above two aspects, the method of patched-conic should be used. Many scholars have studied this method for years. In the early 1960 s, J.V. Breakwell, L.M. Perko and V.R. Bond developed the well-known method of patched-conic for designing a single-leg interplanetary matched-conic trajectory, round-trip interplanetary and matched-conic trajectory [5,6]. In the recent work of Y.H. Zhao and S.X. Yu [7,8], we can find the detail procedures for the above patched-conic method. R.H. Battin proposed the method of the design of the interplanetary freereturn, flyby orbits [9]. Some scholars gave the method to compute the hyperbolic excess speed and the angle of the perigee of the departure hyperbola relative to the earth’s orbital velocity vector [10,11], this method can be used in patched-conic orbit. For more advanced orbit design, N Bradley, R.P. Russell developed a method that converts space trajectories from the patched-conic models to the full gravity models [12]. In the design of the transit orbit between the earth and the moon, B.B.K. Reddy et al. used the micro genetic algorithms to determine the initial condition of the patched-conic approximation that achieve the minimum periselenium distance [13]. Based on iterative computation, a new method of patched-conic for interplanetary orbit design is proposed. Compared with traditional methods, the errors caused by the computation of position and velocity at the boundary of sphere of influence
∗ Corresponding author. E-mail address: [email protected] (J. Li). https://doi.org/10.1016/j.ijleo.2017.10.153 0030-4026/© 2017 Elsevier GmbH. All rights reserved.
122
J. Li et al. / Optik 156 (2018) 121–127
Fig. 1. The difference between the new method and the traditional methods.
can be avoided. It is less dependence on the geometric properties and suitable for computer programming. It has the fast convergence speed. The detailed steps are described and can be directly used in orbit computation. This paper is organized as below. The background and foundation are introduced in Section 2. The procedure of the proposed method is presented in Section 3. The simulation of an orbit from earth to mars is designed in Section 4. Finally, the conclusion is in Section 5. 2. Background and foundation 2.1. Improvement in the new method In traditional methods, both in departure orbit and arrival orbit, the velocity vector of the spacecraft at boundary of sphere of influence is assumed parallel to the asymptote of the hyperbolic orbit (as the Vg’ see in Fig. 1). However, the hyperbola infinitely approaches its asymptote but never cross. So the velocity vector along the tangent of the orbit (as the Vg see in Fig. 1) is not strictly parallel to the asymptote. To prove it, the function of the hyperbola is defined as x2 y2 − 2 =1 2 a b
(1)
So the function of the asymptote is b y = ± x a
(2)
The function of the tangent is dy b x = ± a dx 2 x − a2
(3)
When x is large enough, we could consider that (3) approximates to (2). However, it would bring small error. Because the mission time of the interplanetary explore is long, small error could be accumulated and might cause divergence in precise orbit design. Fig. 1 shows the situation of the spacecraft in the departure orbit, it is almost the same situation in the arrival orbit besides the velocity orientation is reverse. The proposed method less depends on the geometric properties. The Vg and rg (See in Fig. 1, both of the departure planet and the arrival planet) can be calculated by iterative computation to satisfy the orbit constraints and the continuous condition. 2.2. Lambert algorithm for transit orbit In the interplanetary orbit, the transit orbit takes up most of the mission time. Lambert algorithm is used to determine this orbit with two given position vectors and transfer time. The departure position is at the boundary of sphere of influence of departure planet. The arrival position is at the boundary of sphere of influence of target planet. The mission time could be
J. Li et al. / Optik 156 (2018) 121–127
123
Fig. 2. The block diagram of the system design principle.
used as the first approximation of the transfer time. The precise determination of these coefficients is the main part of this proposed method. The Lambert algorithm is the fundamental algorithm of the Keplerian orbit. It has been widely studied and used in orbit design [14–16]. R.H. Battin and R. Vaughan developed an elegant Lambert algorithm based on summarized and improved the previous algorithm [17]. It has the fast convergence characteristic and removes the singularity. For universal property and maturity, it has been chosen in this paper. 2.3. Orbit constraints In this paper, the parking orbit and the hyperbolic orbit are coplanar. The transfer between the coplanar orbits is by one tangent impulse at the injection point. The impulse direction is aligned with the velocity vector of the spacecraft [18]. It is the same for both the departure orbit and arrival orbit. In the case of multiple impulse transfer between the non-coplanar orbits, the orbit before the last impulse is also coplanar to the hyperbolic orbit [19]. The parking orbit is defined as a circular orbit, so the injection point is the periapsis of the hyperbolic orbit. Because of the detection requirement, the restriction of measurement network, the launch site and the tangent impulse, there are some constraints on the departure orbit and the arrival orbit. The work of Q.Q. Luo et al. about the free-return lunar flyby trajectory has been referenced [20]. The orbit constraints include: the flight-direction angle is 90 deg at the injection point, the inclination of the orbit is i and the modulus of periapsis radius rp (See in Fig. 1) is Rp . The departure orbit and arrival orbit have respective three constraints. 3. Procedure of the proposed patched method This method is composed of orbit constraints, orbit initialization and iterative process (See in Fig. 2). They are introduced below. The constraints are introduced in Section 2.3. 3.1. Initialization The initialization contains two parts as below. 1) The departure time td (injection point on the parking orbit relative to the departure planet), the arrival time ta (injection point on the parking orbit relative to the target planet), and the transfer direction are determined by mission requirement. The coordinate system is the same as used in ephemeris. At the boundary of sphere of influence, the velocity vector of the spacecraft relative to the planet is defined as Vg = [Vx ; Vy ; Vz ], the position vector relative to the planet is defined as rg = [rx ; ry ; rz ] (See in Fig. 1). The modulus of the rg is defined as Rg which is the radius of the sphere of influence of the planet. 2) The departure time td of the interplanetary orbit is used as the time tds when the spacecraft is at the boundary of sphere of influence of departure planet. The position vector of the departure planet relative to the sun at tds is obtained from ephemeris and used as the position vector of the spacecraft relative to the sun at the same time. The arrival time ta is used as the time tas when the spacecraft is at the boundary of the sphere of influence of the target planet. The position vector of the arrival planet relative to the sun at tas is obtained from ephemeris and used as the position vector of the spacecraft relative to the sun. By above statement, the transit orbit could be determined by solving the Lambert’s algorithm. So the velocity vectors of the spacecraft relative to sun at tds and tas can be calculated. The velocity vector of the departure planet relative to sun at tds can be obtained from ephemeris. This planet velocity vector is subtracted from the spacecraft velocity vector at tds relative to the sun, and then the spacecraft velocity vector Vg relative to the planet at the boundary
124
J. Li et al. / Optik 156 (2018) 121–127
of sphere of influence of the departure planet is obtained. The Vg about the arrival planet can be calculated in the same way. These two Vg are used in first iteration. After initialization is completed, the iterative process starts. 3.2. Iterative process The velocity vector and position vector respectively at the boundary of each planetary sphere of influence are calculated in each iteration until they all satisfy the continuous condition. The tds and tas are also corrected in each iteration. The motion of the spacecraft in the departure orbit and arrival orbit are considered reversely. In iteration, the process of computation about the two orbits are almost the same, so they would not be described respectively and only the differences would be pointed out. However, in each iteration, both the computation about the two orbits should be executed respectively. The iterative process steps are introduced as follows. Step 1: The current iteration is started with the two Vg which are calculated from the last iteration. If this is the first iteration, the two Vg from initialization are used. Step 2: For conversation of energy, the excess speed of the spacecraft relatively to the planet which is at the focus of the hyperbolic orbit is
Vg 2 −
V∞ =
2 Rg
(4)
is the gravity constant of the planet. Semi-major axis of the hyperbolic orbit is a=−
(5)
V∞ 2
The velocity vector modulus of the spacecraft at periapsis of the hyperbolic orbit is
VP =
1 2 − Rp a
(6)
The eccentricity of the hyperbolic orbit is e=1−
Rp a
(7)
Step 3: The angular momentum vector of the spacecraft relative to the planet can be defined as h = [hx ; hy ; hz ]. For the angular momentum vector is constant, the modulus of angular momentum is defined as ||h|| = Rp ||Vp || when the spacecraft is at the periapsis of the hyperbolic orbit. For the inclination is i, hz = ||h||cosi, we have 2
2
h2x + h2y = h (sin i)
(8)
Because of rg × Vg = h, Vg ⊥h can be expressed as Vg h = Vx hx + Vy hy + Vz hz = 0. So, the hx is hx = −
Vy hy + Vz hz Vy hy + Vz h cos i =− Vx Vx
(9)
Substituting Eq. (9) into Eq. (8), we have 2
2
(Vx2 + Vy2 )h2y + 2Vy Vz hy h cos i + (Vz h cos i) − (Vx h sin i) = 0 From Eq. (10), the hy is
hy =
−Vy Vz h cos i ±
2
Vx2 ((Vx h sin i) + Vy h sin i
2
(10)
2
− (Vz h cos i) )
(11)
(Vx2 + Vy2 )
Substituting Eq. (11) into Eq. (9), the hx is
hx = −
Vx2 Vz h cos i ± Vy
2
Vx2 ((Vx h sin i) + Vy h sin i
2
2
− (Vz h cos i) )
(12)
Vx (Vx2 + Vy2 )
The ‘±’ in hx and hy are both ‘+’ or ‘−’ at the same time. So h is obtained. By the combination of hx and hy in the departure and arrival orbit, four orbits can be obtained. The orbit is chosen according to the mission demand and then rg can be calculated as follows. Step 4: When the spacecraft is at the boundary of sphere of influence, the angle between rg and Vg is defined as . Because of ||h|| = ||rg ||.||Vp ||sin , the is In departure orbit
= sin−1
h rg · Vg
0<≤
In arrival orbit 2
= − sin−1
h rg · Vg
≤< 2
(13)
J. Li et al. / Optik 156 (2018) 121–127
125
Step 5: By Eq. (13), the relation between rg and Vg is rg Vg = rx Vx + ry Vy + rz Vz = rg · Vg cos =
(14)
The detailed description of rg × Vg = h is ry Vz − rz Vy = hx
(15)
rz Vx − rx Vz = hy
(16)
rx Vy − ry Vx = hz
(17)
From Eq. (16), we have rz = (rx Vz + hy )/Vx . From Eq. (17), we have ry = (rx Vy − hz )/Vx . Substituting ry and rz into Eq. (14), the rx is rx =
Vx + hz Vy − hy Vz Vx2 + Vy2 + Vz2
(18)
From Eq. (15), we have rz = (ry Vz − hx )/Vy . From Eq. (17), we have rx = (ry Vx + hz )/Vy . Substituting rx and rz into Eq. (14), the ry is ry =
Vy + hx Vz − hz Vx Vx2 + Vy2 + Vz2
(19)
From Eq. (16), we have rx = (rz Vx − hy )/Vz . From Eq. (15), we have ry = (rz Vy + hx )/Vz . Substituting rx and ry into Eq. (14), the rz is rz =
Vz + hy Vx − hx Vy Vx2 + Vy2 + Vz2
(20)
From above, the rg of the spacecraft at the each boundary of planetary sphere of influence is obtained. In implementation, the above computation is executed respectively for the departure orbit and arrival orbit in each iteration. Step 6: The hyperbolic eccentric anomaly of the spacecraft at the boundary of planetary sphere of influence is Rg 1 )) H = cosh−1 ( (1 − e a
(21)
In departure orbit, from periapsis to the planetary boundary of sphere of influence, the time is
t0 =
−a2 (e sinh H − H)
(22)
In the arrival orbit, from the boundary of sphere of influence to periapsis, the time is
t1 =
−a2 (e sinh H − H)
(23)
By Eq. (22), the tds is corrected by tds = td + t0 . By Eq. (23), the tas is corrected by tas = ta − t1 . They would be used in the next iteration. Step 7: The position vector of the departure planet relative to the sun at tds and the position vector of the arrival planet relative to the sun at tas could be obtained from the ephemeris. By the rg of the two orbits, the position vector of the spacecraft relative to the sun at tds and tas can be relatively calculated. For the transit orbit, the spacecraft position vector at tds is used as the departure position vector, the spacecraft position vector at tas is used as the arrival position vector. So the transit orbit can be calculated by Lambert algorithm, the velocity vectors of the spacecraft relative to sun at tds and tas can be relatively obtained. Then the velocity vectors of the spacecraft relative to the planet at tds and tas can be relatively calculated. These two velocity vectors are compared with the two Vg in the current iteration. If the differences between them converge to smaller than the ordered threshold, the continuous condition is met and the iterative process is ended. Otherwise the next iteration starts. These two velocity vectors are used as the new Vg in the next iteration. The differences of velocity vectors are stored. If the differences are divergent or become constant but larger than the threshold, these three orbits cannot be patched. The td , ta and orbit constraints need to be adjusted. In this method, the minimum ordered threshold is determined by the accuracy of the Lambert algorithm and it is also the accuracy of this method of patched-conic. From steps 1) to 7), they form the iterative process. By the patched orbit, the rg and Vg of the departure orbit are used to solve the Kepler’s equation. Then the velocity vector and position vector of the spacecraft at the periapsis of departure orbit can be calculated and used as the initial value in precise design of the interplanetary orbit. 4. Numerical simulation The simulation of the proposed method is the design of an interplanetary orbit from Earth to Mars. By energy optimization based on Hohmann transfer orbit, the departure time (injection time) is 00:00:00.1 on January 12, 2014 and the arrival time
126
J. Li et al. / Optik 156 (2018) 121–127
Fig. 3. The departure orbit with detail description near the earth.
Fig. 4. The arrival orbit with detail description near the mars.
Fig. 5. The whole patched orbit in the example.
is 00:00:00.1 on 21 August 2014. The constraints of the departure orbit are: periapsis radius is 300 km, inclination of the orbit is 45 deg and the flight direction angle is 90 deg at the injection point. The constraints of the arrival orbit are: periapsis radius is 800 km, inclination of the orbit is 80 deg and the flight direction angle is 90 deg at the injection point. The threshold for the velocity is defined as 10−10 km/s. The simulation is developed in Matlab with Intel CPU i5 M450, 4GB RAM and Windows7 32bit operating system. The patched-conic orbit is obtained by 16 iterations and the total running time is 7.784263s. The departure orbit (See in Fig. 3) and arrival orbit (See in Fig. 4) can be calculated by Kepler’s equation. Then the whole interplanetary orbit (See in Fig. 5) from Earth to Mars is obtained. From the curves in Fig. 3 and Fig. 4, all orbits are also connected smoothly. From the curves in Fig. 6, the initial differences between the velocities of the two orbits at the injection point tend to convergence within five iterations. The proposed method in this paper can be effective in patched-conic for interplanetary orbit. 5. Conclusions A new method of patched-conic is proposed in this paper and the whole steps are introduced. This method avoids the errors caused by using the asymptote direction as the velocity direction at the boundary of planetary sphere of influence,
J. Li et al. / Optik 156 (2018) 121–127
127
Fig. 6. Differences of the velocities calculation in the iterative process (a) departure orbit (b) arrival orbit.
so the precision of the patched orbit is improved. Moreover, it is suitable for implementation by computer programming by avoiding too much dependence on the geometric properties of conic sections. The simulation of an orbit from earth to mars design is presented to demonstrate this method and it also shows that it is fast on convergence. It has been used to test the navigation algorithm for deep space exploring by the design of test orbit. It also has the potential for testing the algorithm and principle of the optical sensor used in space explorer such as star sensor and infrared earth sensor [21,22], by supply the position, velocity and attitude of the spacecraft and relative position of the stars. Acknowledgment This work was supported by National Natural Science Foundation of China (Grant no. 61233005 and 51574012). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]
J.W. Crenshaw, Sphere of influence in patched-conic methods, AIAA J. 1 (1963) 2168–2170. G.A. Henning, P.J. Edelman, J.M. Longuski, Design and optimization of interplanetary aerogravity-assist tours, J. Spacecr. Rocket 51 (2014) 1849–1856. S. Campagnola, Y. Kawakatsu, Jupiter magnetospheric orbiter: trajectory design in the jovian system, J. Spacecr. Rocket 49 (2012) 318–324. R.P. Russell, N.J. Strange, Cycler trajectories in planetary moon systems, J. Guid. Control Dyn. 32 (2009) 143–157. J.V. Breakwell, L.M. Perko, Matched asymptotic expansions, patched conies, and the computation of interplanetary trajectories, in: AIAA/ION Astrodynamics Specialist Conference, Sep 16–17, Monterey, California, 1965. V.R. Bond, Matched-conic solutions to round-trip interplanetary trajectory problems that insure state-vector continuity at all boundaries, in: NASA Technical Note, 1969 (Report NO. D-4942). Y.H. Zhao, Trajectory Design and Orbital Dynamics of Deep Space Exploration, Nanjing University, 2012. S.X. Yu, Orbit Analysis, Design and Control in Deep Space, Nanjing University, 2013. R.H. Battin, Rev ed., in: An Introduction to the Mathematics and Methods of Astrodynamics, AIAA Inc, Reston, Virginia, 1999, pp. 419–446. H.D. Curtis, Orbital Mechanics for Engineering Students, Elsevier Butterworth-Heinemann, Oxford, 2005, pp. 359–366. H. Schaub, J.L. Junkins, Analytical Mechanics of Space Systems, AIAA Inc, Reston, Virginia, 2002, pp. 455–472. N. Bradley, R.P. Russell, A continuation method for converting trajectories from patched conics to full gravity models, J. Astronaut. Sci. 61 (2014) 227–254. B.B.K. Reddy, B. Kimiaghalam, A. Homaifar, Evolutionary algorithms for parameter determination of patched conic approximation, in: Aerospace Conference, 2005 March 5–12, Big Sky, MT. IEEE, 2005, pp. 501–506. R. Allen, R. Wenzel, Tan root for Battin’s lambert recipe, J. Guid. Control Dyn. 38 (2015) 1826–1827. S.L. Nelson, P. Zarchan, Alternative approach to the solution of Lamber’s problem, J. Guid. Control Dyn. 15 (1992) 1003–1009. J. Ahn, S.-I.I. Lee, Lambert algorithm using analytic gradients, J. Guid. Control Dyn. 36 (2013) 1751–1761. R.H. Battin, R.M. Vaughan, An elegant lambert algorithm, J. Guid. Control Dyn. 7 (1984) 662–670. G. Zhang, X.Y. Zhang, X.B. Cao, Tangent-impulse transfer from elliptic orbit to an excess velocity vector, Chin. J. Aeronaut. 27 (2014) 577–583. D.R. Jones, C. Ocampo, Optimization of impulsive trajectories from circular orbit to an excess velocity vector, J. Guid. Control Dyn. 35 (2012) 234–244. Q.Q. Luo, J.F. Yin, C. Han, Design of earth—moon free-return trajectories, J. Guid. Control Dyn. 36 (2013), 271–263. X.B. Hu, J.H. Zhao, Y. Zhao, Z.J. Tu, Moon/Sun interference analysis, identification and suppression on Dual Cone Earth Sensor, Optik 127 (2016) 1665–1670. X.B. Hu, J.H. Zhao, Z.J. Tu, Research on the method of suppressing sun and moon’s interference on infrared conical earth sensor, Proc. Inst. Mech. Eng. C 229 (2015) 399–406.
Contents lists available at ScienceDirect
Optik journal homepage: www.elsevier.de/ijleo
Original research article
A new method of patched-conic for interplanetary orbit Jin Li ∗ , Jianhui Zhao, Fan Li School of Instrumentation Science and Opto-Electronics Engineering, Beihang University, Beijing 100191, China
a r t i c l e
i n f o
Article history: Received 8 August 2017 Accepted 23 October 2017 Keywords: Patched-conic Sphere of influence Orbit constraints Continuous Iteration
a b s t r a c t In this paper, a new method of patched-conic for designing the interplanetary orbit is proposed. It is based on iterative computation and the equations are derived. The new method avoids the errors caused by assuming the velocity vector at the boundary of planetary sphere of influence parallel to the asymptote of the hyperbolic orbit. Based on the orbit constraint, the spacecraft’s velocity vector and the position vector at the boundary of planetary sphere of influence are calculated during each iteration to meet the continuous condition. This method is suitable for implementation in computer program and has a fast convergence speed. © 2017 Elsevier GmbH. All rights reserved.
1. Introduction The method of patched-conic plays a fundamental role in orbit design of space exploration [1–4]. It provides an appropriately approximation to precision orbit by sequence of two-body orbits. Taking the orbit from earth to mars as an example, the departure orbit is a hyperbolic orbit relative to earth; the transit orbit is an elliptical orbit relative to sun; the arrival orbit is a hyperbolic orbit relative to mars. For each of the three orbits is within the sphere of influence of the focus celestial body, the assumption is that only one gravitational center is considered. Each orbit can be calculated by classical two-body orbit functions. To connect three orbits, the differences between the velocity vectors of the two orbits at boundary should be within the ordered threshold. This is the continuous condition. There are some orbit constraints such as inclination, periapsis radius and so on. In order to meet the above two aspects, the method of patched-conic should be used. Many scholars have studied this method for years. In the early 1960 s, J.V. Breakwell, L.M. Perko and V.R. Bond developed the well-known method of patched-conic for designing a single-leg interplanetary matched-conic trajectory, round-trip interplanetary and matched-conic trajectory [5,6]. In the recent work of Y.H. Zhao and S.X. Yu [7,8], we can find the detail procedures for the above patched-conic method. R.H. Battin proposed the method of the design of the interplanetary freereturn, flyby orbits [9]. Some scholars gave the method to compute the hyperbolic excess speed and the angle of the perigee of the departure hyperbola relative to the earth’s orbital velocity vector [10,11], this method can be used in patched-conic orbit. For more advanced orbit design, N Bradley, R.P. Russell developed a method that converts space trajectories from the patched-conic models to the full gravity models [12]. In the design of the transit orbit between the earth and the moon, B.B.K. Reddy et al. used the micro genetic algorithms to determine the initial condition of the patched-conic approximation that achieve the minimum periselenium distance [13]. Based on iterative computation, a new method of patched-conic for interplanetary orbit design is proposed. Compared with traditional methods, the errors caused by the computation of position and velocity at the boundary of sphere of influence
∗ Corresponding author. E-mail address: [email protected] (J. Li). https://doi.org/10.1016/j.ijleo.2017.10.153 0030-4026/© 2017 Elsevier GmbH. All rights reserved.
122
J. Li et al. / Optik 156 (2018) 121–127
Fig. 1. The difference between the new method and the traditional methods.
can be avoided. It is less dependence on the geometric properties and suitable for computer programming. It has the fast convergence speed. The detailed steps are described and can be directly used in orbit computation. This paper is organized as below. The background and foundation are introduced in Section 2. The procedure of the proposed method is presented in Section 3. The simulation of an orbit from earth to mars is designed in Section 4. Finally, the conclusion is in Section 5. 2. Background and foundation 2.1. Improvement in the new method In traditional methods, both in departure orbit and arrival orbit, the velocity vector of the spacecraft at boundary of sphere of influence is assumed parallel to the asymptote of the hyperbolic orbit (as the Vg’ see in Fig. 1). However, the hyperbola infinitely approaches its asymptote but never cross. So the velocity vector along the tangent of the orbit (as the Vg see in Fig. 1) is not strictly parallel to the asymptote. To prove it, the function of the hyperbola is defined as x2 y2 − 2 =1 2 a b
(1)
So the function of the asymptote is b y = ± x a
(2)
The function of the tangent is dy b x = ± a dx 2 x − a2
(3)
When x is large enough, we could consider that (3) approximates to (2). However, it would bring small error. Because the mission time of the interplanetary explore is long, small error could be accumulated and might cause divergence in precise orbit design. Fig. 1 shows the situation of the spacecraft in the departure orbit, it is almost the same situation in the arrival orbit besides the velocity orientation is reverse. The proposed method less depends on the geometric properties. The Vg and rg (See in Fig. 1, both of the departure planet and the arrival planet) can be calculated by iterative computation to satisfy the orbit constraints and the continuous condition. 2.2. Lambert algorithm for transit orbit In the interplanetary orbit, the transit orbit takes up most of the mission time. Lambert algorithm is used to determine this orbit with two given position vectors and transfer time. The departure position is at the boundary of sphere of influence of departure planet. The arrival position is at the boundary of sphere of influence of target planet. The mission time could be
J. Li et al. / Optik 156 (2018) 121–127
123
Fig. 2. The block diagram of the system design principle.
used as the first approximation of the transfer time. The precise determination of these coefficients is the main part of this proposed method. The Lambert algorithm is the fundamental algorithm of the Keplerian orbit. It has been widely studied and used in orbit design [14–16]. R.H. Battin and R. Vaughan developed an elegant Lambert algorithm based on summarized and improved the previous algorithm [17]. It has the fast convergence characteristic and removes the singularity. For universal property and maturity, it has been chosen in this paper. 2.3. Orbit constraints In this paper, the parking orbit and the hyperbolic orbit are coplanar. The transfer between the coplanar orbits is by one tangent impulse at the injection point. The impulse direction is aligned with the velocity vector of the spacecraft [18]. It is the same for both the departure orbit and arrival orbit. In the case of multiple impulse transfer between the non-coplanar orbits, the orbit before the last impulse is also coplanar to the hyperbolic orbit [19]. The parking orbit is defined as a circular orbit, so the injection point is the periapsis of the hyperbolic orbit. Because of the detection requirement, the restriction of measurement network, the launch site and the tangent impulse, there are some constraints on the departure orbit and the arrival orbit. The work of Q.Q. Luo et al. about the free-return lunar flyby trajectory has been referenced [20]. The orbit constraints include: the flight-direction angle is 90 deg at the injection point, the inclination of the orbit is i and the modulus of periapsis radius rp (See in Fig. 1) is Rp . The departure orbit and arrival orbit have respective three constraints. 3. Procedure of the proposed patched method This method is composed of orbit constraints, orbit initialization and iterative process (See in Fig. 2). They are introduced below. The constraints are introduced in Section 2.3. 3.1. Initialization The initialization contains two parts as below. 1) The departure time td (injection point on the parking orbit relative to the departure planet), the arrival time ta (injection point on the parking orbit relative to the target planet), and the transfer direction are determined by mission requirement. The coordinate system is the same as used in ephemeris. At the boundary of sphere of influence, the velocity vector of the spacecraft relative to the planet is defined as Vg = [Vx ; Vy ; Vz ], the position vector relative to the planet is defined as rg = [rx ; ry ; rz ] (See in Fig. 1). The modulus of the rg is defined as Rg which is the radius of the sphere of influence of the planet. 2) The departure time td of the interplanetary orbit is used as the time tds when the spacecraft is at the boundary of sphere of influence of departure planet. The position vector of the departure planet relative to the sun at tds is obtained from ephemeris and used as the position vector of the spacecraft relative to the sun at the same time. The arrival time ta is used as the time tas when the spacecraft is at the boundary of the sphere of influence of the target planet. The position vector of the arrival planet relative to the sun at tas is obtained from ephemeris and used as the position vector of the spacecraft relative to the sun. By above statement, the transit orbit could be determined by solving the Lambert’s algorithm. So the velocity vectors of the spacecraft relative to sun at tds and tas can be calculated. The velocity vector of the departure planet relative to sun at tds can be obtained from ephemeris. This planet velocity vector is subtracted from the spacecraft velocity vector at tds relative to the sun, and then the spacecraft velocity vector Vg relative to the planet at the boundary
124
J. Li et al. / Optik 156 (2018) 121–127
of sphere of influence of the departure planet is obtained. The Vg about the arrival planet can be calculated in the same way. These two Vg are used in first iteration. After initialization is completed, the iterative process starts. 3.2. Iterative process The velocity vector and position vector respectively at the boundary of each planetary sphere of influence are calculated in each iteration until they all satisfy the continuous condition. The tds and tas are also corrected in each iteration. The motion of the spacecraft in the departure orbit and arrival orbit are considered reversely. In iteration, the process of computation about the two orbits are almost the same, so they would not be described respectively and only the differences would be pointed out. However, in each iteration, both the computation about the two orbits should be executed respectively. The iterative process steps are introduced as follows. Step 1: The current iteration is started with the two Vg which are calculated from the last iteration. If this is the first iteration, the two Vg from initialization are used. Step 2: For conversation of energy, the excess speed of the spacecraft relatively to the planet which is at the focus of the hyperbolic orbit is
Vg 2 −
V∞ =
2 Rg
(4)
is the gravity constant of the planet. Semi-major axis of the hyperbolic orbit is a=−
(5)
V∞ 2
The velocity vector modulus of the spacecraft at periapsis of the hyperbolic orbit is
VP =
1 2 − Rp a
(6)
The eccentricity of the hyperbolic orbit is e=1−
Rp a
(7)
Step 3: The angular momentum vector of the spacecraft relative to the planet can be defined as h = [hx ; hy ; hz ]. For the angular momentum vector is constant, the modulus of angular momentum is defined as ||h|| = Rp ||Vp || when the spacecraft is at the periapsis of the hyperbolic orbit. For the inclination is i, hz = ||h||cosi, we have 2
2
h2x + h2y = h (sin i)
(8)
Because of rg × Vg = h, Vg ⊥h can be expressed as Vg h = Vx hx + Vy hy + Vz hz = 0. So, the hx is hx = −
Vy hy + Vz hz Vy hy + Vz h cos i =− Vx Vx
(9)
Substituting Eq. (9) into Eq. (8), we have 2
2
(Vx2 + Vy2 )h2y + 2Vy Vz hy h cos i + (Vz h cos i) − (Vx h sin i) = 0 From Eq. (10), the hy is
hy =
−Vy Vz h cos i ±
2
Vx2 ((Vx h sin i) + Vy h sin i
2
(10)
2
− (Vz h cos i) )
(11)
(Vx2 + Vy2 )
Substituting Eq. (11) into Eq. (9), the hx is
hx = −
Vx2 Vz h cos i ± Vy
2
Vx2 ((Vx h sin i) + Vy h sin i
2
2
− (Vz h cos i) )
(12)
Vx (Vx2 + Vy2 )
The ‘±’ in hx and hy are both ‘+’ or ‘−’ at the same time. So h is obtained. By the combination of hx and hy in the departure and arrival orbit, four orbits can be obtained. The orbit is chosen according to the mission demand and then rg can be calculated as follows. Step 4: When the spacecraft is at the boundary of sphere of influence, the angle between rg and Vg is defined as . Because of ||h|| = ||rg ||.||Vp ||sin , the is In departure orbit
= sin−1
h rg · Vg
0<≤
In arrival orbit 2
= − sin−1
h rg · Vg
≤< 2
(13)
J. Li et al. / Optik 156 (2018) 121–127
125
Step 5: By Eq. (13), the relation between rg and Vg is rg Vg = rx Vx + ry Vy + rz Vz = rg · Vg cos =
(14)
The detailed description of rg × Vg = h is ry Vz − rz Vy = hx
(15)
rz Vx − rx Vz = hy
(16)
rx Vy − ry Vx = hz
(17)
From Eq. (16), we have rz = (rx Vz + hy )/Vx . From Eq. (17), we have ry = (rx Vy − hz )/Vx . Substituting ry and rz into Eq. (14), the rx is rx =
Vx + hz Vy − hy Vz Vx2 + Vy2 + Vz2
(18)
From Eq. (15), we have rz = (ry Vz − hx )/Vy . From Eq. (17), we have rx = (ry Vx + hz )/Vy . Substituting rx and rz into Eq. (14), the ry is ry =
Vy + hx Vz − hz Vx Vx2 + Vy2 + Vz2
(19)
From Eq. (16), we have rx = (rz Vx − hy )/Vz . From Eq. (15), we have ry = (rz Vy + hx )/Vz . Substituting rx and ry into Eq. (14), the rz is rz =
Vz + hy Vx − hx Vy Vx2 + Vy2 + Vz2
(20)
From above, the rg of the spacecraft at the each boundary of planetary sphere of influence is obtained. In implementation, the above computation is executed respectively for the departure orbit and arrival orbit in each iteration. Step 6: The hyperbolic eccentric anomaly of the spacecraft at the boundary of planetary sphere of influence is Rg 1 )) H = cosh−1 ( (1 − e a
(21)
In departure orbit, from periapsis to the planetary boundary of sphere of influence, the time is
t0 =
−a2 (e sinh H − H)
(22)
In the arrival orbit, from the boundary of sphere of influence to periapsis, the time is
t1 =
−a2 (e sinh H − H)
(23)
By Eq. (22), the tds is corrected by tds = td + t0 . By Eq. (23), the tas is corrected by tas = ta − t1 . They would be used in the next iteration. Step 7: The position vector of the departure planet relative to the sun at tds and the position vector of the arrival planet relative to the sun at tas could be obtained from the ephemeris. By the rg of the two orbits, the position vector of the spacecraft relative to the sun at tds and tas can be relatively calculated. For the transit orbit, the spacecraft position vector at tds is used as the departure position vector, the spacecraft position vector at tas is used as the arrival position vector. So the transit orbit can be calculated by Lambert algorithm, the velocity vectors of the spacecraft relative to sun at tds and tas can be relatively obtained. Then the velocity vectors of the spacecraft relative to the planet at tds and tas can be relatively calculated. These two velocity vectors are compared with the two Vg in the current iteration. If the differences between them converge to smaller than the ordered threshold, the continuous condition is met and the iterative process is ended. Otherwise the next iteration starts. These two velocity vectors are used as the new Vg in the next iteration. The differences of velocity vectors are stored. If the differences are divergent or become constant but larger than the threshold, these three orbits cannot be patched. The td , ta and orbit constraints need to be adjusted. In this method, the minimum ordered threshold is determined by the accuracy of the Lambert algorithm and it is also the accuracy of this method of patched-conic. From steps 1) to 7), they form the iterative process. By the patched orbit, the rg and Vg of the departure orbit are used to solve the Kepler’s equation. Then the velocity vector and position vector of the spacecraft at the periapsis of departure orbit can be calculated and used as the initial value in precise design of the interplanetary orbit. 4. Numerical simulation The simulation of the proposed method is the design of an interplanetary orbit from Earth to Mars. By energy optimization based on Hohmann transfer orbit, the departure time (injection time) is 00:00:00.1 on January 12, 2014 and the arrival time
126
J. Li et al. / Optik 156 (2018) 121–127
Fig. 3. The departure orbit with detail description near the earth.
Fig. 4. The arrival orbit with detail description near the mars.
Fig. 5. The whole patched orbit in the example.
is 00:00:00.1 on 21 August 2014. The constraints of the departure orbit are: periapsis radius is 300 km, inclination of the orbit is 45 deg and the flight direction angle is 90 deg at the injection point. The constraints of the arrival orbit are: periapsis radius is 800 km, inclination of the orbit is 80 deg and the flight direction angle is 90 deg at the injection point. The threshold for the velocity is defined as 10−10 km/s. The simulation is developed in Matlab with Intel CPU i5 M450, 4GB RAM and Windows7 32bit operating system. The patched-conic orbit is obtained by 16 iterations and the total running time is 7.784263s. The departure orbit (See in Fig. 3) and arrival orbit (See in Fig. 4) can be calculated by Kepler’s equation. Then the whole interplanetary orbit (See in Fig. 5) from Earth to Mars is obtained. From the curves in Fig. 3 and Fig. 4, all orbits are also connected smoothly. From the curves in Fig. 6, the initial differences between the velocities of the two orbits at the injection point tend to convergence within five iterations. The proposed method in this paper can be effective in patched-conic for interplanetary orbit. 5. Conclusions A new method of patched-conic is proposed in this paper and the whole steps are introduced. This method avoids the errors caused by using the asymptote direction as the velocity direction at the boundary of planetary sphere of influence,
J. Li et al. / Optik 156 (2018) 121–127
127
Fig. 6. Differences of the velocities calculation in the iterative process (a) departure orbit (b) arrival orbit.
so the precision of the patched orbit is improved. Moreover, it is suitable for implementation by computer programming by avoiding too much dependence on the geometric properties of conic sections. The simulation of an orbit from earth to mars design is presented to demonstrate this method and it also shows that it is fast on convergence. It has been used to test the navigation algorithm for deep space exploring by the design of test orbit. It also has the potential for testing the algorithm and principle of the optical sensor used in space explorer such as star sensor and infrared earth sensor [21,22], by supply the position, velocity and attitude of the spacecraft and relative position of the stars. Acknowledgment This work was supported by National Natural Science Foundation of China (Grant no. 61233005 and 51574012). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]
J.W. Crenshaw, Sphere of influence in patched-conic methods, AIAA J. 1 (1963) 2168–2170. G.A. Henning, P.J. Edelman, J.M. Longuski, Design and optimization of interplanetary aerogravity-assist tours, J. Spacecr. Rocket 51 (2014) 1849–1856. S. Campagnola, Y. Kawakatsu, Jupiter magnetospheric orbiter: trajectory design in the jovian system, J. Spacecr. Rocket 49 (2012) 318–324. R.P. Russell, N.J. Strange, Cycler trajectories in planetary moon systems, J. Guid. Control Dyn. 32 (2009) 143–157. J.V. Breakwell, L.M. Perko, Matched asymptotic expansions, patched conies, and the computation of interplanetary trajectories, in: AIAA/ION Astrodynamics Specialist Conference, Sep 16–17, Monterey, California, 1965. V.R. Bond, Matched-conic solutions to round-trip interplanetary trajectory problems that insure state-vector continuity at all boundaries, in: NASA Technical Note, 1969 (Report NO. D-4942). Y.H. Zhao, Trajectory Design and Orbital Dynamics of Deep Space Exploration, Nanjing University, 2012. S.X. Yu, Orbit Analysis, Design and Control in Deep Space, Nanjing University, 2013. R.H. Battin, Rev ed., in: An Introduction to the Mathematics and Methods of Astrodynamics, AIAA Inc, Reston, Virginia, 1999, pp. 419–446. H.D. Curtis, Orbital Mechanics for Engineering Students, Elsevier Butterworth-Heinemann, Oxford, 2005, pp. 359–366. H. Schaub, J.L. Junkins, Analytical Mechanics of Space Systems, AIAA Inc, Reston, Virginia, 2002, pp. 455–472. N. Bradley, R.P. Russell, A continuation method for converting trajectories from patched conics to full gravity models, J. Astronaut. Sci. 61 (2014) 227–254. B.B.K. Reddy, B. Kimiaghalam, A. Homaifar, Evolutionary algorithms for parameter determination of patched conic approximation, in: Aerospace Conference, 2005 March 5–12, Big Sky, MT. IEEE, 2005, pp. 501–506. R. Allen, R. Wenzel, Tan root for Battin’s lambert recipe, J. Guid. Control Dyn. 38 (2015) 1826–1827. S.L. Nelson, P. Zarchan, Alternative approach to the solution of Lamber’s problem, J. Guid. Control Dyn. 15 (1992) 1003–1009. J. Ahn, S.-I.I. Lee, Lambert algorithm using analytic gradients, J. Guid. Control Dyn. 36 (2013) 1751–1761. R.H. Battin, R.M. Vaughan, An elegant lambert algorithm, J. Guid. Control Dyn. 7 (1984) 662–670. G. Zhang, X.Y. Zhang, X.B. Cao, Tangent-impulse transfer from elliptic orbit to an excess velocity vector, Chin. J. Aeronaut. 27 (2014) 577–583. D.R. Jones, C. Ocampo, Optimization of impulsive trajectories from circular orbit to an excess velocity vector, J. Guid. Control Dyn. 35 (2012) 234–244. Q.Q. Luo, J.F. Yin, C. Han, Design of earth—moon free-return trajectories, J. Guid. Control Dyn. 36 (2013), 271–263. X.B. Hu, J.H. Zhao, Y. Zhao, Z.J. Tu, Moon/Sun interference analysis, identification and suppression on Dual Cone Earth Sensor, Optik 127 (2016) 1665–1670. X.B. Hu, J.H. Zhao, Z.J. Tu, Research on the method of suppressing sun and moon’s interference on infrared conical earth sensor, Proc. Inst. Mech. Eng. C 229 (2015) 399–406.