Editing Projectile motion (section)

== Trajectory of a projectile with air resistance ==
[[File:Inclinedthrow2.gif|thumb|400px|Trajectories of a mass thrown at an angle of 70°:<br />
{{color box|black}} without [[Drag (physics)|drag]] (a parabole)<br />
{{color box|blue}} with [[Stokes' law|Stokes' drag]]<br />
{{color box|green}} with [[Newtonian drag]]]]
Air resistance creates a force that (for symmetric projectiles) is always directed against the direction of motion in the surrounding medium and has a magnitude that depends on the absolute speed: <math>\mathbf{F_{air}} = -f(v)\cdot\mathbf{\hat v}</math>. The speed-dependence of the friction force is linear (<math>f(v) \propto v</math>) at very low speeds ([[Stokes drag]]) and quadratic (<math>f(v) \propto v^2</math>) at large speeds ([[Drag equation|Newton drag]]).<ref name="ThorntonMarion2007">{{cite book|author1=Stephen T. Thornton|author2=Jerry B. Marion|title=Classical Dynamics of Particles and Systems|url=https://books.google.com/books?id=30rWGAAACAAJ|year=2007|publisher=Brooks/Cole|isbn=978-0-495-55610-7|page=59}}</ref> The transition between these behaviours is determined by the [[Reynolds number]], which depends on object speed and size, density <math display="inline">\rho</math> and [[Viscosity|dynamic viscosity]] <math display="inline">\eta</math> of the medium. For Reynolds numbers below about 1 the dependence is linear, above 1000 ([[Turbulence modeling|turbulent flow]]) it becomes quadratic. In air, which has a [[kinematic viscosity]] <math>\eta/\rho</math> around 0.15&nbsp;cm<sup>2</sup>/s, this means that the drag force becomes quadratic in ''v'' when the product of object speed and diameter is more than about 0.015 m<sup>2</sup>/s, which is typically the case for projectiles.
* Stokes drag: <math>\mathbf{F_{air}} = -k_{\mathrm{Stokes}}\cdot\mathbf{v}\qquad</math> (for <math display="inline">Re\lesssim 1</math>)
* Newton drag: <math>\mathbf{F_{air}} = -k\,|\mathbf{v}|\cdot\mathbf{v}\qquad</math> (for <math display="inline">Re\gtrsim 1000</math>)

[[Image:Free body diagram gravity air resistance.svg|right|thumb|320px|Free body diagram of a body on which only gravity and air resistance act.]]
The [[free body diagram]] on the right is for a projectile that experiences air resistance and the effects of gravity. Here, air resistance is assumed to be in the direction opposite of the projectile's velocity: <math>\mathbf{F_{\mathrm{air}}} = -f(v)\cdot\mathbf{\hat v}</math>

=== Trajectory of a projectile with Stokes drag ===
Stokes drag, where <math>\mathbf{F_{air}} \propto \mathbf{v}</math>, only applies at very low speed in air, and is thus not the typical case for projectiles. However, the linear dependence of <math>F_\mathrm{air}</math> on <math>v</math> causes a very simple differential equation of motion
:<math>\frac{\mathrm{d}}{\mathrm{d}t}\begin{pmatrix}v_x \\ v_y\end{pmatrix} = \begin{pmatrix}-\mu\,v_x \\ -g-\mu\,v_y\end{pmatrix}</math>
in which the 2 cartesian components become completely independent, and it is thus easier to solve.<ref name="AryaArya1997">{{cite book|author1=Atam P. Arya|author2=Atam Parkash Arya|title=Introduction to Classical Mechanics|url=https://books.google.com/books?id=LxWMAAAACAAJ|date=September 1997|publisher=Prentice Hall Internat.|isbn=978-0-13-906686-3|page=227}}</ref> Here, <math>v_0</math>,<math>v_x</math> and <math>v_y</math> will be used to denote the initial velocity, the velocity along the direction of <var>x</var> and the velocity along the direction of <var>y</var>, respectively. The mass of the projectile will be denoted by <var>m</var>, and <math>\mu:=k/m</math>. For the derivation only the case where <math display="inline">0^o \le \theta \le 180^o </math> is considered. Again, the projectile is fired from the origin (0,0).

{{Collapse top|title=Derivation of horizontal position}}
The relationships that represent the motion of the particle are derived by [[Newton's second law]], both in the x and y directions. 
In the x direction <math> \Sigma F = -kv_x = ma_x </math> and in the y direction <math> \Sigma F = -kv_y - mg = ma_y </math>.

This implies that:

<math> a_x = -\mu v_x = \frac{\mathrm{d}v_x}{\mathrm{d}t} </math> (1),

and

<math> a_y = -\mu v_y - g = \frac{\mathrm{d}v_y}{\mathrm{d}t} </math> (2)

Solving (1) is an elementary [[differential equation]], thus the steps leading to a unique solution for <var>v<sub>x</sub></var> and, subsequently, <var>x</var> will not be enumerated. Given the initial conditions <math> v_x = v_{x0} </math> (where <var>v<sub>x0</sub></var> is understood to be the x component of the initial velocity) and <math> x=0 </math> for <math> t=0 </math>:

<math> v_x = v_{x0} e^{-\mu t} </math> (1a)
{{Collapse bottom}}
:<math> x(t) = \frac{v_{x0}}{\mu}\left(1-e^{-\mu t}\right) </math> (1b)
{{Collapse top|title=Derivation of vertical position}}
While (1) is solved much in the same way, (2) is of distinct interest because of its non-homogeneous nature. Hence, we will be extensively solving (2). Note that in this case the initial conditions are used <math> v_y=v_{y0} </math> and <math> y=0 </math> when <math> t=0 </math>.

<math> \frac{\mathrm{d}v_y}{\mathrm{d}t} = -\mu v_y - g </math> (2)

<math> \frac{\mathrm{d}v_y}{\mathrm{d}t} + \mu v_y = - g </math> (2a)

This first order, linear, non-homogeneous differential equation may be solved a number of ways; however, in this instance, it will be quicker to approach the solution via an [[integrating factor]] <math> e^{\int \mu \,\mathrm{d}t} </math>.

<math> e^{\mu t}(\frac{\mathrm{d}v_y}{\mathrm{d}t} + \mu v_y) = e^{\mu t}(-g) </math> (2c)

<math> (e^{\mu t}v_y)^\prime = e^{\mu t}(-g) </math> (2d)

<math> \int{(e^{\mu t}v_y)^\prime \,\mathrm{d}t} = e^{\mu t}v_y = \int{ e^{\mu t}(-g) \,\mathrm{d}t} </math> (2e)

<math> e^{\mu t}v_y = \frac{1}{\mu} e^{\mu t}(-g) + C </math>(2f)

<math> v_y = \frac{-g}{\mu} + Ce^{-\mu t} </math> (2g)

And by integration we find:

<math> y = -\frac{g}{\mu}t - \frac{1}{\mu}(v_{y0} + \frac{g}{\mu})e^{-\mu t} + C </math> (3)

Solving for our initial conditions:

<math> v_y(t) = -\frac{g}{\mu} + (v_{y0} + \frac{g}{\mu})e^{-\mu t} </math> (2h)

<math> y(t) = -\frac{g}{\mu}t - \frac{1}{\mu}(v_{y0} + \frac{g}{\mu})e^{-\mu t} + \frac{1}{\mu}(v_{y0} + \frac{g}{\mu}) </math> (3a)<br />
<br />
With a bit of algebra to simplify (3a):
{{Collapse bottom}}
:<math> y(t) = -\frac{g}{\mu} t + \frac{1}{\mu} \left(v_{y0} + \frac{g}{\mu} \right) (1 - e^{-\mu t} ) </math> (3b)

{{Collapse top|title=Derivation of the time of flight}}
The total time of the journey in the presence of air resistance (more specifically, when <math>F_{air}=-kv</math>) can be calculated by the same strategy as above, namely, we solve the equation <math>y(t)=0</math>. While in the case of zero air resistance this equation can be solved elementarily, here we shall need the [[Lambert W function]]. The equation
<small><math>y(t)= -\frac{g}{\mu}t + \frac{1}{\mu}(v_{y0} + \frac{g}{\mu})(1 - e^{-\mu t}) = 0</math></small>
is of the form <small><math>c_1t+c_2+c_3e^{c_4t}=0</math></small>, and such an equation can be transformed into an equation solvable by the <small><math>W</math></small> function (see an example of such a transformation [[Lambert W function#Solving equations|here]]). Some algebra shows that the total time of flight, in closed form, is given as<ref name = 'Bernardo'>{{Cite journal | last1 = Reginald Cristian | first1 = Bernardo | last2 = Jose Perico | first2 = Esguerra | last3 = Jazmine Day | first3 = Vallejos | last4 = Jeff Jerard | first4 = Canda | title = Wind-influenced projectile motion | year = 2015 | volume = 36 | issue = 2 | journal = European Journal of Physics | page = 025016 | doi=10.1088/0143-0807/36/2/025016| bibcode = 2015EJPh...36b5016B | s2cid = 119601402 }}</ref>
{{Collapse bottom}}
:<math>t=\frac{1}{\mu} \left(1+ \frac{\mu}{g} v_{y0} + W\bigl(-(1+ \frac{\mu}{g} v_{y0}) e^{-(1+ \frac{\mu}{g} v_{y0})} \bigr)\right)</math>.

=== Trajectory of a projectile with Newton drag ===
[[File:Mplwp skydive trajectory.svg|thumb|250px|Trajectories of a [[skydiver]] in air with Newton drag: para-flight; bomber near(ly above) target. <math display="inline"> v.d=10^3.\nu=150 </math> cm<sup>2</sup>/s?]]
The most typical case of [[air resistance]], in case of [[Reynolds number]]s above about 1000, is Newton drag with a drag force proportional to the speed squared, <math>F_{\mathrm{air}} = -k v^2</math>. In air, which has a [[kinematic viscosity]] around  0.15&nbsp;cm<sup>2</sup>/s, this means that the product of object speed and diameter must be more than about 0.015 m<sup>2</sup>/s.

Unfortunately, the equations of motion can '''''not''''' be easily solved analytically for this case. Therefore, a numerical solution will be examined.

The following assumptions are made:
* Constant [[gravitational acceleration]] <math display="inline">g</math>
* Air resistance is given by the following [[drag equation|drag formula]],
::<math>\mathbf{F_D} = -\tfrac{1}{2} c \rho A\, v\,\mathbf{v}</math>

::Where:

:*''F<sub>D</sub>'' is the drag force,
:*''c'' is the [[drag coefficient]],
:*ρ is the [[air density]],
:*''A'' is the [[cross sectional area]] of the projectile. Again <math>\mu=k/m :=</math> <math>c\rho A/(2m)</math> <math>= c \rho/(2\rho_p l)</math>. Compare this with theory/practice of the [[ballistic coefficient]].

==== Special cases ====
Even though the general case of a projectile with Newton drag cannot be solved analytically, some special cases can. Here we denote the [[terminal velocity]] in free-fall as <math display="inline">v_\infty=\sqrt{g/\mu}</math> and the characteristic ''settling time'' constant <math display="inline">t_f\equiv1/(v_\infty \mu) =1/\sqrt{g\mu}</math>. (Dimension of <math>g</math> [m/s<sup>2</sup>], <math>\mu</math> [1/m])
*Near-'''horizontal''' motion: In case the motion is almost horizontal, <math>|v_x|\gg|v_y|</math>, such as a flying bullet. The vertical velocity component has very little influence on the horizontal motion. In this case:<ref name="Greiner2003">{{cite book|author=Walter Greiner|title=Classical Mechanics: Point Particles and Relativity|url=https://books.google.com/books?id=W7Rwdc6JMb8C|year=2004|publisher=Springer Science & Business Media|isbn=0-387-95586-0|page=181}}</ref>
::<math>\dot{v}_x(t) = -\mu\,v_x^2(t)</math>
::<math>v_x(t) = \frac{1}{1/v_{x,0}+\mu\,t}</math>
::<math>x(t) = \frac{1}{\mu}\ln(1+\mu\,v_{x,0}\cdot t)</math>
:The same pattern applies for motion with friction along a line in any direction, when gravity is negligible (relatively small <math>g</math>). It also applies when vertical motion is prevented, such as for a moving car with its engine off.

*Vertical motion '''upward''':<ref name="Greiner2003"/>
::<math>\dot{v}_y(t) = -g-\mu\,v_y^2(t)</math>
::<math>v_y(t) = v_\infty \tan\frac{t_{\mathrm{peak}}-t}{t_f}</math>
::<math>y(t) = y_{\mathrm{peak}} + \frac{1}{\mu} \ln \bigl(\cos\frac{t_\mathrm{peak}-t}{t_f} \bigr)</math>
:Here
::<math>v_\infty \equiv\sqrt{\frac{g}{\mu}}</math> and <math>t_f=\frac{1}{\sqrt{\mu g}},</math>
::<math>t_{\mathrm{peak}} \equiv t_f \arctan{\frac{v_{y,0}}{v_\infty}} = \frac{1}{\sqrt{\mu g}} \arctan{\left(\sqrt\frac{\mu}{g}v_{y,0}\right)},</math>
:and
::<math>y_{\mathrm{peak}} \equiv -\frac{1}{\mu}\ln{\cos{\frac{t_\mathrm{peak}}{t_f}}} = \frac{1}{2\mu} \ln{ \bigl(1+\frac{\mu}{g} v_{y,0}^2 \bigr)}</math>
:where <math>v_{y,0}</math> is the initial upward velocity at <math>t = 0</math> and the initial position is <math>y(0) = 0</math>.
:A projectile cannot rise longer than <math>t_\mathrm{rise}=\frac{\pi}{2}t_f</math> <!-- tan(t_peak/t_f) = v_y,0 sqrt(μ/g) = ? --> in the vertical direction, when it reaches the peak (0 m, y<sub>peak</sub>) at 0 m/s.

*Vertical motion '''downward''':<ref name="Greiner2003"/>
::<math>\dot{v}_y(t) = -g+\mu\,v_y^2(t)</math>
::<math>v_y(t) = -v_\infty \tanh\frac{t-t_{\mathrm{peak}}}{t_f}</math>
::<math>y(t) = y_{\mathrm{peak}} - \frac{1}{\mu} \ln(\cosh\biggl(\frac{t-t_\mathrm{peak}}{t_f}\biggr))</math> With ''hyperbolic'' functions
:After a time <math>t_f</math> at y=0, the projectile reaches almost terminal velocity <math>-v_\infty</math>.<!-- \tanh\frac{t-t_{\mathrm{peak}} }{t_f} is ca. 1 -->

==== Numerical solution ====
A projectile motion with drag can be computed generically by [[numerical methods for ordinary differential equations|numerical integration]] of the [[ordinary differential equation]], for instance by applying a [[Ordinary differential equation#Reduction to a first-order system|reduction to a first-order system]]. The equation to be solved is

:<math>\frac{\mathrm{d}}{\mathrm{d}t}\begin{pmatrix}x \\ y \\ v_x \\ v_y\end{pmatrix} = \begin{pmatrix}v_x \\ v_y \\ -\mu\,v_x\sqrt{v_x^2+v_y^2} \\ -g-\mu\,v_y\sqrt{v_x^2+v_y^2}\end{pmatrix}</math>.

This approach also allows to add the effects of speed-dependent drag coefficient, altitude-dependent air density (in product <math>c(v)\rho(y)</math>) and position-dependent gravity field <math display="inline">g(y)=g_0/(1+y/R)^2 </math> (when <math display="inline">y \ll R: g \lesssim g_0/(1+2y/R) \approx g_0(1-2y/R)</math>, is linear decrease).