# Fourier Series Approximations toJ

Authors: Wei Wang; Jianping Yuan; Yanbin Zhao; Zheng Chen; Changchun Chen
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The current paper offers a comprehensive dynamical analysis and Fourier series approximations ofJ
2-bounded equatorial orbits. The initial conditions of heterogeneous families of J
2-perturbed equatorial orbits are determined first. Then the characteristics of two types of J
2-bounded orbits, namely, pseudo-elliptic orbit and critical circular orbit, are studied. Due to the ambiguity of the closed-form solutions which comprise the elliptic integrals and Jacobian elliptic functions, showing little physical insight into the problem, a new scheme, termed Fourier series expansion, is adopted for approximation herein. Based on least-squares fitting to the coefficients, the solutions are expressed with arbitrary high-order Fourier series, since the radius and the flight time vary periodically as a function of the polar angle. As a consequence, the solutions can be written in terms of elementary functions such as cosines, rather than complex mathematical functions. Simulations enhance the proposed approximation method, showing bounded and negligible deviations. The approximation results show a promising prospect in preliminary orbits design, determination, and transfers for low-altitude spacecrafts.


---

## Body

## 1. Introduction

The motion of a particle is known to be integrable by quadratures in a central force field, wherein the potential displays a spherical symmetry in the mathematical expression, written in the following form via a finite or infinite expansion:(1)
V
(
r
)
=
∑
n
=
-
∞
+
∞
a
n
r
n
,
n
∈
ℤ
∖
{
0
}
,
thus the motion of the particle depends on its distance r from the center of the planet only. Such potentials are determined by the density of the planets, usually taken to establish the dynamic models in stellar and galactic systems [1]. Several typical ones are the Hernquist-Newton potential, the Plummer potential, the spherical harmonic potential, the power-law potential, the logarithmic potential, and the Kepler potential [2]. Interestingly, whatever formats the potentials are, and despite the differences in mathematical expressions, analogous rosette-shaped bounded orbits always rise under the specific conditions. Moreover, potentials in some special cases are also equivalent to be spherically symmetric, among which are the radial acceleration problems for orbital transfers in spacecrafts mission design [3–14], and the problems of a particle in equatorial orbits considering the J
2 term studied by some researchers recently [15–18].In general, the total gravitational potential, includingJ
2 perturbative potential and Kepler potential, is axially symmetric, and it works significantly on a low-orbiting particle that moves around an ellipsoidal planet. While the particle’s orbit is confined to the equatorial plane, the J
2 potential reduces to be spherically symmetric, and the problem becomes integrable. Due to its integrability in the sense of Liouville [17], the closed-form solutions of J
2-bounded equatorial orbits were given in terms of elliptic integrals [15, 18], based on which the periodic and pseudo-elliptic bounded relative orbits were obtained. Lately, the solutions of J
2-unbounded equatorial orbits were also studied with the help of elliptic functions, and two types of new unbounded orbits in equatorial plane, namely, pseudo-parabolic orbits and pseudo-hyperbolic orbits, were unveiled [16].The dynamical analysis of theJ
2-bounded equatorial orbits, which may be of deep insight to the problem, is essential but scarce in the literature. Besides, although the elliptic integrals and elliptic functions may present compact closed-form solutions according to the previous studies, these mathematical functions complicate the practical use for preliminary mission design due to the lack of physical insight.In view of this, the present paper first offers a comprehensive study on dynamical characteristics and initial conditions ofJ
2-perturbed equatorial orbits. Then we establish a methodology to obtain arbitrary high-order approximations of the J
2-bounded equatorial orbits by means of Fourier series expansions such as trigonometric functions, since the distance varies periodically as a function of time or polar angle. In this paper, we utilize the Fourier series of order 1, and order 2 via least-squares fitting to the coefficients. The deviations are quantifiable and negligible and are guaranteed to be bounded by imposing a constraint that the approximation solutions coincide with the exact solutions at both endpoints of the half period.As a result, the solutions are expressed in terms of elementary functions such as cosines, rather than complex mathematical functions, and the proposed methods are demonstrated to be effective in simulation.


## 2. Problem Formulation and Mathematic Model

The equatorial geopotential comprisingJ
2 perturbative potential and Kepler potential can be expressed as [19]
(2)
V
(
r
)
=
-
μ
r
[
1
+
J
2
2
(
R
e
r
)
2
]
,
where R
e represents the mean equatorial radius of the planet, J
2 is the second zonal harmonics coefficient, μ is the gravitational parameter, and r denotes the modulus of position vector r.Since we are dealing only with gravitational forces in this paper, the orbit of a particle in a given field does not depend on its mass. Hence, we examine the dynamics of a particle of unit mass, so the quantities such as momentum, angular momentum, and energy, and functions such as the Lagrangian and Hamiltonian, are normally written per unit mass. The LagrangianL of the system is formed by subtracting the potential energy V
(
r
) from the kinetic energy T, where polar coordinates (
r
,
θ
) are mostly conveniently used as follows:
(3)
L
=
T
-
V
(
r
)
=
1
2
[
r
˙
2
+
(
r
θ
˙
)
2
]
+
μ
r
[
1
+
J
2
2
(
R
e
r
)
2
]
.Obviously,θ is an ignorable coordinate from (3); the canonical momentum, also referred to as the first integral constant, yields from
(4)
∂
L
∂
θ
˙
=
r
2
θ
˙
=
h
=
const
.The Euler-Lagrange equation forr, therefore, becomes the radial force equation
(5)
r
¨
-
r
θ
˙
2
=
F
(
r
)
,
where F
(
r
)
=
-
μ
(
1
+
3
J
2
R
e
2
/
2
r
2
)
/
r
2, denoting the acceleration per mass that resulted from the center force or imposed by Kepler and J
2 perturbative forces along r. Also the second integral constant is easily obtained from the Hamiltonian
(6)
ε
=
T
+
V
=
1
2
r
˙
2
+
1
2
(
r
θ
˙
)
2
-
μ
r
-
J
2
μ
R
e
2
2
r
3
=
const
,
as the energy is conserved.Denote the effective potential energy of the radial motion as(7)
V
eff
=
1
2
(
r
θ
˙
)
2
+
V
(
r
)
=
h
2
2
r
2
-
μ
r
-
J
2
μ
R
e
2
2
r
3
.
Substituting (6) into (7) yields
(8)
r
˙
=
±
2
(
ε
-
V
eff
)
.


## 3. Initial Conditions ofJ
2-Perturbed Equatorial Orbits

The particle motion in equatorial plane, subject toJ
2 perturbation, is therefore established in Section 2. It may be checked that if we set r
˙
=
0, a unified equation is obtained according to (8):
(9)
ε
r
3
+
μ
r
2
-
h
2
2
r
+
μ
J
2
R
e
2
2
=
0
,
comprising three main cases that need to be discussed. Note that (9) is a cubic equation if ε
≠
0. To have three real roots, the discriminant of (9) satisfies
(10)
Δ
=
(
-
μ
h
2
12
ε
2
-
μ
3
27
ε
3
-
J
2
μ
R
e
2
4
ε
)
2
+
(
-
h
2
6
ε
-
μ
2
9
ε
2
)
3
≤
0
,
where the equality sign indicates the existence of a couple of identical roots.The three roots are given by(11)
x
1
=
-
μ
3
ε
+
-
μ
h
2
12
ε
2
-
μ
3
27
ε
3
-
J
2
μ
R
e
2
4
ε
+
Δ
3
+
-
μ
h
2
12
ε
2
-
μ
3
27
ε
3
-
J
2
μ
R
e
2
4
ε
-
Δ
3
x
2
=
-
μ
3
ε
+
-
1
+
3
i
2
-
μ
h
2
12
ε
2
-
μ
3
27
ε
3
-
J
2
μ
R
e
2
4
ε
+
Δ
3
+
-
1
-
3
i
2
-
μ
h
2
12
ε
2
-
μ
3
27
ε
3
-
J
2
μ
R
e
2
4
ε
-
Δ
3
x
3
=
-
μ
3
ε
+
-
1
-
3
i
2
-
μ
h
2
12
ε
2
-
μ
3
27
ε
3
-
J
2
μ
R
e
2
4
ε
+
Δ
3
+
-
1
+
3
i
2
-
μ
h
2
12
ε
2
-
μ
3
27
ε
3
-
J
2
μ
R
e
2
4
ε
-
Δ
3
.The other special case isε
=
0, and (9) reduces to be a quadratic equation:
(12)
r
2
-
h
2
2
μ
r
+
J
2
R
e
2
2
=
0
.
Accordingly, the discriminant of (9) satisfies
(13)
Δ
=
(
h
2
2
μ
)
2
-
2
J
2
R
e
2
>
0
.The two roots are(14)
x
1,2
=
h
2
/
2
μ
±
h
4
/
4
μ
2
-
2
J
2
R
e
2
2
.
Assume that the J
2 zonal harmonics is taken as J
2
=
1.08263
×
1
0
-
3 herein, r
0 denotes the initial radius, and (
·
~
) represents the normalized variables: r
~
=
r
/
R
e, t
~
=
t
/
R
e
3
/
μ, h
~
=
h
/
R
e
μ, ε
~
=
ε
/
R
e
/
μ. As the total energy ε and angular momentum h are the two integral constant values, (4) and (6) can be translated in a graphical form in (
r
~
0
~
ε
~
) plane for different values of the fixed h
~, as shown in Figure 1.Figure 1
Representation of motion in the plane(
r
~
0
~
ε
~
) with different values of the fixed h
~.Observe that all the contour lines illustrated in Figure1 pass through the horizontal line ε
=
0. In other words, it implies the existence of heterogeneous families of orbits: bounded orbits for ε
<
0 and unbounded obits for ε
=
0 and ε
>
0; the orbital shape is totally determined by the initial conditions. This is confirmed by Figure 2, in which the phase space portrait is depicted, and the corresponding potential well is outlined.Figure 2
Orbital motions in the phase space(
r
,
r
˙
) and potential energy well.Generally, the curve of the potential energy intersects the constant energy line twice forε
≥
0 and three times for ε
<
0. However, the points associated with the minimum roots of (8) are found to be located below the surface of the planet and are of less physical meaning, so they are ignored in Figure 2. If ε
<
0 and ∂
V
eff
/
∂
r
≠
0, the particle will be trapped in the well, oscillating radially, and the radius range is r
min
⁡
≤
r
≤
r
max
⁡. Here r
min
⁡ and r
max
⁡ correspond to r
˙
=
0 in the phase space (
r
,
r
˙
), referred to as the lower limit and upper limit. As the energy decreases until the effective energy arrives at its local minimum, namely, ∂
V
eff
/
∂
r
=
0, the orbital shape degenerates from a pseudo-ellipse to a circle with a radius r
0, corresponding to the fixed point in Figure 2. Otherwise if the energy increases till ε
≥
0, pseudo-parabolic or pseudo-hyperbolic unbounded orbits will occur, which lead to bifurcations in phase space. In this paper, we only focus on the behaviors of the bounded orbits, whereas the unbounded orbits are left for the further research.


## 4. Fundamental Characteristics ofJ
2-Bounded Equatorial Orbits

A comprehensive analysis on initial conditions of heterogeneous families ofJ
2-perturbed equatorial orbits is discussed. In this section, the focus is J
2-bounded equatorial orbits’ fundamental characteristics. To begin, radial and azimuthal kinetic energy are defined for convenience:
(15)
ε
r
=
1
2
r
˙
2
,
ε
θ
=
1
2
(
r
θ
˙
)
2
.


### 4.1. Pseudo-Elliptic Orbit

For pseudo-elliptic bounded orbit,ε
r
=
0 is satisfied for a certain r and the equation r
˙
=
0 normally has two roots r
min
⁡ and r
max
⁡ provided that ∂
V
eff
/
∂
r
≠
0, known as pericenter and apocenter, between which the particle oscillates radially as it revolves. The radial period, and the azimuthal angle increases by an amount, together with the azimuthal period and the distance r can be computed as [2]
(16)
T
r
=
2
-
ε
[
λ
1
F
(
π
2
,
k
)
+
λ
2
E
(
π
2
,
k
)
+
λ
3
Π
(
π
2
,
α
2
,
k
)
]
,
Δ
θ
=
2
h
λ
2
2
-
ε
F
(
π
2
,
k
)
,
T
θ
=
2
π
Δ
θ
T
r
=
(
π
λ
2
[
λ
1
F
(
π
2
,
k
)
+
λ
2
E
(
π
2
,
k
)
m
m
m
m
m
m
m
m
l
l
l
+
λ
3
Π
(
π
2
,
α
2
,
k
)
]
)
×
(
h
F
(
π
2
,
k
)
)
-
1
,
r
(
θ
)
=
(
r
min
⁡
r
max
⁡
)
×
(
[
F
(
π
2
,
k
)
-
λ
2
-
2
ε
2
h
θ
,
k
]
r
min
⁡
+
(
r
max
⁡
-
r
min
⁡
)
s
n
2
m
m
m
l
l
×
[
F
(
π
2
,
k
)
-
λ
2
-
2
ε
2
h
θ
,
k
]
)
-
1
,
where F
(
π
/
2
,
k
),  E
(
π
/
2
,
k
),  Π
(
π
/
2
,
α
2
,
k
) are the complete elliptic integrals of the first, second, and third kind, respectively. The coefficients are as follows:
(17)
λ
1
=
-
r
max
⁡
r
min
⁡
r
max
⁡
-
r
*
,
λ
2
=
r
min
⁡
(
r
max
⁡
-
r
*
)
,
λ
3
=
r
max
⁡
(
r
max
⁡
+
r
min
⁡
+
r
*
)
r
min
⁡
(
r
max
⁡
-
r
*
)
.
Here, r
* is another root of (9), and
(18)
k
=
(
r
max
⁡
-
r
min
⁡
)
r
*
(
r
max
⁡
-
r
*
)
r
min
⁡
,
α
=
1
-
r
max
⁡
r
min
⁡
,
where sn(·) is the Jacobian elliptic function.Such pseudo-elliptic bounded orbit displays a rosette shape in polar frame that bifurcates from the circular orbit and wraps around it, enclosed in between and tangential to two concentric circles of radiir
min
⁡ and r
max
⁡.The flight time less than one-half radial period, as a function of the polar angle, is given by(19)
t
=
∫
r
min
⁡
r
d
t
d
r
d
r
=
∫
r
min
⁡
r
1
2
(
ε
-
V
eff
)
d
r
=
1
-
2
ε
×
[
ε
-
V
eff
r
+
λ
1
F
(
θ
,
k
)
m
m
m
l
+
λ
2
E
(
θ
,
k
)
+
λ
3
Π
(
θ
,
α
2
,
k
)
ε
-
V
eff
r
]
.One should take modulo ift
>
T
/
2, namely,
(20)
t
=
{
n
T
r
+
∫
r
min
⁡
r
d
t
d
r
d
r
,
t
∈
[
n
T
r
,
(
n
+
1
2
)
T
]
(
n
+
1
)
T
r
-
∫
r
r
max
⁡
d
t
d
r
d
r
,
t
∈
[
(
n
+
1
2
)
T
r
,
(
n
+
1
)
T
]
θ
=
{
n
Δ
θ
+
∫
r
min
⁡
r
d
θ
d
r
d
r
,
t

∈

[
n
T
r
,
(
n
+
1
2
)
T
]
(
n
+
1
)
Δ
θ
-
∫
r
r
max
⁡
d
θ
d
r
d
r
,
t
∈
[
(
n
+
1
2
)
T
r
,
(
n
+
1
)
T
]
,
where n
∈
ℕ.


### 4.2. Critical Circular Orbit

Consider a special case in which the curve of the potential energy is tangent to the constant energy line, indicating two equivalent roots that appear in couple with their radii(
r
)
1
,
2
=
r
0. Equivalently, the potential energy arrives at its minimum. Applying ∂
V
eff
/
∂
r
=
0 yields
(21)
(
r
0
)
1
=
h
2
+
h
4
-
6
J
2
μ
2
R
e
2
2
μ
,
(
r
0
)
2
=
h
2
-
h
4
-
6
J
2
μ
2
R
e
2
2
μ
,
which is in agreement with Humi’s conclusion [20], and (
r
0
)
2 is an extraneous root that should be discarded. The angular velocity and the orbital period of the spacecraft are given by
(22)
ω
=
h
r
0
2
,
T
=
2
π
r
0
2
h
.
Substituting r
0 into ∂
2
V
eff
/
∂
r
2 yields
(23)
∂
2
V
eff
∂
r
2
|
r
=
r
0
=
-
2
μ
r
0
3
-
6
J
2
μ
R
e
2
r
0
5
+
3
h
2
r
0
4
>
0
,
indicating that the J
2-circular orbit is stable according to the theorems of the classical mechanics [21].To detect the orbital evolution when small radial perturbations act on the particle, the substitutionu
=
1
/
r is employed. Thus, the orbital dynamic equation is expressed as
(24)
d
2
u
d
θ
2
+
u
=
J
(
u
)
,
where
(25)
J
(
u
)
=
-
1
h
2
d
d
u
V
(
1
u
)
=
1
h
2
(
3
2
J
2
μ
R
e
2
u
2
-
h
2
u
+
μ
)
.Since motion is stable,u is surely bounded and varies within a small region around u
0. To proceed, we define Δ
u
≜
u
-
u
0, denoting the deviation between perturbed and unperturbed motion. Expanding J
(
u
) in a Taylor series at u
0 and neglecting all terms in J
(
u
) of second and higher orders yield
(26)
J
(
u
)
=
u
0
+
(
u
-
u
0
)
d
J
d
u
0
+
O
[
(
u
-
u
0
)
2
]
,
then the equation of orbital motion around the J
2-circular orbit reduces to
(27)
d
2
Δ
u
d
θ
2
+
Δ
u
=
Δ
u
d
J
d
u
0
.For the sake of compactness, denote(28)
χ
2
=
1
-
d
J
d
u
0
=
2
-
3
J
2
μ
R
e
2
u
0
h
2
,
then a simplified solution to (27) is obtained
(29)
Δ
u
=
u
-
u
0
=
A
cos
⁡
(
χ
θ
)
,
where A is related to the external perturbation imposed on the particle. Through transformation, we have
(30)
r
=
r
0
1
+
A
r
0
cos
⁡
(
χ
θ
)
.
Usually, χ is an irrational number. Hence, the particle will come to a pseudo-elliptic orbit and oscillate radially within a small region.


## 4.1. Pseudo-Elliptic Orbit

For pseudo-elliptic bounded orbit,ε
r
=
0 is satisfied for a certain r and the equation r
˙
=
0 normally has two roots r
min
⁡ and r
max
⁡ provided that ∂
V
eff
/
∂
r
≠
0, known as pericenter and apocenter, between which the particle oscillates radially as it revolves. The radial period, and the azimuthal angle increases by an amount, together with the azimuthal period and the distance r can be computed as [2]
(16)
T
r
=
2
-
ε
[
λ
1
F
(
π
2
,
k
)
+
λ
2
E
(
π
2
,
k
)
+
λ
3
Π
(
π
2
,
α
2
,
k
)
]
,
Δ
θ
=
2
h
λ
2
2
-
ε
F
(
π
2
,
k
)
,
T
θ
=
2
π
Δ
θ
T
r
=
(
π
λ
2
[
λ
1
F
(
π
2
,
k
)
+
λ
2
E
(
π
2
,
k
)
m
m
m
m
m
m
m
m
l
l
l
+
λ
3
Π
(
π
2
,
α
2
,
k
)
]
)
×
(
h
F
(
π
2
,
k
)
)
-
1
,
r
(
θ
)
=
(
r
min
⁡
r
max
⁡
)
×
(
[
F
(
π
2
,
k
)
-
λ
2
-
2
ε
2
h
θ
,
k
]
r
min
⁡
+
(
r
max
⁡
-
r
min
⁡
)
s
n
2
m
m
m
l
l
×
[
F
(
π
2
,
k
)
-
λ
2
-
2
ε
2
h
θ
,
k
]
)
-
1
,
where F
(
π
/
2
,
k
),  E
(
π
/
2
,
k
),  Π
(
π
/
2
,
α
2
,
k
) are the complete elliptic integrals of the first, second, and third kind, respectively. The coefficients are as follows:
(17)
λ
1
=
-
r
max
⁡
r
min
⁡
r
max
⁡
-
r
*
,
λ
2
=
r
min
⁡
(
r
max
⁡
-
r
*
)
,
λ
3
=
r
max
⁡
(
r
max
⁡
+
r
min
⁡
+
r
*
)
r
min
⁡
(
r
max
⁡
-
r
*
)
.
Here, r
* is another root of (9), and
(18)
k
=
(
r
max
⁡
-
r
min
⁡
)
r
*
(
r
max
⁡
-
r
*
)
r
min
⁡
,
α
=
1
-
r
max
⁡
r
min
⁡
,
where sn(·) is the Jacobian elliptic function.Such pseudo-elliptic bounded orbit displays a rosette shape in polar frame that bifurcates from the circular orbit and wraps around it, enclosed in between and tangential to two concentric circles of radiir
min
⁡ and r
max
⁡.The flight time less than one-half radial period, as a function of the polar angle, is given by(19)
t
=
∫
r
min
⁡
r
d
t
d
r
d
r
=
∫
r
min
⁡
r
1
2
(
ε
-
V
eff
)
d
r
=
1
-
2
ε
×
[
ε
-
V
eff
r
+
λ
1
F
(
θ
,
k
)
m
m
m
l
+
λ
2
E
(
θ
,
k
)
+
λ
3
Π
(
θ
,
α
2
,
k
)
ε
-
V
eff
r
]
.One should take modulo ift
>
T
/
2, namely,
(20)
t
=
{
n
T
r
+
∫
r
min
⁡
r
d
t
d
r
d
r
,
t
∈
[
n
T
r
,
(
n
+
1
2
)
T
]
(
n
+
1
)
T
r
-
∫
r
r
max
⁡
d
t
d
r
d
r
,
t
∈
[
(
n
+
1
2
)
T
r
,
(
n
+
1
)
T
]
θ
=
{
n
Δ
θ
+
∫
r
min
⁡
r
d
θ
d
r
d
r
,
t

∈

[
n
T
r
,
(
n
+
1
2
)
T
]
(
n
+
1
)
Δ
θ
-
∫
r
r
max
⁡
d
θ
d
r
d
r
,
t
∈
[
(
n
+
1
2
)
T
r
,
(
n
+
1
)
T
]
,
where n
∈
ℕ.


## 4.2. Critical Circular Orbit

Consider a special case in which the curve of the potential energy is tangent to the constant energy line, indicating two equivalent roots that appear in couple with their radii(
r
)
1
,
2
=
r
0. Equivalently, the potential energy arrives at its minimum. Applying ∂
V
eff
/
∂
r
=
0 yields
(21)
(
r
0
)
1
=
h
2
+
h
4
-
6
J
2
μ
2
R
e
2
2
μ
,
(
r
0
)
2
=
h
2
-
h
4
-
6
J
2
μ
2
R
e
2
2
μ
,
which is in agreement with Humi’s conclusion [20], and (
r
0
)
2 is an extraneous root that should be discarded. The angular velocity and the orbital period of the spacecraft are given by
(22)
ω
=
h
r
0
2
,
T
=
2
π
r
0
2
h
.
Substituting r
0 into ∂
2
V
eff
/
∂
r
2 yields
(23)
∂
2
V
eff
∂
r
2
|
r
=
r
0
=
-
2
μ
r
0
3
-
6
J
2
μ
R
e
2
r
0
5
+
3
h
2
r
0
4
>
0
,
indicating that the J
2-circular orbit is stable according to the theorems of the classical mechanics [21].To detect the orbital evolution when small radial perturbations act on the particle, the substitutionu
=
1
/
r is employed. Thus, the orbital dynamic equation is expressed as
(24)
d
2
u
d
θ
2
+
u
=
J
(
u
)
,
where
(25)
J
(
u
)
=
-
1
h
2
d
d
u
V
(
1
u
)
=
1
h
2
(
3
2
J
2
μ
R
e
2
u
2
-
h
2
u
+
μ
)
.Since motion is stable,u is surely bounded and varies within a small region around u
0. To proceed, we define Δ
u
≜
u
-
u
0, denoting the deviation between perturbed and unperturbed motion. Expanding J
(
u
) in a Taylor series at u
0 and neglecting all terms in J
(
u
) of second and higher orders yield
(26)
J
(
u
)
=
u
0
+
(
u
-
u
0
)
d
J
d
u
0
+
O
[
(
u
-
u
0
)
2
]
,
then the equation of orbital motion around the J
2-circular orbit reduces to
(27)
d
2
Δ
u
d
θ
2
+
Δ
u
=
Δ
u
d
J
d
u
0
.For the sake of compactness, denote(28)
χ
2
=
1
-
d
J
d
u
0
=
2
-
3
J
2
μ
R
e
2
u
0
h
2
,
then a simplified solution to (27) is obtained
(29)
Δ
u
=
u
-
u
0
=
A
cos
⁡
(
χ
θ
)
,
where A is related to the external perturbation imposed on the particle. Through transformation, we have
(30)
r
=
r
0
1
+
A
r
0
cos
⁡
(
χ
θ
)
.
Usually, χ is an irrational number. Hence, the particle will come to a pseudo-elliptic orbit and oscillate radially within a small region.


## 5. Fourier Series Approximations

In the previous section, the closed-form analytical solutions forJ
2-bounded equatorial orbits were derived in terms of elliptic integrals and elliptic functions. However, due to the lack of physical insight regarding these mathematical functions, instead of the elementary functions, it is difficult to put them into practical use for mission design. Here, a Fourier series expansion is resorted for analytical approximations.Equation (5) shows that the radius is an even function of θ, with half radial period Δ
θ. Hence, the generalized Fourier series evaluated on the interval of [
0
,
Δ
θ
] can be written by means of cosines alone [3]
(31)
r
app
(
∞
)
(
θ
)
=
r
min
⁡
1
+
∑
i
=
0
∞
C
i
cos
⁡
(
i
π
(
θ
/
Δ
θ
)
)
,
where C
i are the Fourier coefficients, expressed as
(32)
C
0
=
1
2
Δ
θ
∫
0
Δ
θ
[
r
min
⁡
r
(
θ
)
-
1
]
d
θ
,
C
i
=
1
Δ
θ
∫
0
Δ
θ
[
r
min
⁡
r
(
θ
)
-
1
]
cos
⁡
(
i
π
θ
Δ
θ
)
d
θ
.Consider a limiting case in which the first-order approximation is merely retained in (31) under the assumption that e
J
2
≪
1. As such, the truncated expression of (31) is given by
(33)
r
app
(
1
)
(
θ
)
=
r
min
⁡
1
+
B
0
+
B
1
cos
⁡
(
π
(
θ
/
Δ
θ
)
)
,
where B
0 and B
1 are constants and different from the Fourier coefficients C
i in (31). Assume that the approximation solution coincides with the exact solution at both endpoints of the half period [
0
,
Δ
θ
]; thus, the maximum deviation will remain bounded.Suppose the particle starts at pericenter and ends at apocenter on a short timescale[
0
,
Δ
θ
]; then B
0 and B
1 are calculated as follows:
(34)
B
0
=
-
r
max
⁡
-
r
min
⁡
2
r
max
⁡
,
B
1
=
r
max
⁡
-
r
min
⁡
2
r
max
⁡
.
To seek for a brief format of approximation, denote
(35)
a
J
2
=
r
max
⁡
+
r
min
⁡
2
,
e
J
2
=
r
max
⁡
-
r
min
⁡
r
max
⁡
+
r
min
⁡
,
p
J
2
=
a
J
2
(
1
-
e
J
2
2
)
.
Substituting B
0, B
1, and (35) into (33) yields
(36)
r
app
(
1
)
(
θ
)
=
p
J
2
1
+
e
J
2
cos
⁡
(
π
(
θ
/
Δ
θ
)
)
.
Accordingly, the first-order approximation of flight time as a function of polar angle can be computed by substituting equation (36) into equation d
θ
/
d
t
=
h
/
r
2:
(37)
t
app
(
1
)
(
θ
)
≈
∫
0
θ
[
r
app
(
1
)
(
θ
)
]
2
h
d
θ
=
Δ
θ
p
J
2
2
π
h
(
1
-
e
J
2
2
)
3
/
2
×
{
e
J
2
1
-
e
J
2
2
sin
(
π
θ
/
Δ
θ
)
1
+
e
J
2
cos
⁡
(
π
θ
/
Δ
θ
)
2
arctan
[
1
-
e
J
2
1
+
e
J
2
tan
(
π
θ
2
Δ
θ
)
]
m
m
m
l
-
e
J
2
1
-
e
J
2
2
sin
(
π
θ
/
Δ
θ
)
1
+
e
J
2
cos
⁡
(
π
θ
/
Δ
θ
)
}
.To evaluate the accuracy of approximation, it is necessary to introduce the distance deviation(38)
ε
Δ
(
1
)
=
|
r
app
(
1
)
(
θ
)
-
r
|
r
.
And the flight time deviation is also constructed as
(39)
ε
Δ
t
(
1
)
=
|
t
app
(
1
)
(
θ
)
-
t
|
t
.
Figure 3 shows that for h
~
=
1.1205, J
2
=
1.08263
×
1
0
-
3, and e
J
2
=
4
×
1
0
-
6
≪
1, ε
Δ
(
1
) varies periodically, and the magnitude of ε
Δ
(
1
) is less than 8
×
1
0
-
6, implying the effectiveness of the first-order Fourier series approximation for small e
J
2.Figure 3
Variations of the deviationε
Δ
(
1
) for e
J
2
=
4
×
1
0
-
6
≪
1, and t
~
∈
[
0,5
].Ase
J
2 increases gradually even if e
J
2
>
0.1, the first-order Fourier series approximation of flight time seems to be still effective; the magnitude of ε
Δ
t
(
1
) is less than 3
×
1
0
-
5, as shown in Figure 4. Nevertheless, the radius approximation becomes to be inadequate, especially for e
J
2
>
0.1, as displayed in Figure 5. Naturally it necessitates us to retain high-order terms in the proceedings of approximating.Figure 4
Variations of the deviationε
Δ
t
(
1
) for e
J
2
>
0.1, and θ
∈
[
0,3
].Figure 5
Variations of the deviationε
Δ
(
1
), and t
~
∈
[
0,5
].To that end, we retainm
-
1 orders (m
>
1) of the series, written by
(40)
r
app
(
m
)
(
θ
)
=
r
min
⁡
1
+
∑
i
=
0
m
B
i
cos
⁡
(
i
π
(
θ
/
Δ
θ
)
)
.Once again, assume that the approximation solution coincides with the exact solution at both endpoints of the half-period[
0
,
Δ
θ
]. Sequentially, the coefficients in (40) satisfy the following constraints:
(41)
∑
i
=
0
m
B
i
=
0
,
∑
i
=
0
m
B
i
(
-
1
)
i
=
r
min
⁡
r
max
⁡
-
1
.It is noteworthy that these constraints enablem
-
1 degree-of-freedom. One rational approach is to apply least-squares fitting. The objective function is constructed as
(42)
min
⁡
J
(
B
0
,
B
1
,
…
B
m
)
=
∑
j
=
1
n
∑
i
=
0
m
[
r
app
(
m
)
(
B
i
,
θ
j
)
-
r
(
θ
j
)
]
2
,
where n denotes the quantity of total discrete points in the process of integration.Table1 shows the coefficients, with large e
J
2, obtained by Fourier series approximation of order 2, and the corresponding objective functions are also given. The fixed angular momentum h
~
=
1.1205. One should not neglect the fact that due to the constraint r
app
(
2
)
(
θ
)
>
R
e, e
J
2
<
0.254 must be satisfied.Table 1
Least-squares fitting coefficientsB
0, B
1, B
2 vary with e
J
2
(
2
).
e
J
2
(
2
)
B
0
B
1
B
2
min
⁡
J
0.1111
-
9.99953
e
-
2
9.99966
e
-
2
-
1.25461
e
-
6
3.47383
e
-
8
0.1428
-
1.24991
e
-
2
1.24994
e
-
2
-
2.99121
e
-
6
1.11865
e
-
7
0.1765
-
1.49988
e
-
1
1.49991
e
-
1
-
3.18501
e
-
6
2.36149
e
-
7
0.2121
-
1.74984
e
-
1
1.74988
e
-
1
-
3.91574
e
-
6
5.03017
e
-
7
0.2500
-
1.99789
e
-
1
1.99984
e
-
1
-
4.71406
e
-
6
9.65574
e
-
7Figure6 displays the Fourier series approximation of order 2, via least-squares fitting. In contrast to the result by first-order Fourier series approximations, the maximum deviation is less than 0.025% of the orbital radius r
(
θ
). In view of this, it is suitable for the high-accuracy required missions.Figure 6
Variations of the deviationε
Δ
(
2
) for e
J
2
>
0.1, and t
~
∈
[
0,5
].Figure7 illustrates the particle’s actual orbit (solid line) and orbit by Fourier series approximations (dashed line) in polar reference frame (
θ
~
,
r
~
) with e
J
2
=
0.25. In Figure 8, the left plane shows the evolution of the actual radius and radius of Fourier series approximation of order 2, via least-squares fitting. The right plane presents a magnified view of a local segment.Figure 7
Actual orbit and orbit by Fourier series approximation fore
J
2
>
0.1 in polar reference frame (
θ
~
,
r
~
).Figure 8
Evolution of actual radius and radius by Fourier series approximation fore
J
2
=
0.25.


## 6. Conclusions

The main contribution of this paper is that a framework for approximatingJ
2-bounded equatorial orbits is established with arbitrary high-order Fourier series expansions. Since the distance and time vary periodically as a function of polar angle, the solutions are expressed in terms of elementary trigonometric functions, rather than Jacobian elliptic function and elliptic integrals that lack physical insight into the problem. For Fourier series expansion of second order or higher, the coefficients can be selected via least-squares fitting, and the deviations are guaranteed to be still bounded by imposing a constraint that the approximation solutions coincide with the exact solutions at both endpoints of the half period. Also, the approximation of closed-form relationship for the flight time as a function of the polar angle is given using an analytical approach.The presented approximation method has a potential for space missions, such as novel orbits design, and computational efficiency improvement for long-term low-altitude orbits.



---

*Source: 568318-2014-02-23.xml*