Halo Orbits in the Photo-Gravitational Restricted Three-Body Problem

Saurav Ghotekar, Ram Krishan Sharma

Department of Aerospace Engineering, Karunya Institute of Technology and Sciences, Coimbatore, India.

DOI: 10.4236/ijaa.2019.93020                

Abstract

We study halo orbits in the circular restricted three-body problem (CRTBP) with both the primaries as the sources of radiation. The positioning of the triangular equilibrium points is discussed in a rotating coordinate system.

Keywords

Halo Orbits, Circular Restricted Three-Body Problem, Lagrangian Points, Radiation Pressure, Lindstedt-Poincare Method

1. Introduction

Halo Orbits have been the topic of interest for the past few decades among the research society. It is a periodic, three-dimensional orbit around the Lagrangian points in a three-body problem. The name “Halo” was first used by Farquhar [1] in his doctoral thesis (1968). It was first proposed for the Apollo mission for placing the satellite in an orbit around L2 Lagrangian point in the Earth-Moon system. An enormous contribution was made by Richardson [2] towards the construction of Halo orbits using Lindstedt-Poincare method. His method explained the general idea about Halo orbits. The first and simplest periodic exact solution to the three-body problem is the motion on collinear ellipses found by Euler [3]. Euler studied the motion of the Moon assuming the Earth and the Sun orbit each other on circular orbits and the Moon to be massless. This approach is known as the restricted three-body problem. The Lagrangian points are an equilibrium solution to the restricted three-body problem. Szebehely [4] has an excellent thesis on the Lagrangian points and the restricted three-body problem. The classical model of the three-bodies is not accurate for studying any system in the solar family as it does not consider the perturbing forces such as oblateness and radiation pressure as given by Dutt and Sharma [5]. Some of the significant works in the Photogravitational Restricted Three-Body Problem were done by Radzievskii [6], Bhatnagar and Chawla [7], Eapen and Sharma [8]. The first person to study about the solar radiation pressure was Radzievskii and he proposed that the maximum force due to radiation pressure that acts in the radial direction can be given by

$F_p=F_s\left(1-q\right),$

where q is a pressure reduction factor which is defined in terms of particle radius a, density δ and radiation pressure efficiency factor x as

$q=1-\frac{5.6\times {10}^{-3}}{a\delta }\cdot x$ (c.g.s. units)

$q=1-\epsilon$.

$\epsilon$ is a variable, depending upon the specification of the third body which is also a radiating body.

A star radiates or shines due to nuclear fusion of hydrogen into helium in its core which then releases energy that traverses through the interior of the star and then radiates into space. The eventual existence of star is directly dependent on its mass. The characteristics of a star like its diameter and temperature change over time while its rotation and movement are dependent on its environment. When two or more stars that are bound by gravity and move around each other in a stable orbit in the same plane is called a binary or multi-star system.

This paper considers that both the primaries as source of radiation and the following calculations and results are made. So, there will be two pressure radiation terms viz. q₁ and q₂ representing the radiation of the larger and smaller primaries, respectively.

Lindstedt-Poincare method is used to find the solution for the equations of motion of CRTBP through successive approximations. Numerical methods are used with this analytical solution as the initial value to generate the required Halo orbits. Changes in the Halo orbits for various values of radiation pressure are also analyzed.

2. Circular Restricted Three-Body Problem

In the Circular Restricted Three-Body Problem, two bodies revolve around their center of masses in circular orbits under the influence of their mutual gravitational attraction and a third body moves in a plane defined by the motion of the two revolving bodies. The problem of the motion of the third body is called the Circular Restricted Three-Body Problem, henceforth referred to as CRTBP. In the CRTBP, the mass of the third body is considered negligible compared to the primaries and its presence does not disturb the circular motion of the primaries (Figure 1).

Figure 1. Circular restricted three-body problem in normalized units.

The restricted three-body problem has five equilibrium points, called Lagrangian or Libration Points. These points are the points of zero velocities and the objects placed in these points remains stable. Out of the five Lagrangian points, three are collinear (L₁, L₂, L₃) and the other two points (L₄, L₅) form the equilateral triangle. Although the Lagrangian point is just a point in empty space, its peculiar characteristic is that it can be orbited.

The circular restricted three-body problem (CRTBP) is a special case of the restricted three-body problem where bigger and smaller primaries move in circular motion around their common center of mass. The five equilibrium points of CRTBP are known as the Lagrangian points where the gravitational forces due to two primaries and the centrifugal force on a spacecraft are balanced.

Another assumption is made that the universe is strictly limited to these three bodies under consideration, so there is no other gravitational influence in the system. Though the limitations that have been introduced seem to restrict the problem to the point of impracticality, there are many situations in the solar system to which the CRTBP applies.

3. Equations of Motion

The analytical approximation of three dimensional, periodic orbits about collinear points is performed. The equations of motion are developed from a Lagrangian approach described in Richardson [2]. The analytical approximation solution for the equations of motion is calculated through computer-based program.

The equations of motion for the RTBP is written as (Szebehely [4], Sharma [9], Simmons [10] )

$\ddot{x}-2\stackrel{\dot{}}{y}=\frac{\partial \Omega }{\partial x}$ (2.1)

$\ddot{y}+2\stackrel{\dot{}}{x}=\frac{\partial \Omega }{\partial y}$ (2.2)

$\ddot{z}=\frac{\partial \Omega }{\partial z}$ (2.3)

where

$\Omega =\frac{x^2+y^2}{2}+\frac{\left(1-\mu \right)q_1}{r_1}+\frac{\mu q_2}{r_2}$,

$r_1={\left({\left(x+\mu \right)}^2+y^2+z^2\right)}^{1/2}$, $r_2={\left({\left(x-1+\mu \right)}^2+y^2+z^2\right)}^{1/2}$.

Here the oblateness is neglected so and q₁ and q₂ are the radiation pressure terms of larger and smaller primary,

$\mu =\frac{m_2}{m_1+m_2}$

where m₁ and m₂ are masses of larger and smaller primary respectively.

Here $m_{\text{2}}=\mu$ and $m_1=1-\mu$.

Location of Lagrangian Points

These are the equations for the location of the Lagrangian points including the perturbing force of radiation pressure of both the primaries. The equations for the location of Lagrangian points have been referred from Nishanth and Sharma [11].

For L₁

$\begin{array}{l}{x}^5+2\left(\mu -l\right)x^4+\left(l^2-4l\mu +{\mu }^2\right)x^3+\left(2\mu l\left(l-\mu \right)+\mu q_2-l{q}_1\right)x^2\\ +\left({\mu }^2{l}^2+2\left(l^2{q}_1+{\mu }^2{q}_2\right)\right)x+{\mu }^3{q}_2-l^3{q}_1=0\end{array}$

For L₂

$\begin{array}{l}{x}^5+2\left(\mu -l\right)x^4+\left(l^2-4l\mu +{\mu }^2\right)x^3+\left(2\mu l\left(l-\mu \right)-\mu q_2-l{q}_1\right)x^2\\ +\left({\mu }^2{l}^2+2\left(l^2{q}_1-{\mu }^2{q}_2\right)\right)x-{\mu }^3{q}_2-l^3{q}_1=0\end{array}$

For L₃

$\begin{array}{l}{x}^5+2\left(l-\mu \right)x^4+\left({\mu }^2-4\mu l+l^2\right)x^3+\left(2\mu l\left(\mu -l\right)-\mu q_2-l{q}_1\right)x^2\\ +\left({\mu }^2{l}^2-2\left(l^2{q}_1+{\mu }^2{q}_2\right)\right)x-{\mu }^3{q}_2-l^3{q}_1=0\end{array}$

where $l=1-\mu$.

For L₄

$x=-\mu +\frac{1}{2}$, $y=\frac{+3^{1/2}}{2}$,

For L₅

$x=-\mu +\frac{1}{2}$, $y=\frac{-3^{1/2}}{2}$.

The location of Lagrangian points varies with the variation of radiation pressures. This can be observed after solving the equations. The equations will give the positions of Lagrangian points.

4. Analytical Computation of Halo Orbits

For the computation of the Halo orbits, the origin should be transferred to the Lagrangian points L1 and L2. The transformation is given by

$X=x+\mu \pm \xi -1$,

$Y=y$,

$Z=z$.

Hence the equation of motion can be again written as

$\xi \left(\ddot{X}-2n\stackrel{\dot{}}{Y}\right)=\frac{\partial \Omega }{\partial X},$

$\xi \left(\ddot{Y}+2n\stackrel{\dot{}}{X}\right)=\frac{\partial \Omega }{\partial Y},$

$\xi \ddot{Z}=\frac{\partial \Omega }{\partial Z}.$

where

$\Omega =\frac{x^2+y^2}{2}+\frac{\left(1-\mu \right)q_1}{R_1}+\frac{\mu q_2}{R_2}$,

$R_1=\sqrt{{\left(X\xi +1\mp \xi \right)}^2+{\left(Y\xi \right)}^2+{\left(Z\xi \right)}^2}$,

$R_2=\sqrt{{\left(X\xi \mp \xi \right)}^2+{\left(Y\xi \right)}^2+{\left(Z\xi \right)}^2}$.

upper sign in the above equation depict L1 and lower sign depict L2. The distance between smaller primary and larger primary is considered to be unity.

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 as given by Koon [12].

$\frac{1}{\sqrt{{\left(x-A\right)}^2+{\left(y-B\right)}^2+{\left(z-C\right)}^2}}=\frac{1}{D}\underset{m=0}{\overset{\infty }{{\displaystyle \sum }}}{\left(\frac{\rho }{D}\right)}^m{P}_m\left(\frac{Ax+By+Cz}{D\rho }\right),$

where ${\rho }^2=x^2+y^2+z^2$, $D=A^2+B^2+C^2$.

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 steps by defining a new variable, cm after expanding up to m = 2, we get

$\ddot{X}-2n\stackrel{\dot{}}{Y}-\left(1+2c_2\right)X=\frac{\partial }{\partial X}\underset{m\ge 3}{\overset{\infty }{{\displaystyle \sum }}}{c}_m{\rho }^m{P}_m\left(\frac{X}{\rho }\right),$ (3.1)

$\ddot{Y}+2n\stackrel{\dot{}}{X}+\left(c_2-1\right)Y=\frac{\partial }{\partial Y}\underset{m\ge 3}{\overset{\infty }{{\displaystyle \sum }}}{c}_m{\rho }^m{P}_m\left(\frac{X}{\rho }\right),$ (3.2)

$\ddot{Z}+c_{2}Z=\frac{\partial }{\partial Z}\underset{m\ge 3}{\overset{\infty }{{\displaystyle \sum }}}{c}_m{\rho }^m{P}_m\left(\frac{X}{\rho }\right),$ (3.3)

with

$P_m=\underset{k=0}{\overset{\left[\left(m-2\right)/2\right]}{{\displaystyle \sum }}}\left(3+4k-2\right)P_{m-2k-2}\left(\frac{X}{\rho }\right),$

$c_m=\frac{1}{{\xi }^3}\left[\frac{{\left(-1\right)}^{m}q\left(1-\mu \right){\xi }^{m+1}}{{\left(1\mp \xi \right)}^{m+1}}+{\left(\pm 1\right)}^m\left[\mu +\frac{3\mu A_2}{2{\xi }^2}-\frac{9\mu A_2}{2{\xi }^3}\right]\right].$

Neglecting the non-linear higher order terms in Equations (3.1)-(3.3), we get

$\ddot{X}-2n\stackrel{\dot{}}{Y}-\left(1+2c_2\right)X=0,$ (3.4)

$\ddot{Y}+2n\stackrel{\dot{}}{X}+\left(c_2-1\right)Y=0,$ (3.5)

$\ddot{Z}+c_{2}Z=0.$ (3.6)

It is clear that the z-axis solution, obtained by putting X = Y = 0, does not depend upon X and Y and c₂ > 0. Hence, we can conclude that the motion in Z-direction is simple-harmonic. Meanwhile, the motion in XY-plane is coupled.

A fourth-degree polynomial is obtained which gives two real and two imaginary roots as Eigen values.

$\alpha =\sqrt{\frac{\left(c_2-1\right)+\sqrt{9c_2^2-8n^2{c}_2}}{2}},$

$\lambda =\sqrt{\frac{\left(c_2-2\right)-\sqrt{9c_2^2-8n^2{c}_2}}{2}}.$

Thus, the solution to the linearized Equations (3.4)-(3.6) will be of the given form as derived by Thurman and Worfolk [13] as

$X\left(t\right)=A_{1_a}{\text{e}}^{\alpha t}+A_{2_a}{\text{e}}^{-\alpha t}+A_{3_a}\cos \lambda t+A_{4_a}\sin \lambda t,$

$Y\left(t\right)=-k_1{A}_{1_a}{\text{e}}^{\alpha t}+k_1{A}_{2_a}{\text{e}}^{-\alpha t}-k_2{A}_{3_a}\sin \lambda t+k_2{A}_{4_a}\cos \lambda t,$

$Z\left(t\right)=A_{5_a}\cos \sqrt{c_2}t+A_{6_a}\sin \sqrt{c_2}t,$

where

$k_1=\frac{2c_2+1-{\alpha }^2}{2\alpha }$ ; $k_2=\frac{2c_2+1+{\lambda }^2}{2\lambda }$.

Here $A_{1_a},A_{2_a},\cdots ,A_{6_a}$ are arbitrary constants. Since we are interested in periodic solutions we take $A_{1_a}=A_{2_a}=0$.

As mentioned earlier, there is a necessity to introduce frequency and amplitude terms to perform the Lindstedt-Poincare method. Then the solution to the linear problem can be written in terms of amplitudes, phases and the frequencies (λ and $\sqrt{c_2}$), as

$X\left(t\right)=-A_x\cos \left(\lambda t+\varphi \right),$

$Y\left(t\right)=k_2{A}_x\sin \left(\lambda t+\varphi \right),$

$Z\left(t\right)=A_z\sin \left(\lambda t+\psi \right).$

Here an assumption is made that $\lambda =\sqrt{c_2}$.

4.1. Amplitude Constraint

When we consider the periodic 3D orbits, the amplitudes A_x and A_z are related by non-linear algebraic relationship given by Richardson [2].

$l_1{A}_x^2+l_2{A}_z^2+\Delta =0.$

where l₁ and l₂ depend upon the roots of the characteristic equation of the linear equation. The correction term, $\Delta ={\lambda }^{\text{2}}-c^{\text{2}}$ arises due to the addition of frequency term in Equation (3.6).

4.2. Phase Constraint

For periodic 3D orbits, the in-plane phase (ϕ) and out-of-plane phase (ψ) are related as

$\psi -\varphi =\frac{m\text{π}}{2},m=1,3$

When A_x 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, A_z is positive and we have the northern Halo (z > 0) and for m = 3, A_z is negative and we have the northern Halo (z < 0).

4.3. Lindstedt-Poincare Method

We can find better approximations to the non-linear problem in a neighborhood of the equilibrium point solutions by using the perturbation techniques of Lindstedt-Poincare. The equations of motion with non-linear terms up to third order approximation as in Richardson [2] and Thurman and Worfolk [13] are

$\begin{array}{l}\ddot{X}-2n\stackrel{\dot{}}{Y}-\left(1+2c_2\right)X\\ =\frac{3}{2}{c}_3\left(2X^2-Y^2-Z^2\right)+2c_4\left(2X^2-3Y^2-3Z^2\right)X+О\left(4\right),\end{array}$

$\ddot{Y}+2n\stackrel{\dot{}}{X}+\left(c_2-1\right)Y=-3c_{3}XY-\frac{3}{2}{c}_4\left(4X^2-Y^2-Z^2\right)Y+О\left(4\right),$

$\ddot{Z}+c_{2}Z=-3c_{3}XZ-\frac{3}{2}{c}_4\left(4X^2-Y^2-Z^2\right)Z+О\left(4\right).$

A new independent variable τ = ωt is introduced as discussed in Equation (5.1). Here on, $\text{'}$ refers to $\frac{\text{d}}{\text{d}\tau }$.

$\omega =1+{\displaystyle \underset{n>1}{\sum }{\nu }^n{\omega }_n},{\omega }_n<1$

We choose the value of ${\omega }_n$ such that the secular terms get removed with successive approximations.

Under the following assumptions ${\omega }_1=0$ and ${\omega }_2=s_1{A}_x^2+s_2{A}_z^2$, most of the secular terms are removed.

The coefficients $s_1$ and $s_2$ are given in Appendix. The equations of motion changes again.

$\begin{array}{l}{\omega }^2{X}^{″}-2\omega n{Y}^{\prime }-\left(1+2c_2\right)X\\ =\frac{3}{2}{c}_3\left(2X^2-Y^2-Z^2\right)+2c_4\left(2X^2-3Y^2-3Z^2\right)X+О\left(4\right),\end{array}$ (3.7)

${\omega }^2{Y}^{″}+2\omega X^{\prime }+\left(c_2-1\right)Y=-3c_{3}XY-\frac{3}{2}{c}_4\left(4X^2-Y^2-Z^2\right)Y+О\left(4\right),$ (3.8)

${\omega }^2{Z}^{″}+c_{2}Z=-3c_{3}XZ-\frac{3}{2}{c}_4\left(4X^2-Y^2-Z^2\right)Z+О\left(4\right).$ (3.9)

We continue the perturbation analysis by assuming the solutions of the form, where $\nu$ is a small parameter which takes care of the non-linearity.

$X\left(\tau \right)=\nu X_1\left(\tau \right)+{\nu }^2{X}_2\left(\tau \right)+{\nu }^3{X}_3\left(\tau \right)+\cdots ,$ (3.10)

$Y\left(\tau \right)=\nu Y\left(\tau \right)+{\nu }^2{Y}_2\left(\tau \right)+{\nu }^3{Y}_3\left(\tau \right)+\cdots ,$ (3.11)

$Z\left(\tau \right)=\nu Z_1\left(\tau \right)+{\nu }^2{Z}_2\left(\tau \right)+{\nu }^3{Z}_3\left(\tau \right)+\cdots .$ (3.12)

We substitute these Equations (3.10)-(3.12) in Equations (3.7)-(3.9) and equating the coefficients of the same order of $\nu$, we get the first, second and third order equations respectively.

4.3.1. First-Order Equations

The analytical solution can be found using mathematical software like Maxima and Mathematica. The first-order equation is obtained by taking the coefficients of the term $\nu$. It is given as

${X^{″}}_1-2{Y^{\prime }}_1-\left(1+2c_2\right)X_1=0,$

${Y^{″}}_1+2{X^{\prime }}_1+\left(c_2-1\right)Y_1=0,$

${Z^{″}}_1+c_2{Z}_1=0.$

The periodic solution to the equations above is

$X_1\left(t\right)=-A_x\cos \left(\lambda t+\varphi \right),$ (3.13)

$Y_1\left(t\right)=k{A}_x\sin \left(\lambda t+\varphi \right),$ (3.14)

$Z_1\left(t\right)=A_z\sin \left(\lambda t+\psi \right).$ (3.15)

where, k = k₂.

4.3.2. Second-Order Equations

To get the second-order equations we segregate the terms of ν²,

${X^{″}}_2-2{Y^{\prime }}_2-\left(1+2c_2\right)X_2={\gamma }_{21},$

${Y^{″}}_2+2{X^{\prime }}_2+\left(c_2-1\right)Y_2={\gamma }_{22},$

${Z^{″}}_2+c_2{Z}_2={\gamma }_{23}.$

where

${\tau }_1=\lambda \tau +\varphi ;{\tau }_2=\lambda \tau +\psi .$

${\gamma }_{21}={\gamma }_1\cos 2{\tau }_1+{\gamma }_2\cos 2{\tau }_2+{\alpha }_1+\left[2A_x{\omega }_1{\lambda }^2+2A_x{\omega }_{1}nk\lambda \right],$

${\gamma }_{22}={\gamma }_3\sin {\tau }_1+{\gamma }_4\sin 2{\tau }_1,$

${\gamma }_{23}={\gamma }_5\sin {\tau }_2+{\gamma }_6[\sin \left({\tau }_1+{\tau }_2\right)+\sin ({\tau }_2-{\tau }_1)].$

We need to set the value ${\omega }_1=0$, to remove the secular terms. The particular solution of the second-order equations is found using mathematical software like Maxima and Mathematica.

$X_2\left(\tau \right)={\rho }_{20}+{\rho }_{21}\cos 2{\tau }_1+{\rho }_{22}\cos 2{\tau }_2,$ (3.16)

$Y_2\left(\tau \right)={\sigma }_{21}\sin 2{\tau }_1+{\sigma }_{22}\sin 2{\tau }_2,$ (3.17)

$Z_2\left(\tau \right)={\kappa }_{21}\sin \left({\tau }_1+{\tau }_2\right)+{\kappa }_{22}\sin \left({\tau }_2-{\tau }_1\right).$ (3.18)

The coefficients are given in Appendix I.

4.3.3. Third-Order Equations

To get the third-order equations we segregate the terms of ν³,

${X^{″}}_3-2{Y^{\prime }}_3-\left(1+2c_2\right)X_3={\gamma }_{31},$

${Y^{″}}_3+2{X^{\prime }}_3+\left(c_2-1\right)Y_3={\gamma }_{32},$

${Z^{″}}_3+c_2{Z}_3={\gamma }_{33}.$

where,

$\begin{array}{c}{\gamma }_{31}=\left[{\nu }_1+2A_x{\omega }_2\lambda \left(k-\lambda \right)\right]\cos 2{\tau }_1+{\gamma }_7\cos {\tau }_1\\ +{\gamma }_8\cos \left(2{\tau }_2+{\tau }_1\right)+{\gamma }_9\cos \left(2{\tau }_2-{\tau }_1\right),\end{array}$

$\begin{array}{c}{\gamma }_{32}=\left[{\nu }_2+2A_x{\omega }_2\lambda \left(k\lambda -1\right)\right]\sin {\tau }_1+{\beta }_3\sin 3{\tau }_1\\ +{\beta }_4\sin \left(2{\tau }_2+{\tau }_1\right)+{\beta }_5\sin \left(2{\tau }_2-{\tau }_1\right),\end{array}$

$\begin{array}{c}{\gamma }_{33}=\left[{\nu }_3+A_z\left(2{\omega }_2{\lambda }^2\right)\right]\sin {\tau }_2+{\delta }_3\sin 3{\tau }_2\\ +{\delta }_4\sin \left(2{\tau }_1+{\tau }_2\right)+{\delta }_5\cos \left(2{\tau }_1-{\tau }_2\right).\end{array}$

From the above equations, it is not possible to remove the secular terms by setting ${\omega }_{\text{2}}=0$. Hence, the phase relation is used here to remove the secular terms.

$\psi =\varphi +p\frac{\text{π}}{2},p=0,1,2,3$

Hence, the solution of the third-order equation will be

$X_3\left(\tau \right)={\rho }_{31}\cos 3{\tau }_1,$ (3.19)

$Y_3\left(\tau \right)={\sigma }_{31}\sin 3{\tau }_1+{\sigma }_{32}\sin {\tau }_1,$ (3.20)

$Z_3\left(\tau \right)=\{\begin{array}{l}{\left(-1\right)}^{\frac{p}{2}}{\kappa }_{31}\sin 3{\tau }_1,p=0,2\\ {\left(-1\right)}^{\frac{p-1}{2}}{\kappa }_{32}\cos 3{\tau }_1,p=1,3\end{array}$ (3.21)

The coefficients are given in Appendix I. Thus, the third-order analytical solution is developed.

4.3.4. Final Approximation

Applying the mapping, $A_x\to \frac{A_x}{\nu }$ and $A_z\to \frac{A_z}{\nu }$, will remove $\nu$ from all

equations. We now combine the solutions up to third order to get the final solution.

$X\left(\tau \right)={\rho }_{20}+{\rho }_{21}\cos 2{\tau }_1+{\rho }_{22}\cos 2{\tau }_2+{\rho }_{31}\cos 3{\tau }_1-A_x\cos \left(\lambda t+\varphi \right)$

$Y\left(\tau \right)={\sigma }_{21}\sin 2{\tau }_1+{\sigma }_{22}\sin 2{\tau }_2+{\sigma }_{31}\sin 3{\tau }_1+{\sigma }_{32}\sin {\tau }_1+k{A}_x\sin \left(\lambda t+\varphi \right)$

$Z\left(\tau \right)=\{\begin{array}{l}{A}_z\sin \left(\lambda t+\psi \right)+{\kappa }_{21}\sin \left({\tau }_1+{\tau }_2\right)+{\kappa }_{22}\sin \left({\tau }_2-{\tau }_1\right)\\ +{\left(-1\right)}^{\frac{p}{2}}{\kappa }_{31}\sin 3{\tau }_1,p=0,2\\ A_z\sin \left(\lambda t+\psi \right)+{\kappa }_{21}\sin \left({\tau }_1+{\tau }_2\right)+{\kappa }_{22}\sin \left({\tau }_2-{\tau }_1\right)\\ +{\left(-1\right)}^{\frac{p-1}{2}}{\kappa }_{32}\cos 3{\tau }_1,p=1,3\end{array}$

The coefficients are given in the Appendix I.

5. Analytical Construction of Halo Orbits

The output from the analytical computation of Halo orbit is used to construct Halo orbits through computer-based codes and is plotted using software such as MATLAB. Input parameters are given such as mass ratio (µ), mean distance between the two primaries, radiation pressures (q₁ and q₂), amplitude in Z-direction (A_z).

Halo orbits are plotted for the classical case of Antares-Sun-Proxima Centauri system around L₁ Lagrangian point. µ = 0.0011314, q₁ = q₂ = 1, A_z = 110,000 km (Figures 2-5).

Figure 2. X vs Y—classical.

Figure 3. X vs Z—classical.

Figure 4. Y vs Z—classical.

Figure 5. 3D halo orbit—classical.

6. Numerical Computation

Differential Correction (DC) schemes use the STM for targeting purposes. One application is the iterative process to isolate a trajectory arc that connects two points in solution space. Differential Correction techniques can be used to quickly obtain a solution with the desired parameters in a wide range of problems. A sufficiently accurate guess for the initial state is always required. Differential corrections are often used to obtain periodic solutions to the non-linear differential equations in the CRTBP. A common assumption, making use of the symmetry in this problem, is that the desired solution is symmetric about the x-z plane.

The initial guess for finding the solution is taken from the analytic solution. Since the halo orbits are symmetric about x-z plane, (y = 0), and they intersect this plane perpendicularly i.e. ( $\stackrel{\dot{}}{x}=\stackrel{\dot{}}{y}=0$) , the initial state vector takes the form,

$X_0\left(t_0\right)={\left(x_0,0,z_0,0,{\stackrel{\dot{}}{y}}_0,0\right)}^{\text{T}}$

The final state vector which again lies on the same plane, takes the form as given below and crosses the x-z plane perpendicularly.

$X_{1/2}\left(t_{1/2}\right)={\left(x_0,0,z_0,0,{\stackrel{\dot{}}{y}}_0,0\right)}^{\text{T}}$

Then the orbit will be periodic with period $T=2t_{1/2}$. The transition matrix at $t_{1/2}$ can be used to adjust the initial values of a nearby periodic orbit. Using Runge-Kutta fourth order method, 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 as $t_{1/2}$. The orbit is considered periodic if $\left|\stackrel{\dot{}}{x}\right|$ and $\left|\stackrel{\dot{}}{z}\right|$ are nearly zero at $t_{1/2}$. If this is not the case, $\stackrel{\dot{}}{x}$ and $\stackrel{\dot{}}{z}$ can be reduced by correcting two of the three initial conditions and integrate again.

6.1. Numerical Computation of Halo Orbits

The following plots are 2D and 3D plots of the halo orbits for perturbed models. The variation in size, shape of the halo orbits is observed by comparing it with the plots of classical case and the effects of radiation pressures on the halo orbit is observed. The numerical computation of halo orbits is referred from Chidambararaj P. and Sharma R. K. [14].

The input parameters are µ = 0.0011314, q₁ = 0.98, q₂ = 0.97 (Figures 6-10).

Figure 6. X vs Y.

Figure 7. Y vs Z.

Figure 8. X vs Z.

Figure 9. 3D halo orbit (q₁ > q₂).

Figure 10. 3D halo orbit (q₂ > q₁).

6.2. Perturbing Effects on the Orbit

The effects due to radiation pressure can be seen from Figure 11 below and it is seen that with an increase in radiation pressure, the orbit moves towards the larger source of radiation irrespective of the size of the primary. The oblateness coefficient are taken as A₂ = 0, μ = 0.0011314 and the two cases of radiation pressures are studied when q₁ = 0.98, q₂ = 0.97 and when q₁ = 0.96, q₂ = 0.97. With increasing perturbing forces, it is seen that the orbit will move towards the primary with larger perturbing force and thus it doesn’t remain Halo orbit. Hence it is to be seen that the value of A_z should be taken care such that the mission objectives are optimally obtained (Table 1).

It has been observed that the orbit length shortens and the width increases closer to the primary whose radiation pressure is larger. The size factor also comes into effect on the Halo orbit. This can be observed from the 3D plots above as both the plots are not mirror images of each other. This tells that the orbit doesn’t vary with the radiation pressure only. The variation of the orbit shape and size is affected by the mass of the primary as we increase the radiation pressure.

Case: When the radiation pressure for both the primaries is equal and not equal to 1.

Around the L1 Lagrangian point for the Sun-Proxima system.

For µ = 0.109509, q₁ = q₂ = 0.98 Halo period = 2.4662.

Figure 11. 3D halo orbit.

Table 1. Radiation pressure parameters, mass ratio and halo period.

7. Summary and Conclusions

The radiation pressures of larger as well as smaller primary are added. The Lindstedt-Poincare method gives approximate good initial values for the numerical solution. The numerical solution makes use of the adaptive Runge-Kutta fourth-order method as integrator. In the vicinity of L1 point, it has been seen that with increase in radiation pressure, the Halo orbit increases in orbital period. It has been seen that the Halo orbits are very sensitive to initial conditions.

Appendix I

Coefficients for the second and third-order equations and solution [10] [11]

${\gamma }_1=\frac{3A_x^2{c}_3}{4}\left[k^2+2\right],$

${\gamma }_2=\frac{3A_z^2{c}_3}{4},$

${\gamma }_3=2A_x{\omega }_1\left(k{\lambda }^2-n\lambda \right),$

${\gamma }_4=\frac{3A_x^2{c}_{3}k}{2},$

${\gamma }_5=2A_z{\omega }_1{\lambda }^2,$

${\gamma }_6=\frac{3A_x{A}_z{c}_3}{2},$

${\alpha }_1=\frac{3c_3}{4}\left[\left(2-k^2\right)A_x^2-A_z^2\right],$

${\rho }_{20}=-\frac{{\alpha }_1}{\left(n^2+2c_2\right)},$

${\rho }_{21}=\frac{4n\lambda {\gamma }_4-{\gamma }_1\left(4{\lambda }^2+n^2-c_2\right)}{{\left(n^2-4{\lambda }^2\right)}^2+\left(n^2+4{\lambda }^2-2c_2\right)c_2},$

${\rho }_{22}=-\frac{{\gamma }_2\left(4{\lambda }^2+n^2-c_2\right)}{{\left(n^2-4{\lambda }^2\right)}^2+\left(n^2+4{\lambda }^2-2c_2\right)c_2},$

${\sigma }_{21}=\frac{4n\lambda {\gamma }_1-{\gamma }_4\left(4{\lambda }^2+n^2+2c_2\right)}{{\left(n^2-4{\lambda }^2\right)}^2+\left(n^2+4{\lambda }^2-2c_2\right)c_2},$

${\sigma }_{22}=\frac{4n\lambda {\gamma }_2}{{\left(n^2-4{\lambda }^2\right)}^2+\left(n^2+4{\lambda }^2-2c_2\right)c_2},$

${\kappa }_{21}=-\frac{{\gamma }_6}{3{\lambda }^2},$

${\kappa }_{22}=\frac{{\gamma }_6}{{\lambda }^2},$

$\begin{array}{c}{\nu }_1=-\frac{3c_3}{2}\left[A_x\left(4{\rho }_{20}+2{\rho }_{21}+k{\sigma }_{21}\right)+A_z\left({\kappa }_{21}+{\kappa }_{22}\right)\right]\\ +\frac{3c_4{A}_x}{2}\left[\left(k^2-2\right)A_x^2+2A_z^2\right],\end{array}$

${\nu }_2=-\frac{3c_3}{2}\left[A_x\left(2k{\rho }_{20}-k{\rho }_{21}-{\sigma }_{21}\right)\right]+\frac{3c_4{A}_{x}k}{2}\left[\left(\frac{3}{4}{k}^2-1\right)A_x^2+\frac{A_z^2}{2}\right],$

${\nu }_3=\frac{3c_3}{2}\left[A_z\left({\rho }_{22}-{\rho }_{20}\right)+A_x\left({\kappa }_{21}+{\kappa }_{22}\right)\right]+\frac{3c_4{A}_z}{2}\left[\left(\frac{1}{2}{k}^2-2\right)A_x^2+\frac{3A_z^2}{4}\right],$

${\gamma }_7=-\frac{3}{2}{c}_3{A}_x\left(2{\rho }_{21}-k{\sigma }_{21}\right)-\frac{1}{2}{c}_4{A}_x^3\left(2+3k^2\right),$

${\gamma }_8=\frac{3}{2}{c}_3\left(k{A}_x{\sigma }_{22}+A_z{\kappa }_{21}-2A_x{\rho }_{22}\right)-\frac{3}{2}{c}_4{A}_x{A}_z^2,$

${\gamma }_9=-\frac{3}{2}{c}_3\left(k{A}_x{\sigma }_{22}-A_z{\kappa }_{22}+2A_x{\rho }_{22}\right)-\frac{3}{2}{c}_4{A}_x{A}_z^2,$

${\beta }_3=\frac{3}{2}{c}_3{A}_x\left({\sigma }_{21}-k{\rho }_{21}\right)-\frac{3}{8}{c}_{4}k{A}_x^3\left(k^2+4\right),$

${\beta }_4=\frac{3}{2}{c}_3{A}_x\left({\sigma }_{22}-k{\rho }_{22}\right)-\frac{3}{8}{c}_{4}k{A}_x{A}_z^2,$

${\beta }_5=\frac{3}{2}{c}_3{A}_x\left({\sigma }_{22}+k{\rho }_{22}\right)+\frac{3}{8}{c}_{4}k{A}_x{A}_z^2,$

${\delta }_3=-\frac{3}{8}\left(A_z^3{c}_4+4A_z{c}_3{\rho }_{22}\right),$

${\delta }_4=-\frac{3}{2}{c}_3\left(A_z{\rho }_{21}-A_x{k}_{21}\right)-\frac{3}{8}{c}_4{A}_x^2{A}_z\left(k^2+4\right),$

${\delta }_5=\frac{3}{2}{c}_3\left(A_z{\rho }_{21}-A_x{k}_{22}\right)+\frac{3}{8}{c}_4{A}_x^2{A}_z\left(k^2+4\right),$

${\rho }_{31}=\frac{6n\lambda \left({\beta }_3+\zeta {\beta }_4\right)-\left(9{\lambda }^2+n^2-c_2\right)\left({\gamma }_3+\zeta {\gamma }_4\right)}{{\left(n^2-9{\lambda }^2\right)}^2+\left(n^2+9{\lambda }^2-2c_2\right)c_2},$

${\sigma }_{31}=\frac{6n\lambda \left({\gamma }_3+\zeta {\gamma }_4\right)-\left(9{\lambda }^2+n^2+2c_2\right)\left({\beta }_3+\zeta {\beta }_4\right)}{{\left(n^2-9{\lambda }^2\right)}^2+\left(n^2+9{\lambda }^2-2c_2\right)c_2},$

${\sigma }_{32}=-\frac{k{\beta }_5}{2n\lambda }.$

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Farquhar, R. (1968) The Control and Use of Libration Point Satellites.
[2]	Richardson, D. (1980) Analytical Construction of Periodic Orbits about the Collinear Points. Celestial Mechanics, 22, 241-253. https://doi.org/10.1007/BF01229511
[3]	Euler, L. (1767) Themoturectilineorum corpore se mutuoattrahentium, Novi commentarii academium scientistum Petropolitanæ 11. Oeuvres, Seria Secunda tome XXV Commentaries Astronomy, 144-151, 286.
[4]	Szebeheley, V. (1967) Theory of Orbits. Academic Press, New York.
[5]	Dutt, P. and Sharma, R.K. (2011) Evolution of Periodic Orbits in the Sun-Mars System. Journal of Guidance, Control and Dynamics, 34, 635-644. https://doi.org/10.2514/1.51101
[6]	Radzievskii, V. (1950) The Restricted Problem of Three Bodies Taking Account of Light Pressure. The Astronomical Journal, 27, 250-256.
[7]	Bhatnagar, K.B. and Chawla, J.M. (1979) A Study of the Lagrangian Points in the Photo-Gravitational Restricted Three-Body Problem. Indian Journal of Pure and Applied Mathematics, 10, 1443-1451.
[8]	Eapen, R.T. and Sharma, R.K. (2014) A Study of Halo Orbits at the Sun-Mars L1 Lagrangian Point in the Photogravitational Restricted Three-Body Problem. Springer, Berlin, 437-441. https://doi.org/10.1007/s10509-014-1951-6
[9]	Sharma, R.K. and Rao, P.S. (1975) Collinear Equilibria and Their Characteristics Exponents in the Restricted Three-Body Problem When the Primaries Are Oblate Spheroids. Celestial Mechanics, 12, 189-201. https://doi.org/10.1007/BF01230211
[10]	Simmons, J.F.L., Mcdonald, A.J.C. and Brown, J.C. (1985) The 3-Body Problem with Radiation Pressure. Celestial Mechanics, 35, 145-187. https://doi.org/10.1007/BF01227667
[11]	Nishanth, P. and Sharma, R.K. (2017) Mars Interplanetary Trajectory Design via Lagrangian Points in the Restricted Three-Body Problem. Technical Report, ISRO.
[12]	Koon, W.S., Lo, M.W. and Marsden, J.E. (2011) Dynamical Systems, the Three-Body Problem and Mission Space Design. Springer, Berlin.
[13]	Thurman, R. and Worfolk, P.A. (1996) The Geometry of Halo Orbits in the Circular Restricted Three Body Problems. Tech. Rep. GCG 95, Univ. Minnesota, Minneapolis.
[14]	Chidambararaj, P. and Sharma, R.K. (2016) Halo Orbits around Sun-Earth L1 in Photogravitational Restricted Three-Body Problem with Oblateness of Smaller Primary. International Journal of Astronomy and Astrophysics, 6, 293-311. https://doi.org/10.4236/ijaa.2016.63025