Air Drag Effects on the Missile Trajectories

Abstract

The equations of motion of a missile under the air drag effects are constructed. The modified TD88 is surveyed. Using Lagrange's planetary equations in Gauss form, the perturbations, due to the air drag in the orbital elements, are computed between the eccentric anomalies of the burn out and the reentry points [𝐸_𝑏𝑜,2𝜋−𝐸_𝑏𝑜], respectively. The range equation is expressed as an infinite series in terms of the eccentricity e and the eccentric anomaly E. The different errors in the missile-free range due to the drag perturbations in the missile trajectory are obtained.

1. Introduction

The ballistic missile problem is typically represented in astrodynamics as a simplified planar two-body problem with a nonrotating Earth. The missile is only guided during the relatively brief initial and final powered phases of flight, while most of its course, the free flight range, is subsequently governed by the laws of orbital mechanics and ballistics. From this information we can usually determine the range of a projectile, the initial flight path angle, the orbit eccentricity and semimajor axis, and its true/eccentric anomaly angles at launch and impact points. It is of interest to expand upon this basic model in order to add realism to the problem.

In a previous work, Abd El-Salam and Abd El-Bar [1] computed the different errors in the free range of the ballistic missiles taking into account an oblate Earth model retaining the zonal harmonics of the geopotential up to 𝐽₄. They obtained explicit expressions for the errors in the missile range due to the in-orbit plane and out-of-orbit plane changes in the missile range. In this work, the author aims to compute the effect of the atmospheric drag perturbations solely based on the usage of the TD88 density model on the missile trajectory. On one hand, the ballistic missile trajectories are very sensitive to any kind of perturbations, but on the other hand the missile trajectory takes approximately 30 minutes (nearly half an only one orbital period of a low Earth orbit, in brief LEO). Therefore, these perturbations can be safely added up linearly to the gravity perturbations on the missile trajectories obtained previously by the authors mentioned above.

The problem of the orbital motion of celestial bodies in a resistive medium can be traced back to Newton's Principle Newton, (bk. II, Sections I/IV) [2] and Laplace's analysis of a possibly finite velocity of gravitation Laplace [3]. However, with the advent of the space age, there arose the need of precisely predicting the orbital behavior of the newly launched bodies. Compared to the gravitational perturbations, the adequate allowance for atmospheric drag turned out to be a problem for unforeseen complexity Brouwer [4]. Besides the obstacles in the analytical treatment of the drag based on the idealized spherical nonrotating density model, the upper atmosphere revealed itself to behave unexpectedly. The irregular density variations derived from orbital analyses of the first launched satellites were attributed to the influence of the Sun, as they displayed periodical changes of approximately 27 days. Jacchia [5]; King-Hele and Walker [6] introduced the overall dependence of the drag on the angle between the perigee and the Sun. In the subsequent years, together with the advances in theory King-Hele [7], direct evidence has grown indicating that the solar and geomagnetic effects strongly influence the thermospheric density both periodically and randomly Jacchia [8–10], Jacchia and Slowey [11, 12].

The literature on the atmospheric drag effects on artificial satellites is very rich; see for example, King-Hele [7], Hoots and Roehrich [13] Gooding [14], Helali [15] Hoots and France [16, 17], Danielson et al. [18], Fonte et al. [19], Sehnal [20–22], Sehnal and Pospíšilová, [23, 24], Rossi et al. [25], and Abd El-Salam and Sehnal [26].

Although the literature on the perturbed missile trajectories is very limited, we can apply the concepts of atmospheric drag effects on artificial satellites on the perturbed missile trajectories with a modification on the integrating limits.

Isaacson and Vaughan [27] described a method of estimating and predicting ballistic missile trajectories using a Kalman Filter over a spherical, nonrotating Earth. They determined uncertainties in the missile launch point and missile position during flight. Vinh et al. [28] obtained a minimum-fuel interception of a satellite, or a ballistic missile, in elliptic trajectory in a Newtonian central force field, via Lawden's theory of primer vector. McFarland [29] treated ballistic missile problem by modeling spherical Earth, rotation, and addition of atmospheric drag using the state transition matrix. Akgül and Karasoy [30] developed a trajectory prediction program to predict the full trajectory of a tactical ballistic missile. Kamal [31] developed an algorithm that includes detection of cross-range error using Lambert scheme in free space, in the absence of atmospheric drag. Bhowmik and Chandrani [32] investigated the advantages and performance of extended Kalman Filter for the estimation of nonlinear system where linearization takes place about a trajectory that was continually updated with the state estimates resulting from the measurement. Forden [33] described the integration of the three degrees of freedom equations of motion, and approximations are made to the aerodynamic for simulating ballistic missiles. Harlin and Cicci [34] developed a method for the determination of the trajectory of a ballistic missile over a rotating, spherical Earth given only the launch position and impact point. The iterative solution presented uses a state transition matrix to correct the initial conditions of the ballistic missile state vector based upon deviations from a desired set of final conditions. Kamal [35] presented an innovative adaptive scheme which was called “the Multistage Lambert Scheme”. Liu and Chen [36] presented a novel tracking algorithm by integrating input estimation and modified probabilistic data association filter to identify warhead among objects separation from the reentry vehicle in a clear environment.

2. Air Drag Force Model

As an artificial Earth satellite moves through the Earth’s atmosphere, its orbital motion is perturbed by the resisting force. The main effects of the force are to bring out large secular decreases in the elements 𝑎, the semimajor axis, and 𝑒, the eccentricity, that cause the satellite plunges toward Earth. This perturbing influence is one of the most important perturbing forces in the altitude regime from 150 to 600 km, which is called the drag force. The drag force, ^→𝐹_𝐷, per unit mass of the satellite can be represented by^→𝐹_𝐷1=−2𝐶_𝐷𝐴𝑚𝜌𝑉^→𝑉^→, (2.1) where 𝐶_𝐷 is a dimensionless drag coefficient, 𝑉^→ is the velocity of the satellite relative to the atmosphere (called the ambient velocity), 𝐴 represents the effective cross-sectional area of the satellite, which is to be found by averaging all possible projected areas of the satellite onto a plane perpendicular to 𝑉^→, 𝑚 is the mass of the satellite, and 𝜌 is the density of the atmosphere. The density 𝜌 depends on position and time in a very complicated manner. Its variation with height depends on the form of the atmosphere, since it decreases rapidly with slight growing in the altitude, the effect of the flattening of the atmosphere is rather significant.

3. Rotation of the Atmosphere

This effect is so small, so a simple approximation is justified. The differential rotation of the atmosphere complicates the numerical calculation, and it gives no sense in this very small effect. So we will assume that the atmosphere rotates with the same angular velocity of the Earth 𝜎, then𝑉^→=^→𝑣−^→𝜎×^→𝑟, (3.1) where ^→𝑣 is the orbital velocity of the satellite and ^→𝑟 its radius vector. If (𝚤,̂𝚥,̂𝑘) are the base vectors of the rectangular geocentric equatorial frame associated with the spherical coordinates (𝑟,𝜆,𝛿), then ^→𝜎=𝜎×𝑘,𝜎=7.27205217×10⁻⁵ rad./sec. Thus 𝑉²=𝑣²−2^→𝑣⋅^→𝜎×^→𝑟+^→𝜎×^→𝑟²=𝑣²𝑘_𝑅, (3.2) where 𝐼 is the orbital inclination, 𝛿 is the satellite's declination, 𝐡 is the integral of areas, and 𝑘_𝑅 is given by 𝑘_𝑅=1−2𝜎𝐡cos𝐼𝑣²+𝜎²𝑟²cos²𝛿𝑣². (3.3) The third term in 𝑘_𝑅 may usually be neglected, being less than about 1/250, and the second term which is of order of 1/15 may be evaluated at perigee, where the density is much greater than elsewhere on the orbit, thus we may effectively regard 𝑘_𝑅 as a constant.

4. Drag Force Components

Let (𝑟,̂𝑠,̂ℎ) be unit vectors along the radius vector, normal to the radius vector in the orbital plane, and normal to the orbital plane, that is, along to the angular momentum axis, respectively; then the drag force given by (2.1) can be written in the components form as^→𝐹_𝐷=𝑆𝑟̂+𝑇𝑠̂+𝑊ℎ, (4.1) where1𝑆=−2𝐶_𝐷𝐴𝑚𝑘_𝑅𝜌𝑣𝜇𝑝1𝑒sin𝑓,𝑇=−2𝐶_𝐷𝐴𝑚𝑘_𝑅𝜌𝑣𝜇𝑝𝑝𝑟,1−𝜎𝑟cos𝐼𝑊=−2𝐶_𝐷𝐴𝑚𝑘_𝑅𝜌𝑣𝜎𝑟sin𝐼cos𝑓+𝜔. (4.2)

4.1. Modified TD88 Density Model

To integrate the equations of motion, an analytical expression is required for the density of the atmosphere. This density has many irregular and complex variation both in position and time. it is largely affected by solar activity and by the heating or cooling of the atmosphere. This time dependence is very difficult to be introduced in an analytical expression. Since the atmosphere is not actually spherically symmetric but tends to be oblate, we have to allow for this oblateness in any expression for the air density. The complications which come from the accuracy demands may cause the impossibility to use these classical aeronomic empirical models, for example, CIRA 86 Hedien [37], CIRA 72 Oliver [38], and DTM Barlier et al., [39] for an analytical work. In this work, we will use the modified version of the TD88 model first proposed and introduced by Sehnal [20, 22] and Sehnal and Pospíšilová [23] and then developed by Abd El-Salam and Sehnal [26] in terms of the eccentric anomaly 𝐸 as𝜌=𝑓₀𝑓_𝑥𝑘₀⁷³_𝑛=1⁶_𝑗=0⁶_𝜈=0𝑙=0𝐾_𝑛,𝑗𝐴_𝑗exp𝐵_𝑗×𝐺cos𝐸_{𝑐𝜈,𝑛,𝑙}cos(𝜈+𝑙)𝐸+𝐺^𝑐_{𝜈,𝑛,𝑙}cos(𝜈−𝑙)𝐸+𝐺^𝑠_{𝜈,𝑛,𝑙}sin(𝜈+𝑙)𝐸+𝐺^𝑠_{𝜈,𝑛,𝑙}sin(𝜈−𝑙)𝐸, (4.3) where (𝑓₀,𝑓_𝑥,𝑘₀), the main solar and geomagnetic activity parameters affecting the angular density, are given by𝑓₀=𝑎₂+𝑓_𝑚,𝑓_𝑥=1+𝑎₁𝐹_𝑥−𝐹_𝑏,𝑓_𝑚=𝐹_𝑏−60160,𝑘₀=1+𝑎₃𝐾_𝑝−3 (4.4) with 𝐹_𝑥 solar flux measured on 10.7 cm wavelength for a previous day, 𝐹_𝑏 solar flux averaged over three solar rotations (81 days), 𝐾_𝑝 geomagnetic index 3 hours before local time 𝑡, 𝐻_𝑗 the reciprocal of the density scale height 𝐻_𝑗=𝑗(|𝑗−1|!)²/29(𝑗!)², 𝑅_⊕ the equatorial radius of the Earth, 𝜀 the Earth's flattening, 𝐴_𝑗=exp[𝐻_𝑗(120−𝑎+𝑅_⊕)], 𝐵_𝑗=𝐻_𝑗𝑎𝑒, and 𝑎_𝑖,𝐾_𝑛,𝑗,𝑝_𝑛 numerical constants of the model given by Tables 1, 2, and 3.

Table 1 

Constants 𝑎_𝑖.

Table 2 

Constants 𝐾_𝑛,𝑗.

Table 3 

Phases 𝑝_𝑛.

The required nonvanishing coefficients for the model are𝐺^𝑐_{𝜈,𝑛,𝑙}=12𝐶_𝜈𝑔^𝑐_𝑛,𝑙−𝑆_𝜈𝑔^𝑠_𝑛,𝑙,𝐺^𝑐_{𝜈,𝑛,𝑙}=12𝐶_𝜈𝑔^𝑐_𝑛,𝑙+𝑆_𝜈𝑔^𝑠_𝑛,𝑙,𝐺^𝑠_{𝜈,𝑛,𝑙}=12𝐶_𝜈𝑔^𝑠_𝑛,𝑙+𝑆_𝜈𝑔^𝑐_𝑛,𝑙,𝐺^𝑠_{𝜈,𝑛,𝑙}=12𝑆_𝜈𝑔^𝑐_𝑛,𝑙−𝐶_𝜈𝑔^𝑠_𝑛,𝑙, (4.5) where the nonvanishing coefficients 𝑔^𝑐_𝑛,𝑙, 𝑔^𝑠_𝑛,𝑙, 𝐶_𝜈, and 𝑆_𝜈 were computed previously by the author, see Abd El-Salam and Sehnal [26].

5. Integration of the Equations of Motion

The Lagrange's planetary equations are a very well-known system of six first-order differential equations. They are used to evaluate the perturbations in the usual Keplerian orbital elements (𝑎,𝑒,𝐼,Ω,𝜔,𝐸). Gauss modified this system so as to be applicable for all kinds of force especially nonconservative force, see Roy [40],𝑑𝑎=2𝑑𝑡𝑛√1−𝑒²𝑝𝑆𝑒sin𝑓+𝑇𝑟,𝑑𝑒=√𝑑𝑡1−𝑒²,𝑛𝑎𝑆sin𝑓+𝑇cos𝑓+cos𝐸𝑑𝐼=𝑑𝑡𝑟cos𝑓+𝜔𝑛𝑎²√1−𝑒²𝑊,𝑑Ω=𝑑𝑡𝑟sin𝑓+𝜔𝑛𝑎²sin𝐼√1−𝑒²𝑊,𝑑𝜔𝑑𝑡=−cos𝐼𝑑Ω+√𝑑𝑡1−𝑒²𝑟𝑛𝑎𝑒−𝑆cos𝑓+𝑇1+𝑝sin𝑓,𝑑𝑀√𝑑𝑡=𝑛−1−𝑒²𝑑𝜔𝑑𝑡+cos𝐼𝑑Ω𝑑𝑡−𝑆2𝑟𝑛𝑎² (5) with 𝑝=𝑎(1−𝑒²), 𝑛²𝑎³=𝜇, where 𝑛,𝑝,𝜇 are the mean motion in the orbit, the semi-latus rectum of the orbit, and the gravitational parameter, respectively, and (𝑆,𝑇,𝑊) are the drag force components given by (4.2).

To integrate the equations of motion (5), we first do the following simplifying steps:(1)introducing the density given by the modified TD88 model to the drag force components (𝑆,𝑇,𝑊);(2)using 𝑑𝑡=((1−𝑒cos𝐸)/𝑛)𝑑𝐸, to replace the explicit time with the eccentric anomaly 𝐸 as an independent variable;(3)using the binomial theorem and neglecting orders of 𝒪(𝑒⁵) and 𝒪(𝜎𝑒⁴);(4)since the series expansion of exp(𝑥)forall𝑥∈ℝ, the real numbers, is convergent, then we replace exp(𝐵_𝑗cos𝐸) appeared in the density expression with its equivalent series expansion;(5)let us definecos(𝛼±𝛽)𝐸=cos(𝛼+𝛽)𝐸+cos(𝛼−𝛽)𝐸.

5.1. The Increment in the Semimajor Axis

Using the first equation from the equations set (5) and following the above-mentioned simplifying steps, then after some lengthy algebraic manipulations we obtain the change in the semimajor axis as𝑑𝑎𝑑𝐸=−𝐶_𝐷𝐴𝑚𝑘_𝑅𝑎²𝑓₀𝑓_𝑥𝑘₀⁷³_𝑛=1⁶_𝑗=0⁶_𝜈=0⁴_𝑙=0^∞_𝑚=0𝜅=012𝐾_𝑛,𝑗𝐴_𝑗𝑏^𝑎_𝑚×𝐵^𝜅_𝑗cos^𝜅𝐸𝐺𝜅!^{𝑐𝜈,𝑛,𝑙}cos(𝜈+𝑙±𝑚)𝐸+𝐺^𝑐_{𝜈,𝑛,𝑙}cos(𝜈−𝑙±𝑚)𝐸+𝐺^𝑠_{𝜈,𝑛,𝑙}sin(𝜈+𝑙±𝑚)𝐸+𝐺^𝑠_{𝜈,𝑛,𝑙}sin(𝜈−𝑙±𝑚)𝐸, (5.2) where the nonvanishing coefficients are𝑏^𝑎₀=31+4𝑒²+21𝑒64⁴−𝜎cos𝐼𝑛1−𝑒²11−4𝑒²−3𝑒64⁴,𝑏^𝑎₁=32𝑒+4𝑒³,𝑏^𝑎₂=34𝑒²+7𝑒16⁴+𝜎cos𝐼𝑛1−𝑒²14𝑒²+1𝑒16⁴,𝑏^𝑎₃=14𝑒³,𝑏^𝑎₄=7𝑒64⁴+164𝜎cos𝐼𝑛1−𝑒²𝑒⁴ (5.3) using the trigonometric identitycos^𝜅𝐸=𝛿_[]^𝜅/2_𝜅/212^𝜅𝐶^𝜅_[]𝜅/2+12^{)]𝜅−1[(𝜅−1/2}_𝑟=0𝐶^𝜅_𝑟cos(𝜅−2𝑟)𝐸, (5.4) where 𝛿^𝜅/2_[𝜅/2] is the usual Kronecker Delta symbol and [𝜅/2] is the greatest integer function, 𝜅∈ℤ⁺, and 𝐶^𝜅_𝑟=𝜅!/𝑟!(𝜅−𝑟)! is the binomial coefficients. Now (5.4) can be written in the form𝑑𝑎𝑑𝐸=−𝐾_𝐷𝑎²𝑓₀𝑓_𝑥𝑘₀⁷³_𝑛=1⁶_𝑗=0⁶_𝜈=0⁴_𝑙=0^∞_𝑚=0𝜅=0𝐾_𝑛,𝑗𝐴_𝑗𝑏^𝑎_𝑚𝐵^𝜅_𝑗2^𝜅+1×𝛿𝜅!_{[]𝜅/2𝜅/2}𝐶^𝜅_[]𝜅/2𝐺^𝑐_{𝜈,𝑛,𝑙}cos(𝜈+𝑙±𝑚)𝐸+𝐺^𝑐_{𝜈,𝑛,𝑙}cos(𝜈−𝑙±𝑚)𝐸+𝐺^𝑠_{𝜈,𝑛,𝑙}sin(𝜈+𝑙±𝑚)𝐸+𝐺^𝑠_{𝜈,𝑛,𝑙}()𝐸+sin𝜈−𝑙±𝑚_{)][(𝜅−1/2𝑟=0}𝐶^𝜅_𝑟𝐺^𝑐_{𝜈,𝑛,𝑙}((𝐸+𝐺cos𝜅−2𝑟±𝜈+𝑙±𝑚))_{𝑐𝜈,𝑛,𝑙}cos(𝜅−2𝑟±(𝜈−𝑙±𝑚))𝐸+𝐺^𝑠_{𝜈,𝑛,𝑙})±(𝐸+𝐺sin((𝜈+𝑙±𝑚𝜅−2𝑟))_{𝑠𝜈,𝑛,𝑙}sin((𝜈−𝑙±𝑚)±(𝜅−2𝑟))𝐸. (5.5) Now (5.5) can be integrated between the eccentric anomalies of the burn out and the reentry points [𝐸_𝑏𝑜,2𝜋−𝐸_𝑏𝑜], respectively, then we have integrals of the formℐ^𝑐𝜅,𝑟,𝜇=^{2𝜋−𝐸𝑏𝑜}_𝐸𝑏𝑜^∞_𝜅=0𝐵^𝜅_𝑗𝛿𝜅!^{[]𝜅/2𝜅/2}12^𝜅𝐶^𝜅_[]𝜅/21cos𝜇𝐸+2^{𝜅)][(𝜅−1/2}_𝑟=0𝐶^𝜅_𝑟cos𝜅−2𝑟±𝜇𝐸𝑑𝐸, (5.6) from which we can getℐ^𝑐𝜅,𝑟,𝜇=2𝜋^∞_𝜅=0(1/2)𝐵_𝑗^𝜅𝛿𝜅!^{[]𝜅/2𝜅/2}𝐶^𝜅_[]𝜅/21−sin𝜇𝐸_𝑏𝑜𝜇_𝜇≠0+^{)][(𝜅−1/2}_𝑟=0𝐶^𝜅_𝑟𝐸1−sin𝜅−2𝑟±𝜇_𝑏𝑜𝜅−2𝑟±𝜇_{(𝜅−2𝑟±𝜇)≠0}. (5.7) Similarly, we can evaluateℐ^𝑠𝜅,𝑟,𝜇=^{2𝜋−𝐸𝑏𝑜}_{𝐸𝑏𝑜,}^∞_𝜅=0𝐵^𝜅_𝑗2^𝜅𝛿𝜅!^{[]𝜅/2𝜅/2}𝐶^𝜅_[]𝜅/2+2^{)][(𝜅−1/2}_𝑟=0𝐶^𝜅_𝑟()𝐸=cos𝜅−2𝑟sin𝜇𝐸_{2𝜋−𝐸𝑏𝑜𝐸𝑏𝑜,∞𝜅=0}𝐵^𝜅_𝑗2^𝜅𝛿𝜅!^{[]𝜅/2𝜅/2}𝐶^𝜅_[]𝜅/2sin𝜇𝐸+^{)][(𝜅−1/2}_𝑟=0𝐶^𝜅_𝑟sin𝜅−2𝑟±𝜇𝐸, (5.8) which yieldsℐ^𝑠𝜅,𝑟,𝜇=^∞_𝜅=0𝐵^𝜅_𝑗2^𝜅𝜅!−𝛿_[]^𝜅/2_𝜅/2𝐶^𝜅_[]𝜅/2cos𝜇𝐸𝜇^{2𝜋−𝐸𝑏𝑜}_{𝐸𝑏𝑜,,𝜇≠0}−^{)][(𝜅−1/2}_𝑟=0𝐶^𝜅_𝑟cos𝜅−2𝑟±𝜇𝐸𝜅−2𝑟±𝜇^{2𝜋−𝐸𝑏𝑜}_{𝐸𝑏𝑜,,𝜅−2𝑟±𝜇≠0}=0, (5.9) where terms involving the integrals of the type ℐ^𝑠(𝜅,𝑟,𝜇) within the expression (5.5) have no contribution to the problem. as usual to the drag theory of artificial satellites. Now substituting the relevant integrals ℐ^𝑐(𝜅,𝑟,𝜇) from (5.7) into (5.5) yields the solution as(Δ𝑎)_𝐷=−𝐾_𝐷𝑎²𝑓₀𝑓_𝑥𝑘₀⁷³_𝑛=1⁶_𝑗=0⁶_𝜈=0⁴_𝑙=0^∞_𝑚=0𝜅=012𝐾_𝑛,𝑗𝐴_𝑗𝑏^𝑎_𝑚𝐵^𝜅_𝑗×𝐺𝜅!_{𝑐𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈+𝑙+𝑚)+ℐ^𝑐(𝜅,𝑟,𝜈+𝑙−𝑚)+ℐ^𝑐(𝜅,𝑟,−𝜈−𝑙−𝑚)+ℐ^𝑐(+𝐺𝜅,𝑟,−𝜈−𝑙+𝑚)_{𝑐𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈−𝑙+𝑚)+ℐ^𝑐()+𝐺𝜅,𝑟,𝜈−𝑙−𝑚_{𝑐𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,−𝜈+𝑙+𝑚)+ℐ^𝑐(𝜅,𝑟,−𝜈+𝑙−𝑚). (5.10) Using the similar procedure to integrating the rest of (5) between the eccentric anomalies of the burn out and the re entry points [𝐸_𝑏𝑜,2𝜋−𝐸_𝑏𝑜], respectively, we obtained the following explicit expressions for the other orbital elements as follows.

5.2. The Increment in the Eccentricity

(Δ𝑒)_𝐷=−𝐶_𝐷𝐴𝑚𝑘_𝑅𝑎²𝑓₀𝑓_𝑥𝑘₀⁷³_𝑛=1⁶_𝑗=0⁶_𝜈=0⁵_𝑙=0^∞_𝑚=0𝜅=012𝐾_𝑛,𝑗𝐴_𝑗𝑏^𝑒_𝑚×𝐵^𝜅_𝑗𝐺𝜅!^{𝑐𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈+𝑙+𝑚)+ℐ^𝑐(𝜅,𝑟,𝜈+𝑙−𝑚)+ℐ^𝑐(𝜅,𝑟,−𝜈−𝑙−𝑚)+ℐ^𝑐(+𝐺𝜅,𝑟,−𝜈−𝑙+𝑚)_{𝑐𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈−𝑙+𝑚)+ℐ^𝑐(𝜅,𝑟,𝜈−𝑙−𝑚)+ℐ^𝑐(𝜅,𝑟,−𝜈+𝑙+𝑚)+ℐ^𝑐(𝜅,𝑟,−𝜈+𝑙−𝑚), (5.11) where the nonvanishing coefficients 𝑏^𝑒_𝑚 are:𝑏^𝑒₀=123𝑒+𝑒16³+12𝜎cos𝐼𝑛√1−𝑒²327𝑒−𝑒16³,𝑏^𝑒₁=31+8𝑒²+15𝑒64⁴−12𝜎cos𝐼𝑛√1−𝑒²32−4𝑒²−5𝑒32⁴,𝑏^𝑒₂=121𝑒+4𝑒³+12𝜎cos𝐼𝑛√1−𝑒²121𝑒−2𝑒³,𝑏^𝑒₃=18𝑒²+15𝑒128⁴+12𝜎cos𝐼𝑛√1−𝑒²14𝑒²+5𝑒64⁴,𝑏^𝑒₄=1𝑒16³−132𝜎cos𝐼𝑛√1−𝑒²𝑒³,𝑏^𝑒₅=3𝑒128⁴+1128𝜎cos𝐼𝑛√1−𝑒²𝑒⁴. (5.2)

5.3. The Increment in the Inclination

(Δ𝐼)_𝐷=−𝜋𝐶_𝐷(𝐴/𝑚)𝑘_𝑅𝜎sin𝐼√1−𝑒²𝑎⁵𝜇𝑓₀𝑓_𝑥𝑘₀×⁷³_𝑛=1⁶_𝑗=0⁶_𝜈=0⁶_𝑙=0^∞_𝑚=0𝜅=0𝐾_𝑛,𝑗𝐴_𝑗𝐵^𝜅_𝑗×𝑏𝜅!_𝑚𝐼𝑐𝐺^𝑐_{𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈+𝑙+𝑚)+ℐ^𝑐(𝜅,𝑟,𝜈+𝑙−𝑚)+𝑏_𝑚^𝐼𝑐𝐺^𝑐_{𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈−𝑙+𝑚)+ℐ^𝑐(𝜅,𝑟,𝜈−𝑙−𝑚)+𝑏_𝑚^𝐼𝑠𝐺^𝑐_{𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈+𝑙−𝑚)−ℐ^𝑐(𝜅,𝑟,𝜈+𝑙+𝑚)+𝑏_𝑚^𝐼𝑠𝐺^𝑐_{𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈−𝑙−𝑚)−ℐ^𝑐(𝜅,𝑟,𝜈−𝑙+𝑚), (5.13) where the nonvanishing coefficients 𝑏_𝑚^𝐼𝑐 and 𝑏_𝑚^𝐼𝑠 are𝑏₀^𝐼𝑐=34−9𝑒32²+7𝑒256⁴+14+9𝑒32²−7𝑒256⁴𝑏cos2𝜔,_1𝐼𝑐=3𝑒−8𝑒³+3𝑒−8𝑒³,𝑏cos2𝜔_2𝐼𝑐=18𝑒²−5𝑒32⁴+12−38𝑒²−39𝑒256⁴,𝑏cos2𝜔_3𝐼𝑐=18𝑒³(),𝑏1+cos2𝜔_4𝐼𝑐=−1𝑒64⁴−1𝑒16²−5𝑒128⁴,𝑏cos2𝜔_6𝐼𝑐1=−𝑒256⁴𝑏cos2𝜔,_1𝐼𝑠=√1−𝑒²18𝑒³𝑏−𝑒sin2𝜔,_2𝐼𝑠=√1−𝑒²12−18𝑒²−5𝑒256⁴𝑏sin2𝜔,_3𝐼𝑠=18𝑒³√1−𝑒²𝑏sin2𝜔,_4𝐼𝑠=−√1−𝑒²1𝑒16²+1𝑒64⁴𝑏sin2𝜔,_6𝐼𝑠1=−𝑒256⁴√1−𝑒²sin2𝜔. (5.3)

5.4. The Increment in the Ascending Node

(ΔΩ)_𝐷=−𝜋𝐶_𝐷(𝐴/𝑚)𝑘_𝑅𝜎√1−𝑒²𝑎⁵𝜇𝑓₀𝑓_𝑥𝑘₀⁷³_𝑛=1⁶_𝑗=0⁶_𝜈=0⁶_𝑙=0^∞_𝑚=0𝜅=0𝐾_𝑛,𝑗×𝐴_𝑗𝐵^𝜅_𝑗𝑏𝜅!^𝑚Ω𝑐𝐺^𝑐_{𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈+𝑙+𝑚)+ℐ^𝑐(𝜅,𝑟,𝜈+𝑙−𝑚)+𝑏_𝑚^Ω𝑐𝐺^𝑐_{𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈−𝑙+𝑚)+ℐ^𝑐(𝜅,𝑟,𝜈−𝑙−𝑚)+𝑏_𝑚^Ω𝑠𝐺^𝑐_{𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈+𝑙−𝑚)−ℐ^𝑐(𝜅,𝑟,𝜈+𝑙+𝑚)+𝑏_𝑚^Ω𝑠𝐺^𝑐_{𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈−𝑙−𝑚)−ℐ^𝑐(𝜅,𝑟,𝜈−𝑙+𝑚), (5.15) where the nonvanishing coefficients 𝑏_𝑚^Ω𝑐 and 𝑏_𝑚^Ω𝑠 are𝑏₀^Ω𝑐=11𝑒16²−11𝑒64⁴𝑏sin2𝜔,_1Ω𝑐3=−𝑒−8𝑒³𝑏sin2𝜔,_2Ω𝑐=12−38𝑒²−39𝑒256⁴𝑏sin2𝜔,_3Ω𝑐=18𝑒³𝑏sin2𝜔,_4Ω𝑐=−1𝑒16²+3𝑒128⁴𝑏sin2𝜔,_6Ω𝑐1=−𝑒256⁴𝑏sin2𝜔,_1Ω𝑠=√1−𝑒²1−𝑒+8𝑒³𝑏cos2𝜔,_2Ω𝑠=√1−𝑒²12−18𝑒²−5𝑒256⁴𝑏cos2𝜔,_3Ω𝑠=18𝑒³√1−𝑒²𝑏cos2𝜔,_4Ω𝑠=√1−𝑒²1𝑒16²+1𝑒64⁴𝑏cos2𝜔,_6Ω𝑠1=−𝑒256⁴√1−𝑒²cos2𝜔. (5.4)

5.5. The Increment in the Argument of the Perigee

(Δ𝜔)_𝐷=−2𝜋(ΔΩ)cos𝐼−𝜋𝐶_𝐷(𝐴/𝑚)𝑘_𝑅𝜌𝑒𝑓₀𝑓_𝑥𝑘₀⁷³_𝑛=1⁶_𝑗=0⁶_𝜈=0⁵_𝑙=0^∞_𝑚=0𝜅=0𝐾_𝑛,𝑗𝐴_𝑗𝑏_𝑚^𝜔𝑠×𝐵^𝜅_𝑗𝐺𝜅!^{𝑠𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈−𝑙+𝑚)+ℐ^𝑐(𝜅,𝑟,𝜈−𝑙−𝑚)−ℐ^𝑐(𝜅,𝑟,𝜈+𝑙+𝑚)−ℐ^𝑐(+𝐺𝜅,𝑟,𝜈+𝑙−𝑚)_{𝑠𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈−𝑙−𝑚)−ℐ^𝑐(𝜅,𝑟,𝜈−𝑙+𝑚)−ℐ^𝑐(𝜅,𝑟,𝜈+𝑙−𝑚)−ℐ^𝑐(𝜅,𝑟,𝜈+𝑙+𝑚), (5.17) where the nonvanishing coefficients 𝑏_𝑚^𝜔𝑠 are 𝑏₁^𝜔𝑠=𝑎√1−𝑒²12+4𝑒²+3𝑒32⁴−𝜎cos𝐼𝑎⁵𝜇52−4𝑒²+3𝑒32⁴,𝑏₂^𝜔𝑠=𝑎√1−𝑒²1𝑒+4𝑒³+12𝜎cos𝐼𝑎⁵𝜇1𝑒−4𝑒³,𝑏₃^𝜔𝑠=14𝑎√1−𝑒²𝑒²+9𝑒16⁴+14𝜎cos𝐼𝑎⁵𝜇𝑒²−5𝑒16⁴,𝑏₄^𝜔𝑠=18𝑎√1−𝑒²𝑒³−116𝜎cos𝐼𝑎⁵𝜇𝑒³,𝑏₅^𝜔𝑠=3𝑎√641−𝑒²𝑒⁴+164𝜎cos𝐼𝑎⁵𝜇𝑒⁴. (5.18)

5.6. The Increment in the Eccentric Anomaly

(Δ𝐸)_𝐷=^{2𝜋−𝐸𝑏𝑜}_𝐸𝑏𝑜−3Δ𝑎()√2𝑎1−𝑒cos𝐸𝑑𝐸−2𝜋1−𝑒²[Δ𝜔+ΔΩcos𝐼]+𝜋𝐶_𝐷𝐴𝑚𝑘_𝑅𝑎𝑓₀𝑓_𝑥𝑘₀⁷³_𝑛=1⁶_𝑗=0⁶_𝜈=0³_𝑙=0𝑚=1𝐾_𝑛,𝑗𝐴_𝑗𝑏_𝑚^𝐸𝑠×𝐺^𝑠_{𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈−𝑙+𝑚)+ℐ^𝑐(𝜅,𝑟,𝜈−𝑙−𝑚)−ℐ^𝑐(𝜅,𝑟,𝜈+𝑙+𝑚)−ℐ^𝑐(+𝐺𝜅,𝑟,𝜈+𝑙−𝑚)_{𝑠𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈−𝑙−𝑚)−ℐ^𝑐(𝜅,𝑟,𝜈−𝑙+𝑚)−ℐ^𝑐(𝜅,𝑟,𝜈+𝑙−𝑚)−ℐ^𝑐(𝜅,𝑟,𝜈+𝑙+𝑚), (5.19) where the nonvanishing coefficients 𝑏_𝑚^𝐸𝑠 are𝑏₁^𝐸𝑠1=2𝑒−4𝑒³,𝑏₃^𝐸𝑠1=−4𝑒³. (5.20)

6. Drag Effect on the Free Flight Range Ψ

The usual ICBM problem concerns the determination of the free flight range angle. In order to improve the results, we will take into account the perturbations in the orbital elements coming from the atmospheric drag. Following the treatment given by Helali [15] the free flight range equation readsΨcos2=1−𝑄_𝑏𝑜cos²𝜙_𝑏𝑜1+𝑄_𝑏𝑜𝑄_𝑏𝑜−2cos²𝜙_𝑏𝑜, (6.1) where 𝑄_𝑏𝑜=(𝑣/𝑣_𝑐)²=𝑟𝑣²/𝜇 defines the dimensionless parameter, 𝑟 is the magnitude of the position vector of the missile relative to the Earth's center, 𝑣 is the missile speed at any point in its orbit, 𝑣_𝑐 is the corresponding missile circular speed at this point, and 𝜇 is the gravitational parameter. From the orbital mechanics of the two body problem and the symmetry of the free-flight portion (see [41, Figure 6.2.2 and Figure 6.3.3, pages 283–285]) we can express the range equation in terms of the eccentricity 𝑒 and the eccentric anomaly 𝐸 as follows:Ψcos2=(𝑒+cos𝐸)(1−𝑒cos𝐸)⁻¹=(𝑒+cos𝐸)^∞_𝑛=0𝑒^𝑛cos^𝑛𝐸=^∞_𝑛=0𝑒^𝑛+1cos^𝑛𝐸+𝑒^𝑛cos^𝑛+1𝐸=^∞_𝑛=0𝑒^𝑛+1+𝑛|𝑛−1|!𝑒𝑛!^𝑛cos^𝑛𝐸=^∞[]_𝑛=0^𝑛/2_𝑚=0𝐶^𝑛_𝑚2^𝑛−11−𝛿^𝑛_𝑚/2𝑒^𝑛+1+𝑛|𝑛−1|!𝑒^𝑛()𝑛!cos𝑛−2𝑚𝐸. (6.2) Since the range equation is a function the eccentricity and the eccentric anomaly Ψ=Ψ(𝑒,𝐸) and ΔΨ^𝐷_𝑒 and ΔΨ^𝐷_𝐸 denote the variations in the free flight range angle due to the variation of the eccentricity and due to the variation of the eccentric anomaly, respectively, under the effect of the atmospheric drag, therefore, we can write ΔΨ=𝜕Ψ()𝜕𝑒Δ𝑒_𝐷+𝜕Ψ()𝜕𝐸Δ𝐸_𝐷=ΔΨ^𝐷_𝑒+ΔΨ^𝐷_𝐸. (6.3) The involved partial derivatives in (6.3) are given by𝜕Ψ=𝜕𝑒^∞[]_𝑛=0^𝑛/2_𝑚=0𝐶^𝑛_𝑚2^𝑛−11−𝛿^𝑛_𝑚/2(𝑛+1)𝑒^𝑛+𝑛²|𝑛−1|!𝑒^𝑛−1()𝑛!cos𝑛−2𝑚𝐸,𝜕Ψ=𝜕𝐸^∞[]_𝑛=0^𝑛/2_𝑚=0−𝐶^𝑛_𝑚2^𝑛−11−𝛿^𝑛_𝑚/2𝑒^𝑛+1+𝑛|𝑛−1|!𝑒^𝑛)()𝑛!(𝑛−2𝑚sin𝑛−2𝑚𝐸. (6.4)

6.1. Computation of the Quantities ΔΨ^𝐷_𝑒 and ΔΨ^𝐷_𝐸

Substituting from (5.11) and (6.4) into the first term of (6.3) yieldsΔΨ^𝐷_𝑒=^∞[]_𝑛=0^𝑛/2_𝑚=0𝐶^𝑛_𝑚2^𝑛−11−𝛿^𝑛_𝑚/2(𝑛+1)𝑒^𝑛+𝑛²|𝑛−1|!𝑒^𝑛−1𝑛!×𝐶_𝐷𝐴𝑚𝑘_𝑅𝑎²𝑓₀𝑓_𝑥𝑘₀()𝐸×cos𝑛−2𝑚_{73𝑛=16𝑗=06𝜈=05𝑙=0∞𝑚=0𝜅=0}−12𝐾_𝑛,𝑗𝐴_𝑗𝑏^𝑒_𝑚𝐵^𝜅_𝑗×𝐺𝜅!_{𝑐𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈+𝑙+𝑚)+ℐ^𝑐(𝜅,𝑟,𝜈+𝑙−𝑚)+ℐ^𝑐(𝜅,𝑟,−𝜈−𝑙−𝑚)+ℐ^𝑐(+𝐺𝜅,𝑟,−𝜈−𝑙+𝑚)_{𝑐𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈−𝑙+𝑚)+ℐ^𝑐(𝜅,𝑟,𝜈−𝑙−𝑚)+ℐ^𝑐(𝜅,𝑟,−𝜈+𝑙+𝑚)+ℐ^𝑐(𝜅,𝑟,−𝜈+𝑙−𝑚). (6.5)

Substituting again from (5.11) and (6.4) into the second term of (6.3) yieldsΔΨ^𝐷_𝐸=^∞[]_𝑛=0^𝑛/2_𝑚=0−(𝑛−2𝑚)𝐶^𝑛_𝑚2^𝑛−11−𝛿^𝑛_𝑚/2𝑒^𝑛+1+𝑛|𝑛−1|!𝑒^𝑛()𝐸𝑛!×sin𝑛−2𝑚^{2𝜋−𝐸𝑏𝑜}_𝐸𝑏𝑜−32Δ𝑎𝑎()√1−𝑒cos𝐸𝑑𝐸−2𝜋1−𝑒²[Δ𝜔+ΔΩcos𝐼]+𝜋𝐶_𝐷𝐴𝑚𝑘_𝑅×𝑎𝑓₀𝑓_𝑥𝑘₀⁷³_𝑛=1⁶_𝑗=0⁶_𝜈=0³_𝑙=0𝑚=1𝐾_𝑛,𝑗𝐴_𝑗𝑏_𝑚^𝐸𝑠×𝐺^𝑠_{𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈−𝑙+𝑚)+ℐ^𝑐(𝜅,𝑟,𝜈−𝑙−𝑚)−ℐ^𝑐(𝜅,𝑟,𝜈+𝑙+𝑚)−ℐ^𝑐(+𝐺𝜅,𝑟,𝜈+𝑙−𝑚)_{𝑠𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈−𝑙−𝑚)−ℐ^𝑐(𝜅,𝑟,𝜈−𝑙+𝑚)−ℐ^𝑐(𝜅,𝑟,𝜈+𝑙−𝑚)−ℐ^𝑐(𝜅,𝑟,𝜈+𝑙+𝑚). (6.6)

6.2. Errors in Ψ due to (ΔΩ)_𝐷

For the sake of the simplicity, let us take the burnout point on the equator. If the reentry point, for some reason, was displaced by an amount, Δ𝑥, and if all of the orbital elements were kept fixed except the longitude of the ascending node, then this displacement Δ𝑥 could be interpreted as a change in the longitude of the ascending node (ΔΩ)_𝐷. Due to this change, a cross-range error, ΔΨ^𝐷_×Ω, at impact occurs, the value of which is obtained by applying the law of cosines in the spherical trigonometry. Hence, cosΔΨ^𝐷_×Ω=sin²Ψ+cos²Ψcos(ΔΩ)_𝐷. (6.7) Since ΔΨ^𝐷_×Ω and (ΔΩ)_𝐷 are small angles, then we haveΔΨ^𝐷_×Ω≈(ΔΩ)_𝐷cosΨ. (6.8) Substituting (5.15) into (6.8) yieldsΔΨ^𝐷_×Ω≈−2𝜋𝐶_𝐷(𝐴/𝑚)𝑘_𝑅𝜎√1−𝑒²𝑎⁵𝜇𝑓₀𝑓_𝑥𝑘₀×cosΨ_{73𝑛=16𝑗=06𝜈=06𝑙=0∞𝑚=0𝜅=0}12𝐾_𝑛,𝑗𝐴_𝑗𝐵^𝜅_𝑗×𝑏𝜅!_𝑚Ω𝑐𝐺^𝑐_{𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈+𝑙+𝑚)+ℐ^𝑐(𝜅,𝑟,𝜈+𝑙−𝑚)+𝑏_𝑚^Ω𝑐𝐺^𝑐_{𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈−𝑙+𝑚)+ℐ^𝑐(𝜅,𝑟,𝜈−𝑙−𝑚)+𝑏_𝑚^Ω𝑠𝐺^𝑐_{𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈+𝑙−𝑚)−ℐ^𝑐(𝜅,𝑟,𝜈+𝑙+𝑚)+𝑏_𝑚^Ω𝑠𝐺^𝑐_{𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈−𝑙−𝑚)−ℐ^𝑐(𝜅,𝑟,𝜈−𝑙+𝑚). (6.9)

6.3. Errors in Ψ due to (Δ𝐼)_𝐷

Again assume that the burnout point were on the equator and assume also the actual launch azimuth differs from the intended value by an amount, Δ𝛽, and if all of the orbital elements were kept fixed except the orbital inclination, then this displacement Δ𝛽 could be interpreted as a change in the orbital inclination (Δ𝐼)_𝐷. Due to this change a cross-range error, ΔΨ^𝐷_×𝐼, at impact occurs, the value of which is obtained in a similar manner as the previous subsection by applying the law of cosines in the spherical trigonometry. Hence, cosΔΨ^𝐷_×𝐼=cos²Ψ+sin²Ψcos(Δ𝐼)_𝐷. (6.10) Since ΔΨ^𝐷_×𝐼 and (Δ𝐼)_𝐷 are small angles, then we haveΔΨ^𝐷_×𝐼≈(Δ𝐼)_𝐷sinΨ. (6.11) Substituting (5.13) into (6.11) yieldsΔΨ^𝐷_×𝐼≈−2𝜋𝐶_𝐷(𝐴/𝑚)𝑘_𝑅𝜎sin𝐼2√1−𝑒²𝑎⁵𝜇𝑓₀𝑓_𝑥𝑘₀×sinΨ⁷³_𝑛=1⁶_𝑗=0⁶_𝜈=0⁶_𝑙=0^∞_𝑚=0𝜅=0𝐾_𝑛,𝑗𝐴_𝑗𝐵^𝜅_𝑗×𝑏𝜅!_𝑚𝐼𝑐𝐺^𝑐_{𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈+𝑙+𝑚)+ℐ^𝑐(𝜅,𝑟,𝜈+𝑙−𝑚)+𝑏_𝑚^𝐼𝑐𝐺^𝑐_{𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈−𝑙+𝑚)+ℐ^𝑐(𝜅,𝑟,𝜈−𝑙−𝑚)+𝑏_𝑚^𝐼𝑠𝐺^𝑐_{𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈+𝑙−𝑚)−ℐ^𝑐(𝜅,𝑟,𝜈+𝑙+𝑚)+𝑏_𝑚^𝐼𝑠𝐺^𝑐_{𝜈,𝑛,𝑙}ℐ^𝑐(𝜅,𝑟,𝜈−𝑙−𝑚)−ℐ^𝑐(𝜅,𝑟,𝜈−𝑙+𝑚). (6.12)

7. Conclusion

In this work, the equations of motion of a missile are constructed. Using Lagrange's planetary equations in Gauss form, the perturbations due to the air drag in the orbital elements are computed between the eccentric anomalies of the burn out and the re entry points [𝐸_𝑏𝑜,2𝜋−𝐸_𝑏𝑜], respectively. The range equation is expressed in terms of the eccentricity 𝑒 and the eccentric anomaly 𝐸. The different errors in the missile free range due to the drag perturbations in the missile trajectory are obtained. These obtained perturbations can be safely added up linearly to the gravity perturbations on the missile trajectories obtained previously by Abd El-Salam and Abd El-Bar [1] to constitute a more precise theory of the missile trajectories and the errors in the free flight range.

Acknowledgments

The author is deeply indebted to the team work at the deanship of the scientific research—Taibah University for their valuable help and critical guidance and for facilitating many administrative procedures. This research work was financed (supported) by Grant no. (627/1431) from the deanship of the scientific research at Taibah university, Al-Madinah Al- Munawwarah, Saudi Arabia.