The Exact Solution of the Falling Body Problem in Three-Dimensions: Comparative Study

: Very recently, the system of di ﬀ erential equations governing the three-dimensional falling body problem (TDFBP) has been approximately solved. The previously obtained approximate solution was based on the fact that the Earth’s rotation (ER) is quite slow and hence all high order terms of ω in addition to the magnitude ω 2 R were neglected, where ω is the angular velocity and R is the radius of Earth. However, it is shown in this paper that the ignorance of such magnitudes leads, in many cases, to signiﬁcant errors in the estimated falling time and other physical quantities. The current results are based on obtaining the exact solutions of the full TDFBP-system and performing several comparisons with the approximate ones in the relevant literature. The obtained results are of great interest and importance, especially for other planets in the Solar System or exterior planets, in which ω and / or ω 2 R are of considerable amounts and hence cannot be ignored. Therefore, the present analysis is valid in analyzing the TDFBP near to the surface of any spherical celestial body.


Introduction
In the past decades, some attention was given to the study of the falling body problem [1][2][3] and two-dimensional projectile motion [4][5][6][7][8][9][10][11][12][13][14]. In the literature [1][2][3], the falling body problem was modeled as the vertical motion of a particle near the Earth's gravitational field, i.e., as a problem formulated in one dimension. Hence, the results obtained in [1][2][3] and in [4][5][6][7][8][9][10][11][12][13][14] were only valid in a fixed frame, that is, a non-rotating one. By this, the effect of Earth's rotation (ER) on the falling body problem and also on the projectile motion was not considered. In addition, the effect of ER on the three dimensional-falling body problem (TDFBP) was analyzed very recently by El-Zahar et al. [15]. However, the analysis introduced by El-Zahar et al. [15] still has some restrictions where all high order terms of ω in addition to the magnitude ω 2 R were neglected. In this paper, we consider the effect of the ER on the TDFBP without any constraints/restrictions where all of the physical quantities are left arbitrary and therefore the present results generalize the previously published ones in [15]. The physical problem is described as follows. Let → r = (x(t), y(t), z(t)) be the position vector of a particle located at a point Q near the Earth's surface, relative to P (a surface point) as in Figure 1. x y z at P. The x direction points south, the y direction points east, and the z direction points up. The position vector of a point Q relative to P is denoted r  .
The equations of motion in the frame ( , , ) x y z are given by [15][16][17] where g represents the acceleration due to gravity,  denotes the Earth's angular velocity, [0, ]    is the colatitude, and R is the Earth's radius. Moreover, the x , y , and z directions points south, east, and up, respectively. Moreover, assume that a particle is released from rest at the point Q In the literature [15], system (1)- (4) has been approximately solved by neglecting all high order terms of  in addition to the magnitude 2 R  . The present problem may be of practical importance in engineering and space sciences as explained below. Every day, many planes, satellites, and other space objects (such as the Mir Space Station and the International Space Station) orbit the Earth at different altitudes from the surface of the Earth. The Earth's planet is also exposed, annually, to many crashes of some aircrafts and meteorites, which cause a potential danger to vital installations and regions of population density in some countries. Therefore, the accurate identification of the possible location of these falling objects clearly helps in the possibility of finding An inertial frame (X,Y,Z) attached to the centre of the Earth, and a point P on the surface, described by colatitude λ and longitude φ. The local Cartesian frame (x,y,z) at P. The x direction points south, the y direction points east, and the z direction points up. The position vector of a point Q relative to P is denoted → r .
where g represents the acceleration due to gravity, ω denotes the Earth's angular velocity, λ ∈ [0, π] is the colatitude, and R is the Earth's radius. Moreover, the x, y, and z directions points south, east, and up, respectively. Moreover, assume that a particle is released from rest at the point Q = (x 0 , y 0 , z 0 ); hence, the initial conditions (ICs) are: .
In the literature [15], system (1)-(4) has been approximately solved by neglecting all high order terms of ω in addition to the magnitude ω 2 R. The present problem may be of practical importance in engineering and space sciences as explained below. Every day, many planes, satellites, and other space objects (such as the Mir Space Station and the International Space Station) orbit the Earth at different altitudes from the surface of the Earth. The Earth's planet is also exposed, annually, to many crashes of some aircrafts and meteorites, which cause a potential danger to vital installations and regions of population density in some countries. Therefore, the accurate identification of the possible location of these falling objects clearly helps in the possibility of finding them successfully; hence, the importance of the current study. So, the objective of this paper is to solve the present system, taking into account the contributions of all terms without any constraints/restrictions on any of the physical quantities.

Coordinates of the Falling Point
Let G be the point at which the falling body hits the xy-plane, see Figure 2. At such point, the vertical displacement z(t) vanishes. Suppose that T is the time required such that z(T) = 0; then, the Cartesian coordinates of the falling point G in the local frame (x, y, z) can be expressed as (x G , y G , z G ) = (x(T), y(T), 0). In addition, the Polar coordinates (ρ, δ) of G are: In Appendix C, the coordinates of the falling point at various special cases are introduced. In Appendix C, the coordinates of the falling point at various special cases are introduced.

Characteristics of the Falling Point
In this section, we highlight some of the characteristics of the particle's falling point G. In particular, G lies on the Earth's surface when some conditions are satisfied. As 0   , G belongs to the Earth's surface, as shown from Figure 2. However, 0 The following lemmas address these issues.
At such point G, the component ( ) z t vanishes at a particular value = t T . Hence, G x and G y are, respectively, given from (27) and (28) by The particle, at rest, released from Q in the local Cartesian frame (x, y, z) hits the xyplane at a point G. The angle between PG and the x-axis is denoted δ.

Characteristics of the Falling Point
In this section, we highlight some of the characteristics of the particle's falling point G. In particular, G lies on the Earth's surface when some conditions are satisfied. As ρ → 0 , G belongs to the Earth's surface, as shown from Figure 2. However, ρ → 0 implies x G = y G = 0. The following lemmas address these issues. Lemma 1. If the particle is released from a point (0, 0, z 0 ) in the system (x, y, z) such that λ ∈ {0, π}, then G belongs to the Earth's surface.
Proof of Lemma 1. Substituting x 0 = 0 and y 0 = 0 into (19)-(21) and implementing the corresponding k i (i = 1,2,3,4,5) in (14)-(18), we have: At such point G, the component z(t) vanishes at a particular value t = T. Hence, x G and y G are, respectively, given from (27) and (28) by where T is a non-trivial solution of the following equation (z(T) = 0): It is clear from Equations (30) and (31) that x G = y G = 0 when sin 2λ = 0 and sin λ = 0. These two equations imply that λ ∈ {0, π} which completes the proof. In addition, Equation (29) reduces to: As a confirmation of this Lemma, we have from Appendix C (1 and 2) that ρ = When (x 0 , y 0 ) → (0, 0) , we observe that ρ → 0 ; hence, G lies on the Earth's surface. Proof of Lemma 2. Let us first introduce the relation between the systems (x, y, z) and (X, Y, Z) described in Figure 1. It was shown in [16,17] that the following transformations hold after a time t: From Lemma 1, it was proved that T = 2z 0 g when x 0 = y 0 = 0 and λ ∈ {0, π}, where G = (x G , y G , z G )= (0,0,0) in the local frame (x, y, z). Accordingly, the Cartesian coordinates (X G , Y G , Z G ) of G in the inertial frame (X, Y, Z) are: Substituting x G = y G = z G = 0 into Equations (37)-(39) yields: At λ = 0, we obtain: which are the Cartesian coordinates of the North pole in the frame (X, Y, Z). Moreover, we have from Equations (40)-(42) at λ = π that: which are equivalent to the Cartesian coordinates of the South Pole.

Results and Discussions
This section presents the main difference between the current exact results and those approximately obtained by El-Zahar et al. [15]. Besides, the two Lemmas presented in the previous section will be numerically and graphically validated in this discussion. First of all, the falling time T is approximately obtained in [15] by T Approx. = 2z 0 g ∀ λ ∈ [0, π]. At λ = 0, π, it is noted from Appendix C (1 and 2) that the falling time T is obtained by T = 2z 0 g which agrees with the corresponding expression in [15]. In such special cases, there is no difference between the current exact formula of the falling time T and the corresponding approximate one in [15]. However, when λ takes other values, i.e., λ ∈ (0, π), there are differences in the calculations of T. In these cases, the exact values of T (T Exact ) can be calculated by solving Equation (21) as z(T Exact ) = 0; hence, the error formula of T is defined as: Table 1 presents the comparisons between the current values of the falling time T and the corresponding results in [15] at twenty different cases of the colatitude λ and the release point Q. The results reveal that the obtained error increases as the colatitude λ increases, especially at the equator of the Earth (λ = π/2). In [15], the approximate Cartesian coordinates (x G ) Approx. and (y G ) Approx. of the falling point G were given by: and the approximate polar coordinates (ρ) Approx. and (δ) Approx. of G were obtained as [15]: Therefore, the errors in x G , y G , ρ, and δ at λ = 0 are given by: Similarly, at λ = π, we have: (56) Table 2 presents the obtained errors in the Cartesian coordinates of the falling point G(x G , y G , 0) in the frame (x, y, z) and the corresponding errors in the Polar coordinates (ρ, δ) at several points Q(x 0 , y 0 , z 0 ) when λ = 0, π. The results of Table 2 indicate that the obtained errors in ρ and δ are of the amount at λ = 0 (North Pole) and λ = π (South Pole). In addition, the errors increase as the height z 0 increases above the Earth's surface. While the obtained errors in ρ and δ in Table 2 seem very small at these particular values of the colatitude λ, it will be declared that the errors in x G and y G , ρ, and δ are significant at other values of λ ∈ (0, π). Repeating the analysis above at λ = π/6, λ = π/4, λ = π/3, and λ = π/2, we calculated the errors of x G and y G , as introduced in Table 3. Table 2. The obtained errors of the Cartesian coordinates of the falling point G(x G , y G , 0) in the frame (x, y, z) and the corresponding errors in the polar coordinates (ρ, δ) at several points Q(x 0 , y 0 , z 0 ) when λ = 0, π.  Table 3. The obtained errors of x G and y G at several release points Q(x 0 , y 0 , z 0 ) when λ = π/6, π/4, π/3, π/2.

Case Colatitude (λ) (Radian) Cartesian Coordinates of the Release Point Q Error of Cartesian Coordinates of G
x 0 (km) y 0 (km) z 0 (km) Error of x G (m) Error of y G (m) The results show that the errors in x G and y G increases as λ, x 0 , y 0 , and z 0 increase. For example, when the release point Q(x 0 , y 0 , z 0 ) takes the values Q(10,10,10), the errors in x G are found as 14.9922 m at λ = π/6, 17.3067 m at λ = π/4, 75.9343 m at λ = π/3. However, the obtained maximum errors in y G are 1.4469 m when λ = π/3. Moreover, at λ = π/2 (the equator of Earth), the obtained errors in x G are small for the five cases of Q(x 0 , y 0 , z 0 ), while the obtained error in y G becomes significant and considerable, equaling 246.6514 m for Q(50,50,50). Of course, the significant errors in the Cartesian coordinates x G and in y G imply significant errors in the polar coordinates ρ and in δ. The above calculations implemented the real data of the Earth's planet and can also be extended to cover other planets in the Solar System. Finally, the advantage of the present analysis over the previously published one is observable, especially when applied to specific planets of greater radius R and greater angular velocity ω than those of the Earth.
In order to confirm such conclusion, the current exact solutions are applied for all planets in our Solar System, Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, and Neptune. Astronomers usually divide the planets into two groups, which are known as the inner planets (Mercury, Venus, Earth, and Mars) and the outer planets (Jupiter, Saturn, Uranus, and Neptune) [18]. It will be indicated below that the obtained errors for the falling time T and the landing positions (x G , y G , 0) are bigger for the outer planets than those of the inner planets.
In Table 4, the physical data for all planets are introduced. Before launching to the main results, it can be seen from Table 4 that the value of angular velocity ω of any outer planet (Jupiter, Saturn, Uranus, and Neptune) is greater than the value of any inner planet. This means that the outer planets rotate more speedily about their axes than the inner planets. In another word, the Coriolis effect is expected to affect the calculated errors for outer planets more than it does for the inner planets. In addition, the radii of the outer planets are greater than the radii of the inner planets. So, it is expected that the differences between the present exact solutions and the approximate ones [15] will greatly appear in the cases of outer planets. Our point of view can be confirmed through the calculations in Tables 5 and 6 for λ = π/6 and λ = π/4, respectively. Such calculations are based on assuming that a projectile, at rest, is released from the point Q(50,50,50). The results in Tables 5 and 6 reveal that the obtained errors of T for outer planets are greater than those of inner ones. In addition, for the obtained errors of y G , although small for all planets, we find that the corresponding errors of x G are significant, especially, for those outer planets. Therefore, the current exact analysis may be more helpful than the approximate one, particularly when studying the falling body problem near to the surface of a planet with greater values of ω and R than those of Earth, as it already appeared in cases of Jupiter, Saturn, Uranus, and Neptune. Finally, we are in a period of new planets, new space bodies. For this reason, the current results may also be valid for new planets, which may be discovered later.

Conclusions
In this paper, the three dimensional model of the falling body problem near the Earth's surface was analyzed, taking into account the contribution of all physical quantities. Moreover, unlike the previous published paper [15], the values of the involved parameters were left arbitrary without any restrictions/constraints. Thus, the exact solutions of the full-three coupled equations of motion were obtained in terms of the release point Q(x 0 , y 0 , z 0 ), the colatitude λ, the angular velocity ω, and the radius R of Earth. The present analysis was based on the Laplace transform as a method of solution. Furthermore, the properties of the falling point in the rotating frame and the original inertial frame were given by means of two lemmas. In addition, the conditions at which the falling point belongs G= (x G , y G , z G ) to the Earth's surface were theoretically proven. The effectiveness and validity of the present analysis over the published one in the literature [15] arises when applied to other planets of greater radius R and angular velocity ω than those of the Earth. Hence, the present analysis is applicable to analyzing the TDFBP near enough to the surface of any spherical celestial body.  Acknowledgments: The authors would like to thank the referees for their valuable comments and suggestions which helped to improve the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Laplace Transform and Its Properties
The LT of a function u(t) is defined as [19] L u(t) = ∞ 0 e −st u(t)dt = U(s), where the LT of u(t) is said to exist if the integral (A1) converges. The sufficient conditions for existence of the LT are detailed in [8], see pages (1)(2). Moreover, the LT of the n-derivative of u(t) is given by [19].