Received 20 July 2016; accepted 11 September 2016; published 14 September 2016
The subject of halo orbits in restricted three-body problem (RTBP) has received considerable attention in the last four decades. A halo orbit is a periodic, three-dimensional orbitnear the collinear Lagrangian points in the three-body problem. Analytical solutions to generate the halo orbits around any of the collinear points can be obtained. Farquhar first used the name “halo” in his doctoral thesis  and proposed an idea for placing a satellite in an orbit around L2 of the Earth-Moon RTBP. Had this idea been implemented, one could continuously view the Earth and the dark side of the moon at the same time. A communication link satellite is not necessary, if satellites are placed around Lagrangian points of the respective systems. Breakwell & Brown  generated halo orbits for the Earth-Moon system. Richardson  contributed significantly to obtain the halo orbits, by finding an analytical solution using Lindstedt-Poincaré method. A good amount of research has been carried out till date in this interesting area of halo orbits. Some of the important contributions are by Howell  , Howell and Pernicka  , Folta and Richon  and Howell et al.  .
The classical model of the restricted three-body problem are relatively less accurate for studying the motion of any Sun-planet system as they do not account for the effect of the perturbing forces such as oblateness of planets, cosmic rays, magnetic field, radiation pressure etc. Cosmic rays are immensely high-energy radiation, mainly originating outside the solar system. They produce showers of secondary particles that penetrate and impact the Earth’s atmosphere and sometimes even reach the surface. They interact with gaseous and other matter at high altitudes and produce secondary radiation. The combination of both contributes to the space radiation environment. In addition to the particles originating from the Sun are particles from other stars and heavy ion sources such as nova and supernova in our galaxy and beyond. In interplanetary space these ionizing particles constitute the major radiation threat. These particles are influenced by planetary or Earth’s magnetic field to form radiation belts, which in Earth’s case are known as Van Allen Radiation belts, containing trapped electrons in the outer belt and protons in the inner belt. The composition and intensity of the radiation varies significantly with the trajectory of a space vehicle. Anomalies in communication satellite operation have been caused by the unexpected triggering of digital circuits by the cosmic rays. In our present study, we restrict ourselves with the solar radiation and oblateness of the planet.
Radzievskii  was the first one to study the effect of solar radiation pressure. He found out that the maximum force due to the radiation pressure acts in the radial direction, given by
where q is defined in terms of particle radius a, density δ and radiation pressure efficiency factor x as
is a variable, depending upon the nature of the third body (satellite). The value of q can be considered as a constant, if the fluctuations in the beam of solar radiation and the effect of planet’s shadow are neglected. Using the model of  , Dutt and Sharma  studied periodic orbits in the Sun-Mars system using the numerical technique of Poincare surface of sections and found out more than 74 periodic orbits.
Sharma and Subba Rao  introduced the oblateness of the more massive primary in the three-dimensional restricted three-body problem. It has two terms, one with a z term in the numerator. Sharma  studied the periodic orbits around the Lagrangian points in the planar RTBP by considering Sun as source of radiation and smaller primary as an oblate spheroid with its equatorial plane coincident with the plane of motion. Tiwary and Kushvah  followed the model of Sharma  to study the halo orbits around the Lagrangian points L1 and L2 analytically. However, they did not consider the z term in oblateness in their study. In the present work, we have considered both the terms due to oblateness of the smaller primary in the photogravitational restricted three-body problem to study the halo orbits analytically as well as numerically around L1.
The first satellite placed in the halo orbit at Sun-Earth L1 point was International Sun-Earth Explorer-3 (ISEE-3), launched in 1978. Solar Heliospheric Observatory (SOHO), launched in 1975 by NASA succeeded ISEE-3.We have taken data of the path of the SOHO mission from the mission website over a period of January-June 2008, to validate our analytical and numerical solutions.
2. Circular Restricted Three-Body Problem
The circular RTBP consists of two primary masses revolving in circular orbits around their centre of mass under the influence of their mutual gravity. The third body of infinitesimal mass moves under the gravitational effect of these two primaries (Figure 1). RTBP has five equilibrium points, called Lagrangian or libration points. These points are the points of zero velocities and an object placed in these points remains there. Out of the five Lagrangian points, three are collinear (L1, L2, L3) and the other two points (L4, L5) form equilateral triangles with the primaries. Although the Lagrangian point is just a point in empty space, its peculiar characteristic is that
Figure 1. Three-dimensional restricted three-body problem.
it can be orbited. Halo orbits are three-dimensional orbits around the collinear points.
The equations of motion for the RTBP (Szebehely  ) including radiation pressure and oblateness of the smaller primary is written in accordance with Sharma & Subba Rao  , Sharma  as
, where m1 and m2 are masses of larger and smaller primary respectively.
The perturbed mean motion, n of the primaries due to oblateness is given by
AE, AP being the dimensional equatorial and polar radii of the smaller primary and R is the distance between the primaries. The two terms occurring in Ω due to oblateness of the smaller primary were introduced by Sharma and Subba Rao  .
3. Computation of Halo Orbits
For the computation of the halo orbits, the origin is transferred to the Lagrangian points L1 and L2. The transformation is given by
The equation of motion can be written as
The upper sign in the above equations depicts the Lagrangian point L1 and the lower sign corresponds to L2. The distance between these Lagrangian points and the smaller primary is considered to be the normalized unit as in Koon et al.  and Tiwary and Kushvah  .
The usage of Legendre polynomials can result in some computational advantages, when non-linear terms are considered. The non-linear terms are expanded by using the following formula given in Koon et al.  :
The above formula is used for expanding the non-linear terms in the equations of motion. The equations of motion after substituting the values of the non-linear terms and by some algebraic manipulations by defining a new variable cm after expanding up to m = 2 become
Neglecting the higher-order terms in Equations ((4)-(6)), we get
It is clear that the z-axis solution obtained by putting X = Y = 0 does not depend upon X and Y and c2 > 0. Hence we can conclude that the motion in Z-direction is simple harmonic. The motion in XY-plane is coupled.
A fourth degree polynomial is obtained which gives two real and two imaginary roots as eigenvalues:
The solution of the linearized Equations ((7)-(9)), as in Thurman and Worfolk  , is
are arbitrary constants. Since we are concentrating on constructing a halo orbit, which is periodic, we consider.
Frequency and amplitude terms are introduced to find solution through Lindstedt-Poincaré method. The solution of the linearized equations is written in terms of the amplitudes (Ax and Az) and the phases (In-plane phase, ϕ and out-of-plane phase, ψ) and the frequencies (λ and), with an assumption that, as
3.1. Amplitude Constraint
For halo orbits, the amplitudes and are constrained by a non-linear algebraic relationship given by Richardson  :
where l1 and l2 depend upon the roots of the characteristic equation of the linear equation. The correction term arises due to the addition of frequency term in Equation (9).
Hence, any halo orbit can be characterized by specifying a particular out-of-plane amplitude of the solution to linearized equations of motion. Both analytical and numerical methods employ this scheme. From the above expression, we can find the minimum permissible value of Ax to form the halo orbit (Az > 0).
3.2. Phase Constraint
For halo orbits, the phases ϕ and ψ are related as
When Ax is greater than certain value, the third-order solution bifurcates. This bifurcation is manifested through the phase-angle constraint. The solution branches are obtained according to the value of m. For m = 1, Az is positive and we have the northern halo (z > 0). For m = 3, Az is negative and we have the southern halo (z < 0).
3.3. Lindstedt-Poincaré Method
Lindstedt-Poincaré method involves successive adjustments of the frequencies to avoid secular terms and provides approximate periodic solution. The equations of motion with non-linear terms up to third-order approximation as in Richardson  and Thurman and Worfolk  are
A new independent variable τ = ωt is introduced, where refers to.
The values of are chosen in such a way that the secular terms are removed with successive approximations. Most of the secular terms are removed by the following assumptions:
The coefficients and are given in Appendix. The equations of motion are
We continue the perturbation analysis by assuming the solutions of the form
where is a small parameter.
Substituting Equations ((13)-(15)) in Equations ((10)-(12)), and equating the coefficients of the same order of, we get the first-, second- and the third-order equations, respectively.
3.3.1. First-Order Equations
The first-order equations are obtained by taking the coefficients of the term. These are
The periodic solution to the above equations is
3.3.2. Second-Order Equations
Coefficients of the term provide the second-order equations
To remove the secular terms, we set ω1 = 0. The particular solution of the second-order equations are obtained with the help of Maxima software as
The coefficients are given in the Appendix.
3.3.3. Third-Order Equations
By collecting the coefficients of, we get the third-order equations
From the above equations, it is not possible to remove the secular terms by setting ω2 = 0.
Hence, the phase relation is used to remove the secular terms.
The solution of the third-order equation is obtained as
The coefficients are given in the Appendix. Thus, the third-order analytical solution is obtained.
3.3.4. Final Approximation
The mapping, and , will remove from all the equations. We now combine the solutions up to third-order to get the final solution as
The coefficients are given in the Appendix.
4. Analytical Construction of Orbit
4.1. Input Parameters
The input parameters to generate a halo orbit include mass ratio (μ), mean distance between the two primaries, radiation pressure (q or ε), oblateness coefficient (A2), amplitude in z-direction (Az). Once all the input parameters are given, the co-ordinates for the halo orbit are computed.
4.2. Halo Orbit for Classical Case
The halo orbits are generated with the following data for Sun-Earth L1:
Mass ratio = 0.000003, mean distance = 149,600,000 km, q = 1, A2 = 0.
The amplitude in z-direction is taken to be Az= 110,000 km (Amplitude of the ISEE-3 mission around the Sun-Earth L1 point). Figures 2-5 show the different views of the halo orbits for the classical case.
4.3. Halo Orbit for Different Values of Radiation Pressure
The halo orbits are generated for different values of radiation pressurein the vicinity of L1. Figure 6 shows the variation in orbit for different values of radiation pressure. The orbital period increases and the orbits increase in size with the increase in the radiation pressure of the more massive primary.
Table 1 shows the change in the distance of the Lagrangian point L1 due to radiation pressure. It is noticed that the Lagrangian point L1 moves towards the more massive primary with increase in radiation pressure. It is also noticed that L1 moves further towards the more massive primary for q < 0.8.
5. Numerical Computation of Halo Orbits
The third-order analytical solution provides a good initial estimate to compute the halo orbits numerically. Analytical approximation must be combined with a suitable numerical method to obtain accurate halo orbits around the Lagrangian points. The method of differential correction is a powerful application of Newton’s method that employs the State Transition Matrix (STM) to solve various kinds of boundary value problems. In order to propagate an orbit in RTBP, the equations of motion are numerically integrated by using adaptive fourth-order Runge-Kutta method, until the desired orbit is obtained.
The initial guess for finding the numerical solution is taken from the analytic solution. Since the halo orbits are symmetric about xz-plane (y = 0), and they intersect this plane perpendicularly i.e., the initial state vector takes the form
The final state vector which lies on the same plane, takes the form as given below and crosses the xz-plane perpendicularly:
Figure 2. X vs Y for Sun-Earth system.
Figure 3. Y vs Z for Sun-Earth system.
Then the orbit will be periodic with period The state transition matrix at can be used to adjust the initial values of a nearby periodic orbit. The equations of motion are integrated until y changes sign. Then the step size is reduced and the integration goes forward again. This is repeated until y becomes almost zero, and the time at this point is defined to be. The tolerance is considered to be around 10−12. The orbit is considered periodic if and are nearly zero at. If this is not the case, and can be reduced by correcting two of the three initial conditions and then integrate the equations of motion again.
5.1. Numerically Generated Halo Orbits
The numerically generated orbits which are more accurate, are usually considered for mission purposes. Figure 7
Figure 4. X vs Z for sun-earth system.
Figure 5. 3D view for sun-earth system.
Figure 6. Analytically generated halo orbits [A2 = 0].
Figure 7. Numerically generated halo orbits for different values of q [A2 = 0].
Table 1. Variation in the location of Lagrangian point L1 with radiation pressure q.
shows the numerically generated orbit due to variation in radiation pressure. The numerically generated obits show a close resemblance to the analytically generated orbits. Figure 8 shows the comparison between them.
The analytical solution which includes the radiation pressure of more massive primary and oblateness of the smaller primary can improve the accuracy over the existing solution as given by Tiwary and Kushvah  .
The time period of SOHO orbit was around 179 days and its dimensions were approximately 4.3 × 2.7 × 3.7 m. The computed value of radiation pressure is roughly 0.99997. The flight data of SOHO mission from January to June 2008 was taken from the mission website and is plotted. From Figure 9 and Figure 10, it is clear that the addition of the extra term due to oblateness shows improvement in the orbit computation. However, the deviation in the orbit is due to the actual scenario, where the real satellite in space undergoes the effect of various other perturbations.
Figure 8. Comparison of numerically and analytically generated halo orbits.
Figure 9. Comparison of numerically generated halo orbit with SOHO Mission orbit.
Figure 10. Comparison of numerically generated halo orbit with SOHO mission orbit (enlarged view).
Figure 11. Effect of adding the new term on halo orbit at L1 with q = 1, A2 = 0.000001.
5.2. Improvement over the Existing Results
It can be seen from the Figure 11 that the results show improvement in orbit computation due to the addition of the new term due to oblateness.
The time period of the orbit increases by 30 minutes for q = 1 and A2 = 0.000001 as given in Table 2. Thus, it is clear that the addition of the new term in the potential function has helped in predicting the orbit closer to the real time condition.
5.3. Perturbing Effects on the Orbit
The effects due to radiation pressure and oblateness can be seen from Figure 12. It is seen that with increase in radiation pressure as well as oblateness, the orbit moves towards the source of radiation i.e. the more massive primary. The oblateness coefficient are taken as A2 = 0, 0.000001, 0.000002 and the radiation pressure, q = 1, 0.99, 0.98. With increase in perturbing forces, it is seen that the orbit moves towards the more massive primary.
Since the amplitude Az is bounded by amplitude constraint, the variation in it also affects the size and time period of the orbit, which can be seen in Figure 13. The time period of the orbit also changes with the perturbations and the amplitude Az. For the halo around L1, with increase in radiation pressure and oblateness, the time period increases. It can be noticed in Figure 14 and Table 3.
The amplitudes about the z-axis, Az was taken in a range of 80,000 to 140,000 km. There is a minimum Ax for the feasibility of a halo orbit, when we fix the amplitude Az.
The analytical solution obtained using Lindstedt-Poincaré method provides good initial estimate to generate halo orbits around the collinear point L1 with numerical integration in the photogravitational restricted three-body
Table 2. Effect of new term on time period of the halo orbit.
Figure 12. Variation in orbit due to perturbations around L1 (μ = 0.000003).
Figure 13. Variation in orbit due to Az (μ = 0.000003).
Figure 14. Time period due to q and Amplitude, Az around L1 (μ = 0.000003).
problem, when the smaller primary is considered as an oblate spheroid with its equatorial plane coincident with the plane of motion. A comparison with the orbital data of SOHO mission is carried out. The radiation pressure and oblateness are found to play significant role in improving the accuracy of halo orbits computation. In the vicinity of L1, it is noticed that with increase in radiation pressure, the size of halo orbit increases with the increase in the orbital period. The results for L1 coincide with the results obtained by Eapen and Sharma  . The Lagrangian point L1 moves towards the more massive primary with perturbations and the orbit tends to bend towards the smaller primary with the increase in oblateness. The future aspects of the project involve designing
Table 3. Time period of the halo orbit around L1 (μ = 0.000003).
interplanetary trajectories from Earth to Mars with the help of halo orbits. These trajectories can reduce the cost of the interplanetary and deep space missions.
Coefficients for the second and third-order equations and solution
Submit or recommend next manuscript to SCIRP and we will provide best service for you:
Accepting pre-submission inquiries through Email, Facebook, LinkedIn, Twitter, etc.
A wide selection of journals (inclusive of 9 subjects, more than 200 journals)
Providing 24-hour high-quality service
User-friendly online submission system
Fair and swift peer-review system
Efficient typesetting and proofreading procedure
Display of the result of downloads and visits, as well as the number of cited articles
Maximum dissemination of your research work
Submit your manuscript at: http://papersubmission.scirp.org/