# Air Drag Effects on the Missile Trajectories

Authors: F. A. Abd El-Salam
Journal: Journal of Applied Mathematics
 (2011)

---

## 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 [Ebo,2π−Ebo], 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.


---

## Body

## 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 J4. 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 elementsa, the semimajor axis, and e, 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, F→D, per unit mass of the satellite can be represented by(2.1)F→D=-12CD(Am)ρ|V→|V→,
where CD is a dimensionless drag coefficient, V→ is the velocity of the satellite relative to the atmosphere (called the ambient velocity), A 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 V→, m  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)V→=v→-σ→×r→,
where v→ is the orbital velocity of the satellite and r→ its radius vector. If (ı̂,𝚥̂,k̂) are the base vectors of the rectangular geocentric equatorial frame associated with the spherical coordinates (r,λ,δ), then σ→=σ×k̂,σ=7.27205217×10-5 rad./sec. Thus (3.2)V2=v2-2v→⋅(σ→×r→)+(σ→×r→)2=v2kR,
where I is the orbital inclination, δ is the satellite's declination, h is the integral of areas, and kR is given by (3.3)kR=1-2σhcos⁡Iv2+σ2r2cos⁡2δv2.
The third term in kR 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 kR as a constant.


## 4. Drag Force Components

Let(r̂,ŝ,ĥ) 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)F→D=Sr̂+Tŝ+Wĥ,
where(4.2)S=-12CD(Am)kRρvμpesin⁡f,T=-12CD(Am)kRρv(μppr-σrcos⁡I),W=-12CD(Am)kRρvσrsin⁡Icos⁡(f+ω).


### 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 E  as(4.3)ρ=f0fxk0∑n=17∑j=03∑ν=06∑l=06Kn,jAjexp⁡(Bjcos⁡E)×{Gν,n,lccos⁡(ν+l)E+G̃ν,n,lccos⁡(ν-l)E+Gν,n,lssin⁡(ν+l)E+G̃ν,n,lssin⁡(ν-l)E},
where (f0,fx,k0), the main solar and geomagnetic activity parameters affecting the angular density, are given by(4.4)f0=a2+fm,fx=1+a1(Fx-Fb),fm=(Fb-60)160,k0=1+a3(Kp-3)
with Fx solar flux measured on 10.7 cm wavelength for a previous day, Fb  solar flux averaged over three solar rotations (81 days), Kp   geomagnetic index 3 hours before local time t, Hj the reciprocal of the density scale height Hj=j(|j-1|!)2/29(j!)2, R⊕ the equatorial radius of the Earth, ε the Earth's flattening, Aj=exp[Hj(120-a+R⊕)], Bj=Hjae, and ai,Kn,j,pn numerical constants of the model given by Tables 1, 2, and 3.Table 1
Constantsai.
i12345678ai0.0070.28750.047620.04717.07.00.333315.0Table 2
ConstantsKn,j.
nj=0j=1j=2j=312.968×10-157.664×10-091.657×10-103.871×10-1122.815×10-14-4.401×10-093.343×10-109.352×10-113-1.233×10-141.181×10-10-1.478×10-10-1.518×10-124-1.149×10-17-1.597×10-11-6.467×10-12-2.049×10-125-3.901×10-16-2.408×10-10-1.399×10-11-3.059×10-1267.424×10-156.438×10-111.362×10-103.517×10-117-3.416×10-167.447×10-124.542×10-122.080×10-12Table 3
Phasespn.
n34567pn263−263−29.418.091310.0813The required nonvanishing coefficients for the model are(4.5)Gν,n,lc=12(C̃νgn,lc-S̃νgn,ls),G̃ν,n,lc=12(C̃νgn,lc+S̃νgn,ls),Gν,n,ls=12(C̃νgn,ls+S̃νgn,lc),G̃ν,n,ls=12(S̃νgn,lc-C̃νgn,ls),
where the nonvanishing coefficients gn,lc, gn,ls, C̃ν, and S̃ν were computed previously by the author, see Abd El-Salam and Sehnal [26].


## 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 E  as(4.3)ρ=f0fxk0∑n=17∑j=03∑ν=06∑l=06Kn,jAjexp⁡(Bjcos⁡E)×{Gν,n,lccos⁡(ν+l)E+G̃ν,n,lccos⁡(ν-l)E+Gν,n,lssin⁡(ν+l)E+G̃ν,n,lssin⁡(ν-l)E},
where (f0,fx,k0), the main solar and geomagnetic activity parameters affecting the angular density, are given by(4.4)f0=a2+fm,fx=1+a1(Fx-Fb),fm=(Fb-60)160,k0=1+a3(Kp-3)
with Fx solar flux measured on 10.7 cm wavelength for a previous day, Fb  solar flux averaged over three solar rotations (81 days), Kp   geomagnetic index 3 hours before local time t, Hj the reciprocal of the density scale height Hj=j(|j-1|!)2/29(j!)2, R⊕ the equatorial radius of the Earth, ε the Earth's flattening, Aj=exp[Hj(120-a+R⊕)], Bj=Hjae, and ai,Kn,j,pn numerical constants of the model given by Tables 1, 2, and 3.Table 1
Constantsai.
i12345678ai0.0070.28750.047620.04717.07.00.333315.0Table 2
ConstantsKn,j.
nj=0j=1j=2j=312.968×10-157.664×10-091.657×10-103.871×10-1122.815×10-14-4.401×10-093.343×10-109.352×10-113-1.233×10-141.181×10-10-1.478×10-10-1.518×10-124-1.149×10-17-1.597×10-11-6.467×10-12-2.049×10-125-3.901×10-16-2.408×10-10-1.399×10-11-3.059×10-1267.424×10-156.438×10-111.362×10-103.517×10-117-3.416×10-167.447×10-124.542×10-122.080×10-12Table 3
Phasespn.
n34567pn263−263−29.418.091310.0813The required nonvanishing coefficients for the model are(4.5)Gν,n,lc=12(C̃νgn,lc-S̃νgn,ls),G̃ν,n,lc=12(C̃νgn,lc+S̃νgn,ls),Gν,n,ls=12(C̃νgn,ls+S̃νgn,lc),G̃ν,n,ls=12(S̃νgn,lc-C̃νgn,ls),
where the nonvanishing coefficients gn,lc, gn,ls, C̃ν, and S̃ν 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(a,e,I,Ω,ω,E). Gauss modified this system so as to be applicable for all kinds of force especially nonconservative force, see Roy [40],(5)dadt=2n1-e2(Sesin⁡f+Tpr),dedt=1-e2na(Ssin⁡f+T(cos⁡f+cos⁡E)),dIdt=rcos⁡(f+ω)na21-e2W,dΩdt=rsin⁡(f+ω)na2sin⁡I1-e2W,dωdt=-cos⁡IdΩdt+1-e2nae(-Scos⁡f+T(1+rp)sin⁡f),dMdt=n-1-e2(dωdt+cos⁡IdΩdt)-S2rna2
with p=a(1-e2), n2a3=μ, where n,p,μ are the mean motion in the orbit, the semi-latus rectum of the orbit, and the gravitational parameter, respectively, and (S,T,W) 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(S,T,W);(2)
usingdt=((1-ecosE)/n)dE, to replace the explicit time with the eccentric anomaly E as an independent variable;(3)
using the binomial theorem and neglecting orders of𝒪(e5) and 𝒪(σe4);(4)
since the series expansion ofexp(x)for allx∈ℝ, the real numbers, is convergent, then we replace exp(BjcosE) appeared in the density expression with its equivalent series expansion;(5)
let us definecos(α±β)E=cos(α+β)E+cos(α-β)E.


### 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(5.2)dadE=-CD(Am)kRa2f0fxk0∑n=17∑j=03∑ν=06∑l=06∑m=04∑κ=0∞12Kn,jAjbma×Bjκcos⁡κEκ!{Gν,n,lccos⁡(ν+l±m)E+G̃ν,n,lccos⁡(ν-l±m)E+Gν,n,lssin⁡(ν+l±m)E+G̃ν,n,lssin⁡(ν-l±m)E},
where the nonvanishing coefficients are(5.3)b0a=(1+34e2+2164e4)-σcos⁡In(1-e2)(1-14e2-364e4),b1a=(2e+34e3),b2a=(34e2+716e4)+σcos⁡In(1-e2)(14e2+116e4),b3a=14e3,b4a=764e4+164σcos⁡In(1-e2)e4
using the trigonometric identity(5.4)cos⁡κE=[δ[κ/2]κ/212κC[κ/2]κ+12κ-1∑r=0[(κ-1)/2]Crκcos⁡(κ-2r)E],
where δ[κ/2]κ/2 is the usual Kronecker Delta symbol and [κ/2] is the greatest integer function, κ∈ℤ+, and Crκ=κ!/r!(κ-r)! is the binomial coefficients. Now (5.4) can be written in the form(5.5)dadE=-KDa2f0fxk0∑n=17∑j=03∑ν=06∑l=06∑m=04∑κ=0∞Kn,jAjbmaBjκ2κ+1κ!×{δ[κ/2]κ/2C[κ/2]κ[Gν,n,lccos⁡(ν+l±m)E+G̃ν,n,lccos⁡(ν-l±m)E+Gν,n,lssin⁡(ν+l±m)E+G̃ν,n,lssin⁡(ν-l±m)E]+∑r=0[(κ-1)/2]Crκ[Gν,n,lccos⁡(κ-2r±(ν+l±m))E+G̃ν,n,lccos⁡(κ-2r±(ν-l±m))E+Gν,n,lssin⁡((ν+l±m)±(κ-2r))E+G̃ν,n,lssin⁡((ν-l±m)±(κ-2r))E]}.
Now (5.5) can be integrated between the eccentric anomalies of the burn out and the reentry points [Ebo,2π-Ebo], respectively, then we have integrals of the form(5.6)Ic(κ,r,μ)=∫Ebo2π-Ebo∑κ=0∞Bjκκ![δ[κ/2]κ/212κC[κ/2]κcos⁡μE+12κ∑r=0[(κ-1)/2]Crκ{cos⁡(κ-2r±μ)E}]dE,
from which we can get(5.7)Ic(κ,r,μ)=2π∑κ=0∞((1/2)Bj)κκ!{δ[κ/2]κ/2C[κ/2]κ[1-sin⁡(μEbo)μ]μ≠0+∑r=0[(κ-1)/2]Crκ[1-sin⁡((κ-2r±μ)Ebo)(κ-2r±μ)]|(κ-2r±μ)≠0}.
Similarly, we can evaluate(5.8)Is(κ,r,μ)=∫Ebo,2π-Ebo∑κ=0∞Bjκ2κκ![δ[κ/2]κ/2C[κ/2]κ+2∑r=0[(κ-1)/2]Crκcos⁡(κ-2r)E]sin⁡μE=∫Ebo,2π-Ebo∑κ=0∞Bjκ2κκ![δ[κ/2]κ/2C[κ/2]κsin⁡μE+∑r=0[(κ-1)/2]Crκsin⁡(κ-2r±μ)E],
which yields(5.9)Is(κ,r,μ)=∑κ=0∞Bjκ2κκ!{-δ[κ/2]κ/2C[κ/2]κcos⁡(μE)μ|Ebo,,μ≠02π-Ebo-∑r=0[(κ-1)/2]Crκ[cos⁡((κ-2r±μ)E)(κ-2r±μ)]|Ebo,,(κ-2r±μ)≠02π-Ebo}=0,
where terms involving the integrals of the type ℐs(κ,r,μ) 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 ℐc(κ,r,μ) from (5.7) into (5.5) yields the solution as(5.10)(Δa)D=-KDa2f0fxk0∑n=17∑j=03∑ν=06∑l=06∑m=04∑κ=0∞12Kn,jAjbmaBjκκ!×{Gν,n,lc[Ic(κ,r,-ν-l+m)Ic(κ,r,ν+l+m)+Ic(κ,r,ν+l-m)+Ic(κ,r,-ν-l-m)+Ic(κ,r,-ν-l+m)]+G̃ν,n,lc[G̃ν,n,lcIc(κ,r,ν-l+m)+Ic(κ,r,ν-l-m)+G̃ν,n,lcIc(κ,r,-ν+l+m)+Ic(κ,r,-ν+l-m)]Gν,n,lc}.
Using the similar procedure to integrating the rest of (5) between the eccentric anomalies of the burn out and the re entry points [Ebo,2π-Ebo], respectively, we obtained the following explicit expressions for the other orbital elements as follows.


### 5.2. The Increment in the Eccentricity

(5.11)(Δe)D=-CD(Am)kRa2f0fxk0∑n=17∑j=03∑ν=06∑l=06∑m=05∑κ=0∞12Kn,jAjbme×Bjκκ!{Gν,n,lc[IcIc(κ,r,ν+l+m)+Ic(κ,r,ν+l-m)+Ic(κ,r,-ν-l-m)+Ic(κ,r,-ν-l+m)]+G̃ν,n,lc[Ic(κ,r,ν-l+m)+Ic(κ,r,ν-l-m)]+Ic(κ,r,-ν+l+m)+Ic(κ,r,-ν+l-m)]Gν,n,lc},
where the nonvanishing coefficients bme are:(5.2)b0e=(12e+316e3)+12σcos⁡In1-e2(32e-716e3),b1e=(1+38e2+1564e4)-12σcos⁡In1-e2(2-34e2-532e4),b2e=(12e+14e3)+12σcos⁡In1-e2(12e-12e3),b3e=(18e2+15128e4)+12σcos⁡In1-e2(14e2+564e4),b4e=116e3-132σcos⁡In1-e2e3,b5e=3128e4+1128σcos⁡In1-e2e4.


### 5.3. The Increment in the Inclination

(5.13)(ΔI)D=-πCD(A/m)kRσsin⁡I1-e2a5μf0fxk0×∑n=17∑j=03∑ν=06∑l=06∑m=06∑κ=0∞Kn,jAjBjκκ!×{bmIcGν,n,lc[Ic(κ,r,ν+l+m)+Ic(κ,r,ν+l-m)]+bmIcG̃ν,n,lc[Ic(κ,r,ν-l+m)+Ic(κ,r,ν-l-m)]+bmIsGν,n,lc[Ic(κ,r,ν+l-m)-Ic(κ,r,ν+l+m)]+bmIsG̃ν,n,lc[Ic(κ,r,ν-l-m)-Ic(κ,r,ν-l+m)]},
where the nonvanishing coefficients bmIc and bmIs are(5.3)b0Ic=(34-932e2+7256e4)+(14+932e2-7256e4)cos⁡2ω,b1Ic=((e-38e3)+(e-38e3)cos⁡2ω),b2Ic=((18e2-532e4)+(12-38e2-39256e4)cos⁡2ω),b3Ic=18e3(1+cos⁡2ω),b4Ic=(-164e4-(116e2-5128e4)cos⁡2ω),b6Ic=-1256e4cos⁡2ω,b1Is=1-e2(18e3-e)sin⁡2ω,b2Is=1-e2(12-18e2-5256e4)sin⁡2ω,b3Is=18e31-e2sin⁡2ω,b4Is=-1-e2(116e2+164e4)sin⁡2ω,b6Is=-1256e41-e2sin⁡2ω.


### 5.4. The Increment in the Ascending Node

(5.15)(ΔΩ)D=-πCD(A/m)kRσ1-e2a5μf0fxk0∑n=17∑j=03∑ν=06∑l=06∑m=06∑κ=0∞Kn,j×AjBjκκ!{bmΩcG̃ν,n,lcbmΩcGν,n,lc[Ic(κ,r,ν+l+m)+Ic(κ,r,ν+l-m)]+bmΩcG̃ν,n,lc[Ic(κ,r,ν-l+m)+Ic(κ,r,ν-l-m)]+bmΩsGν,n,lc[Ic(κ,r,ν+l-m)-Ic(κ,r,ν+l+m)]+bmΩsG̃ν,n,lc[Ic(κ,r,ν-l-m)-Ic(κ,r,ν-l+m)]},
where the nonvanishing coefficients bmΩc and bmΩs are(5.4)b0Ωc=(1116e2-1164e4)sin⁡2ω,b1Ωc=-(e-38e3)sin⁡2ω,b2Ωc=(12-38e2-39256e4)sin⁡2ω,b3Ωc=18e3sin⁡2ω,b4Ωc=(-116e2+3128e4)sin⁡2ω,b6Ωc=-1256e4sin⁡2ω,b1Ωs=1-e2(-e+18e3)cos⁡2ω,b2Ωs=1-e2(12-18e2-5256e4)cos⁡2ω,b3Ωs=18e31-e2cos⁡2ω,b4Ωs=1-e2(116e2+164e4)cos⁡2ω,b6Ωs=-1256e41-e2cos⁡2ω.


### 5.5. The Increment in the Argument of the Perigee

(5.17)(Δω)D=-2π(ΔΩ)cos⁡I-πCD(A/m)kRρef0fxk0∑n=17∑j=03∑ν=06∑l=06∑m=05∑κ=0∞Kn,jAjbmωs×Bjκκ!{Gν,n,ls[Ic(κ,r,ν-l+m)+Ic(κ,r,ν-l-m)-Ic(κ,r,ν+l+m)-Ic(κ,r,ν+l-m)]+G̃ν,n,ls[Ic(κ,r,ν-l-m)-Ic(κ,r,ν-l+m)-Ic(κ,r,ν+l-m)-Ic(κ,r,ν+l+m)]Gν,n,ls},
where the nonvanishing coefficients bmωs are (5.18)b1ωs=a1-e2(2+14e2+332e4)-σcos⁡Ia5μ(2-54e2+332e4),b2ωs=a1-e2(e+14e3)+12σcos⁡Ia5μ(e-14e3),b3ωs=14a1-e2(e2+916e4)+14σcos⁡Ia5μ(e2-516e4),b4ωs=18a1-e2e3-116σcos⁡Ia5μe3,b5ωs=364a1-e2e4+164σcos⁡Ia5μe4.


### 5.6. The Increment in the Eccentric Anomaly

(5.19)(ΔE)D=∫Ebo2π-Ebo-3Δa2a(1-ecos⁡E)dE-2π1-e2[Δω+ΔΩcos⁡I]+πCD(Am)kRaf0fxk0∑n=17∑j=03∑ν=06∑l=06∑m=13Kn,jAjbmEs×{Gν,n,ls[Ic(κ,r,ν-l+m)+Ic(κ,r,ν-l-m)-Ic(κ,r,ν+l+m)-Ic(κ,r,ν+l-m)]+G̃ν,n,ls[Ic(κ,r,ν-l-m)-Ic(κ,r,ν-l+m)-Ic(κ,r,ν+l-m)-Ic(κ,r,ν+l+m)]},
where the nonvanishing coefficients bmEs are(5.20)b1Es=2e-14e3,b3Es=-14e3.


## 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(5.2)dadE=-CD(Am)kRa2f0fxk0∑n=17∑j=03∑ν=06∑l=06∑m=04∑κ=0∞12Kn,jAjbma×Bjκcos⁡κEκ!{Gν,n,lccos⁡(ν+l±m)E+G̃ν,n,lccos⁡(ν-l±m)E+Gν,n,lssin⁡(ν+l±m)E+G̃ν,n,lssin⁡(ν-l±m)E},
where the nonvanishing coefficients are(5.3)b0a=(1+34e2+2164e4)-σcos⁡In(1-e2)(1-14e2-364e4),b1a=(2e+34e3),b2a=(34e2+716e4)+σcos⁡In(1-e2)(14e2+116e4),b3a=14e3,b4a=764e4+164σcos⁡In(1-e2)e4
using the trigonometric identity(5.4)cos⁡κE=[δ[κ/2]κ/212κC[κ/2]κ+12κ-1∑r=0[(κ-1)/2]Crκcos⁡(κ-2r)E],
where δ[κ/2]κ/2 is the usual Kronecker Delta symbol and [κ/2] is the greatest integer function, κ∈ℤ+, and Crκ=κ!/r!(κ-r)! is the binomial coefficients. Now (5.4) can be written in the form(5.5)dadE=-KDa2f0fxk0∑n=17∑j=03∑ν=06∑l=06∑m=04∑κ=0∞Kn,jAjbmaBjκ2κ+1κ!×{δ[κ/2]κ/2C[κ/2]κ[Gν,n,lccos⁡(ν+l±m)E+G̃ν,n,lccos⁡(ν-l±m)E+Gν,n,lssin⁡(ν+l±m)E+G̃ν,n,lssin⁡(ν-l±m)E]+∑r=0[(κ-1)/2]Crκ[Gν,n,lccos⁡(κ-2r±(ν+l±m))E+G̃ν,n,lccos⁡(κ-2r±(ν-l±m))E+Gν,n,lssin⁡((ν+l±m)±(κ-2r))E+G̃ν,n,lssin⁡((ν-l±m)±(κ-2r))E]}.
Now (5.5) can be integrated between the eccentric anomalies of the burn out and the reentry points [Ebo,2π-Ebo], respectively, then we have integrals of the form(5.6)Ic(κ,r,μ)=∫Ebo2π-Ebo∑κ=0∞Bjκκ![δ[κ/2]κ/212κC[κ/2]κcos⁡μE+12κ∑r=0[(κ-1)/2]Crκ{cos⁡(κ-2r±μ)E}]dE,
from which we can get(5.7)Ic(κ,r,μ)=2π∑κ=0∞((1/2)Bj)κκ!{δ[κ/2]κ/2C[κ/2]κ[1-sin⁡(μEbo)μ]μ≠0+∑r=0[(κ-1)/2]Crκ[1-sin⁡((κ-2r±μ)Ebo)(κ-2r±μ)]|(κ-2r±μ)≠0}.
Similarly, we can evaluate(5.8)Is(κ,r,μ)=∫Ebo,2π-Ebo∑κ=0∞Bjκ2κκ![δ[κ/2]κ/2C[κ/2]κ+2∑r=0[(κ-1)/2]Crκcos⁡(κ-2r)E]sin⁡μE=∫Ebo,2π-Ebo∑κ=0∞Bjκ2κκ![δ[κ/2]κ/2C[κ/2]κsin⁡μE+∑r=0[(κ-1)/2]Crκsin⁡(κ-2r±μ)E],
which yields(5.9)Is(κ,r,μ)=∑κ=0∞Bjκ2κκ!{-δ[κ/2]κ/2C[κ/2]κcos⁡(μE)μ|Ebo,,μ≠02π-Ebo-∑r=0[(κ-1)/2]Crκ[cos⁡((κ-2r±μ)E)(κ-2r±μ)]|Ebo,,(κ-2r±μ)≠02π-Ebo}=0,
where terms involving the integrals of the type ℐs(κ,r,μ) 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 ℐc(κ,r,μ) from (5.7) into (5.5) yields the solution as(5.10)(Δa)D=-KDa2f0fxk0∑n=17∑j=03∑ν=06∑l=06∑m=04∑κ=0∞12Kn,jAjbmaBjκκ!×{Gν,n,lc[Ic(κ,r,-ν-l+m)Ic(κ,r,ν+l+m)+Ic(κ,r,ν+l-m)+Ic(κ,r,-ν-l-m)+Ic(κ,r,-ν-l+m)]+G̃ν,n,lc[G̃ν,n,lcIc(κ,r,ν-l+m)+Ic(κ,r,ν-l-m)+G̃ν,n,lcIc(κ,r,-ν+l+m)+Ic(κ,r,-ν+l-m)]Gν,n,lc}.
Using the similar procedure to integrating the rest of (5) between the eccentric anomalies of the burn out and the re entry points [Ebo,2π-Ebo], respectively, we obtained the following explicit expressions for the other orbital elements as follows.


## 5.2. The Increment in the Eccentricity

(5.11)(Δe)D=-CD(Am)kRa2f0fxk0∑n=17∑j=03∑ν=06∑l=06∑m=05∑κ=0∞12Kn,jAjbme×Bjκκ!{Gν,n,lc[IcIc(κ,r,ν+l+m)+Ic(κ,r,ν+l-m)+Ic(κ,r,-ν-l-m)+Ic(κ,r,-ν-l+m)]+G̃ν,n,lc[Ic(κ,r,ν-l+m)+Ic(κ,r,ν-l-m)]+Ic(κ,r,-ν+l+m)+Ic(κ,r,-ν+l-m)]Gν,n,lc},
where the nonvanishing coefficients bme are:(5.2)b0e=(12e+316e3)+12σcos⁡In1-e2(32e-716e3),b1e=(1+38e2+1564e4)-12σcos⁡In1-e2(2-34e2-532e4),b2e=(12e+14e3)+12σcos⁡In1-e2(12e-12e3),b3e=(18e2+15128e4)+12σcos⁡In1-e2(14e2+564e4),b4e=116e3-132σcos⁡In1-e2e3,b5e=3128e4+1128σcos⁡In1-e2e4.


## 5.3. The Increment in the Inclination

(5.13)(ΔI)D=-πCD(A/m)kRσsin⁡I1-e2a5μf0fxk0×∑n=17∑j=03∑ν=06∑l=06∑m=06∑κ=0∞Kn,jAjBjκκ!×{bmIcGν,n,lc[Ic(κ,r,ν+l+m)+Ic(κ,r,ν+l-m)]+bmIcG̃ν,n,lc[Ic(κ,r,ν-l+m)+Ic(κ,r,ν-l-m)]+bmIsGν,n,lc[Ic(κ,r,ν+l-m)-Ic(κ,r,ν+l+m)]+bmIsG̃ν,n,lc[Ic(κ,r,ν-l-m)-Ic(κ,r,ν-l+m)]},
where the nonvanishing coefficients bmIc and bmIs are(5.3)b0Ic=(34-932e2+7256e4)+(14+932e2-7256e4)cos⁡2ω,b1Ic=((e-38e3)+(e-38e3)cos⁡2ω),b2Ic=((18e2-532e4)+(12-38e2-39256e4)cos⁡2ω),b3Ic=18e3(1+cos⁡2ω),b4Ic=(-164e4-(116e2-5128e4)cos⁡2ω),b6Ic=-1256e4cos⁡2ω,b1Is=1-e2(18e3-e)sin⁡2ω,b2Is=1-e2(12-18e2-5256e4)sin⁡2ω,b3Is=18e31-e2sin⁡2ω,b4Is=-1-e2(116e2+164e4)sin⁡2ω,b6Is=-1256e41-e2sin⁡2ω.


## 5.4. The Increment in the Ascending Node

(5.15)(ΔΩ)D=-πCD(A/m)kRσ1-e2a5μf0fxk0∑n=17∑j=03∑ν=06∑l=06∑m=06∑κ=0∞Kn,j×AjBjκκ!{bmΩcG̃ν,n,lcbmΩcGν,n,lc[Ic(κ,r,ν+l+m)+Ic(κ,r,ν+l-m)]+bmΩcG̃ν,n,lc[Ic(κ,r,ν-l+m)+Ic(κ,r,ν-l-m)]+bmΩsGν,n,lc[Ic(κ,r,ν+l-m)-Ic(κ,r,ν+l+m)]+bmΩsG̃ν,n,lc[Ic(κ,r,ν-l-m)-Ic(κ,r,ν-l+m)]},
where the nonvanishing coefficients bmΩc and bmΩs are(5.4)b0Ωc=(1116e2-1164e4)sin⁡2ω,b1Ωc=-(e-38e3)sin⁡2ω,b2Ωc=(12-38e2-39256e4)sin⁡2ω,b3Ωc=18e3sin⁡2ω,b4Ωc=(-116e2+3128e4)sin⁡2ω,b6Ωc=-1256e4sin⁡2ω,b1Ωs=1-e2(-e+18e3)cos⁡2ω,b2Ωs=1-e2(12-18e2-5256e4)cos⁡2ω,b3Ωs=18e31-e2cos⁡2ω,b4Ωs=1-e2(116e2+164e4)cos⁡2ω,b6Ωs=-1256e41-e2cos⁡2ω.


## 5.5. The Increment in the Argument of the Perigee

(5.17)(Δω)D=-2π(ΔΩ)cos⁡I-πCD(A/m)kRρef0fxk0∑n=17∑j=03∑ν=06∑l=06∑m=05∑κ=0∞Kn,jAjbmωs×Bjκκ!{Gν,n,ls[Ic(κ,r,ν-l+m)+Ic(κ,r,ν-l-m)-Ic(κ,r,ν+l+m)-Ic(κ,r,ν+l-m)]+G̃ν,n,ls[Ic(κ,r,ν-l-m)-Ic(κ,r,ν-l+m)-Ic(κ,r,ν+l-m)-Ic(κ,r,ν+l+m)]Gν,n,ls},
where the nonvanishing coefficients bmωs are (5.18)b1ωs=a1-e2(2+14e2+332e4)-σcos⁡Ia5μ(2-54e2+332e4),b2ωs=a1-e2(e+14e3)+12σcos⁡Ia5μ(e-14e3),b3ωs=14a1-e2(e2+916e4)+14σcos⁡Ia5μ(e2-516e4),b4ωs=18a1-e2e3-116σcos⁡Ia5μe3,b5ωs=364a1-e2e4+164σcos⁡Ia5μe4.


## 5.6. The Increment in the Eccentric Anomaly

(5.19)(ΔE)D=∫Ebo2π-Ebo-3Δa2a(1-ecos⁡E)dE-2π1-e2[Δω+ΔΩcos⁡I]+πCD(Am)kRaf0fxk0∑n=17∑j=03∑ν=06∑l=06∑m=13Kn,jAjbmEs×{Gν,n,ls[Ic(κ,r,ν-l+m)+Ic(κ,r,ν-l-m)-Ic(κ,r,ν+l+m)-Ic(κ,r,ν+l-m)]+G̃ν,n,ls[Ic(κ,r,ν-l-m)-Ic(κ,r,ν-l+m)-Ic(κ,r,ν+l-m)-Ic(κ,r,ν+l+m)]},
where the nonvanishing coefficients bmEs are(5.20)b1Es=2e-14e3,b3Es=-14e3.


## 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(6.1)cos⁡Ψ2=1-Qbocos⁡2ϕbo1+Qbo(Qbo-2)cos⁡2ϕbo,
where Qbo=(v/vc)2=rv2/μ  defines the dimensionless parameter, r  is the magnitude of the position vector of the missile relative to the Earth's center, v is the missile speed at any point in its orbit, vc 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 e and the eccentric anomaly E as follows:(6.2)cos⁡Ψ2=(e+cos⁡E)(1-ecos⁡E)-1=(e+cos⁡E)∑n=0∞encos⁡nE=∑n=0∞{en+1cos⁡nE+encos⁡n+1E}=∑n=0∞{en+1+n|n-1|!n!en}cos⁡nE=∑n=0∞∑m=0[n/2]Cmn2n-1(1-δm/2n){en+1+n|n-1|!enn!}cos⁡(n-2m)E.
Since the range equation is a function the eccentricity and the eccentric anomaly Ψ=Ψ(e,E) and ΔΨeD and ΔΨED 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)ΔΨ=∂Ψ∂e(Δe)D+∂Ψ∂E(ΔE)D=ΔΨeD+ΔΨED.
The involved partial derivatives in (6.3) are given by(6.4)∂Ψ∂e=∑n=0∞∑m=0[n/2]Cmn2n-1(1-δm/2n){(n+1)en+n2|n-1|!en-1n!}cos⁡(n-2m)E,∂Ψ∂E=∑n=0∞∑m=0[n/2]-Cmn2n-1(1-δm/2n){en+1+n|n-1|!enn!}(n-2m)sin⁡(n-2m)E.


### 6.1. Computation of the QuantitiesΔΨeD and ΔΨED

Substituting from (5.11) and (6.4) into the first term of (6.3) yields(6.5)ΔΨeD=∑n=0∞∑m=0[n/2]Cmn2n-1(1-δm/2n){(n+1)en+n2|n-1|!en-1n!}×CD(Am)kRa2f0fxk0cos⁡(n-2m)E×∑n=17∑j=03∑ν=06∑l=06∑m=05∑κ=0∞-12Kn,jAjbmeBjκκ!×{Gν,n,lc[Ic(κ,r,ν+l+m)+Ic(κ,r,ν+l-m)+Ic(κ,r,-ν-l-m)+Ic(κ,r,-ν-l+m)]+G̃ν,n,lc[Ic(κ,r,ν-l+m)+Ic(κ,r,ν-l-m)+Ic(κ,r,-ν+l+m)+Ic(κ,r,-ν+l-m)]Gν,n,lc}.Substituting again from (5.11) and (6.4) into the second term of (6.3) yields(6.6)ΔΨED=∑n=0∞∑m=0[n/2]-(n-2m)Cmn2n-1(1-δm/2n){en+1+n|n-1|!enn!}×sin⁡(n-2m)E{∫Ebo2π-Ebo-32Δaa(1-ecos⁡E)dE-2π1-e2[Δω+ΔΩcos⁡I]+πCD(Am)kR×af0fxk0∑n=17∑j=03∑ν=06∑l=06∑m=13Kn,jAjbmEs×{Gν,n,ls[Ic(κ,r,ν-l+m)+Ic(κ,r,ν-l-m)-Ic(κ,r,ν+l+m)-Ic(κ,r,ν+l-m)]+G̃ν,n,ls[Ic(κ,r,ν-l-m)-Ic(κ,r,ν-l+m)-Ic(κ,r,ν+l-m)-Ic(κ,r,ν+l+m)]Gν,n,ls}∫Ebo2π-Ebo-32Δaa(1-ecos⁡E)dE}.


### 6.2. Errors inΨ due to (ΔΩ)D

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,Δx, and if all of the orbital elements were kept fixed except the longitude of the ascending node, then this displacement Δx could be interpreted as a change in the longitude of the ascending node (ΔΩ)D. Due to this change, a cross-range error, ΔΨ×ΩD, at impact occurs, the value of which is obtained by applying the law of cosines in the spherical trigonometry. Hence, (6.7)cos⁡ΔΨ×ΩD=sin⁡2Ψ+cos⁡2Ψcos⁡(ΔΩ)D.
Since ΔΨ×ΩD and (ΔΩ)D are small angles, then we have(6.8)ΔΨ×ΩD≈(ΔΩ)Dcos⁡Ψ.
Substituting (5.15) into (6.8) yields(6.9)ΔΨ×ΩD≈-2πCD(A/m)kRσ1-e2a5μf0fxk0cos⁡Ψ×∑n=17∑j=03∑ν=06∑l=06∑m=06∑κ=0∞12Kn,jAjBjκκ!×{bmΩcG̃ν,n,lcbmΩcGν,n,lc[Ic(κ,r,ν+l+m)+Ic(κ,r,ν+l-m)]+bmΩcG̃ν,n,lc[Ic(κ,r,ν-l+m)+Ic(κ,r,ν-l-m)]+bmΩsG̃ν,n,lc[Ic(κ,r,ν+l-m)-Ic(κ,r,ν+l+m)]+bmΩsG̃ν,n,lc[Ic(κ,r,ν-l-m)-Ic(κ,r,ν-l+m)]}.


### 6.3. Errors inΨ due to (ΔI)D

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 (ΔI)D. Due to this change a cross-range error, ΔΨ×ID, 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, (6.10)cos⁡ΔΨ×ID=cos⁡2Ψ+sin⁡2Ψcos⁡(ΔI)D.
Since ΔΨ×ID and (ΔI)D are small angles, then we have(6.11)ΔΨ×ID≈(ΔI)Dsin⁡Ψ.
Substituting (5.13) into (6.11) yields(6.12)ΔΨ×ID≈-2πCD(A/m)kRσsin⁡I21-e2a5μf0fxk0×sin⁡Ψ∑n=17∑j=03∑ν=06∑l=06∑m=06∑κ=0∞Kn,jAjBjκκ!×{bmIcGν,n,lc[Ic(κ,r,ν+l+m)+Ic(κ,r,ν+l-m)]+bmIcG̃ν,n,lc[Ic(κ,r,ν-l+m)+Ic(κ,r,ν-l-m)]+bmIsGν,n,lc[Ic(κ,r,ν+l-m)-Ic(κ,r,ν+l+m)]+bmIsG̃ν,n,lc[Ic(κ,r,ν-l-m)-Ic(κ,r,ν-l+m)]}.


## 6.1. Computation of the QuantitiesΔΨeD and ΔΨED

Substituting from (5.11) and (6.4) into the first term of (6.3) yields(6.5)ΔΨeD=∑n=0∞∑m=0[n/2]Cmn2n-1(1-δm/2n){(n+1)en+n2|n-1|!en-1n!}×CD(Am)kRa2f0fxk0cos⁡(n-2m)E×∑n=17∑j=03∑ν=06∑l=06∑m=05∑κ=0∞-12Kn,jAjbmeBjκκ!×{Gν,n,lc[Ic(κ,r,ν+l+m)+Ic(κ,r,ν+l-m)+Ic(κ,r,-ν-l-m)+Ic(κ,r,-ν-l+m)]+G̃ν,n,lc[Ic(κ,r,ν-l+m)+Ic(κ,r,ν-l-m)+Ic(κ,r,-ν+l+m)+Ic(κ,r,-ν+l-m)]Gν,n,lc}.Substituting again from (5.11) and (6.4) into the second term of (6.3) yields(6.6)ΔΨED=∑n=0∞∑m=0[n/2]-(n-2m)Cmn2n-1(1-δm/2n){en+1+n|n-1|!enn!}×sin⁡(n-2m)E{∫Ebo2π-Ebo-32Δaa(1-ecos⁡E)dE-2π1-e2[Δω+ΔΩcos⁡I]+πCD(Am)kR×af0fxk0∑n=17∑j=03∑ν=06∑l=06∑m=13Kn,jAjbmEs×{Gν,n,ls[Ic(κ,r,ν-l+m)+Ic(κ,r,ν-l-m)-Ic(κ,r,ν+l+m)-Ic(κ,r,ν+l-m)]+G̃ν,n,ls[Ic(κ,r,ν-l-m)-Ic(κ,r,ν-l+m)-Ic(κ,r,ν+l-m)-Ic(κ,r,ν+l+m)]Gν,n,ls}∫Ebo2π-Ebo-32Δaa(1-ecos⁡E)dE}.


## 6.2. Errors inΨ due to (ΔΩ)D

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,Δx, and if all of the orbital elements were kept fixed except the longitude of the ascending node, then this displacement Δx could be interpreted as a change in the longitude of the ascending node (ΔΩ)D. Due to this change, a cross-range error, ΔΨ×ΩD, at impact occurs, the value of which is obtained by applying the law of cosines in the spherical trigonometry. Hence, (6.7)cos⁡ΔΨ×ΩD=sin⁡2Ψ+cos⁡2Ψcos⁡(ΔΩ)D.
Since ΔΨ×ΩD and (ΔΩ)D are small angles, then we have(6.8)ΔΨ×ΩD≈(ΔΩ)Dcos⁡Ψ.
Substituting (5.15) into (6.8) yields(6.9)ΔΨ×ΩD≈-2πCD(A/m)kRσ1-e2a5μf0fxk0cos⁡Ψ×∑n=17∑j=03∑ν=06∑l=06∑m=06∑κ=0∞12Kn,jAjBjκκ!×{bmΩcG̃ν,n,lcbmΩcGν,n,lc[Ic(κ,r,ν+l+m)+Ic(κ,r,ν+l-m)]+bmΩcG̃ν,n,lc[Ic(κ,r,ν-l+m)+Ic(κ,r,ν-l-m)]+bmΩsG̃ν,n,lc[Ic(κ,r,ν+l-m)-Ic(κ,r,ν+l+m)]+bmΩsG̃ν,n,lc[Ic(κ,r,ν-l-m)-Ic(κ,r,ν-l+m)]}.


## 6.3. Errors inΨ due to (ΔI)D

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 (ΔI)D. Due to this change a cross-range error, ΔΨ×ID, 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, (6.10)cos⁡ΔΨ×ID=cos⁡2Ψ+sin⁡2Ψcos⁡(ΔI)D.
Since ΔΨ×ID and (ΔI)D are small angles, then we have(6.11)ΔΨ×ID≈(ΔI)Dsin⁡Ψ.
Substituting (5.13) into (6.11) yields(6.12)ΔΨ×ID≈-2πCD(A/m)kRσsin⁡I21-e2a5μf0fxk0×sin⁡Ψ∑n=17∑j=03∑ν=06∑l=06∑m=06∑κ=0∞Kn,jAjBjκκ!×{bmIcGν,n,lc[Ic(κ,r,ν+l+m)+Ic(κ,r,ν+l-m)]+bmIcG̃ν,n,lc[Ic(κ,r,ν-l+m)+Ic(κ,r,ν-l-m)]+bmIsGν,n,lc[Ic(κ,r,ν+l-m)-Ic(κ,r,ν+l+m)]+bmIsG̃ν,n,lc[Ic(κ,r,ν-l-m)-Ic(κ,r,ν-l+m)]}.


## 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[Ebo,2π-Ebo], respectively. The range equation is expressed 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. 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.



---

*Source: 871304-2011-10-11.xml*