Fuerza de rozamiento proporcional al cuadrado de la velocidad

Un proyectil de masa m que se dispara desde el origen con velocidad v0, haciendo un ángulo θ0 con la horizontal

Si despreciamos el empuje, las fuerzas que actúan sobre el cuerpo de masa m son:

Las ecuaciones del movimiento del cuerpo serán por tanto.

m d v x dt =mb v x v x 2 + v y 2 m d v y dt =mgmb v y v x 2 + v y 2

Las condiciones iniciales son t=0, velocidad inicial: v0x=v0·cosθ0, v0y=v0·sinθ0, posición inicial: x=0, y=0.

Este sistema de ecuaciones diferenciales se resuelve aplicando procedimientos numéricos.

b=0.0025; %parámetro
v0=60;  %velocidad incial

x0=zeros(1,4);
x0(1)=0;  % x
x0(3)=0; % y
tspan=[0,10];
fg=@(t,x)[x(2);-b*sqrt(x(2)^2+x(4)^2)*x(2); 
x(4);-9.8-b*sqrt(x(2)^2+x(4)^2)*x(4)];
opts=odeset('events',@stop_proyectil_v2);
hold on
for angulo=(10:5:45)*pi/180; %ángulo de tiro
    x0(2)=v0*cos(angulo); % vx
    x0(4)=v0*sin(angulo); % vy
    [t,x]=ode45(fg,tspan,x0,opts);
    plot(x(:,1),x(:,3)) %trayectoria
end
grid on
xlabel('x(m)')
ylabel('y(m)');
title('tiro con rozamiento')
hold off

Definimos una función para que se detenga el proceso de integración cuando el proyectil llega al suelo, ordenada y=0, en el código x(3)

function [detect,stopin,direction]=stop_proyectil_v2(t,x)
    detect=x(3); 
    stopin=1; 
    direction=-1; 
end

Actividades

Introducimos:

Se pulsa el botón titulado Nuevo

El programa interactivo traza y calcula el alcance de los proyectiles disparados con ángulos de 10, 15, 20, 25, 30, 35, 40 y 45º (en color rojo).

Compara estas trayectorias con la que seguiría el mismo proyectil disparado con un ángulo de 45º en el vacío (en color azul).

En la parte superior derecha, se muestra el alcance (aproximado) de cada uno de los proyectiles. Observamos que el máximo alcance del proyectil no se obtiene para el ángulo de disparo de 45º sino para un ángulo inferior. Y como cabía esperar, el alcance del proyectil disparado con 45º es inferior en un medio como el aire que en el vacío.


Soluciones analíticas aproximadas

Las ecuaciones del movimiento del cuerpo son:

d v x dt =bv· v x d v y dt =gbv· v y v= v x 2 + v y 2

Vamos a considerar tres casos:

En este apartado, vamos a obtener soluciones analíticas para estos tres casos

Angulos de tiro θ pequeños

Para ángulos de tiro θ pequeños, vx>>vy, v=vx. Se cumple en todos los puntos de la trayectoria. Las ecuaciones del movimiento se escriben

d v x dt =b v x 2 d v y dt =gb v x · v y

Integramos la primera ecuación diferencial, con la condición inicial de que en el instante t=0, vx=v0x=v0·cosθ

v 0x v x d v x v x 2 = 0 t bdt v x = v 0x 1+ v 0x bt

Integramos la segunda ecuación diferencial, con la condición inicial de que en el instante t=0, vy=v0y=v0·sinθ

d v y dt +b v 0x 1+ v 0x bt v y =g

Llamando z=vy y a=b·v0x. Escribimos la ecuación diferencial de la forma equivalente

dz dt + a 1+at z=g

Probamos la solución z=u·v

du dt v+ dv dt u+ a 1+at u·v=g v( du dt + a 1+at u )+ dv dt u=g

>> syms z a t g v0y;
>> dsolve('Dz+a*z/(1+a*t)=-g','z(0)=v0y')
ans =v0y/(a*t + 1) - (g*t*(a*t + 2))/(2*(a*t + 1))
>> simplify(ans)
ans =-(a*g*t^2 + 2*g*t - 2*v0y)/(2*(a*t + 1))

Conocidas las componentes vx=dx/dt y vy=dy/dt de la velocidad en función del tiempo, integramos para obtener la posición x e y de la partícula en función del tiempo t, sabiendo que parte del origen en el instante t=0

Ecuación de la trayectoria. Despejamos el tiempo en la primera ecuación y sustituimos en la segunda

t= 1 b v 0x ( e bx 1 ) y=( v 0y + g 2b v 0x ) x v 0x g 4 b 2 v 0x 2 ( e 2bx 1 )

El rango R, se obtiene cuando y=0.

( v 0y + g 2b v 0x ) R v 0x g 4 b 2 v 0x 2 ( e 2bR 1 )=0 e 2bR =2b( 1+ 2b v 0x v 0y g )R+1

Expresamos el rango en términos de la función W de Lambert, solución de la ecuación trascendente ab·x+c=d·x+f con a→e, b→2b, c→0, d2b( 1+ 2b v 0x v 0y g ) , f→1

R= 1 2b { W( 1 1+β exp( 1 1+β ) )+ 1 (1+β) } β= 2b v 0x v 0y g

Altura máxima, se obtiene dy/dx=0. La abscisa del máximo de y es

h=( v 0y + g 2b v 0x ) 1 2b v 0x ln( 2b v 0x v 0y g +1 ) v 0y 2b v 0x h= v 0y 2b v 0x { ( β+1 β )ln( β+1 )1 }

Poniendo x=R, obtenemos el tiempo de vuelo

t= 1 b v 0x ( e bx 1 ) T= 1 b v 0x { exp( 1 2 { W( 1 1+β exp( 1 1+β ) )+ 1 (1+β) } )1 }

Lanzamos un proyectil con velocidad inicial v0=9.8 m/s, con un ángulo de tiro θ=20°. Supondremos que el parámetro b=0.1 m-1

b=0.1; %parámetro
v0=9.8;  %velocidad incial
angulo=20*pi/180; %ángulo de tiro
%solución analítica aproximada 
beta=2*b*v0*cos(angulo)*v0*sin(angulo)/9.8;
z=-exp(-1/(1+beta))/(1+beta);
R=-(lambertw(-1,z)+1/(1+beta))/(2*b); %alcance
f=@(x) (v0*sin(angulo)+9.8/(2*b*v0*cos(angulo)))*x/(v0*cos(angulo))-
9.8*(exp(2*b*x)-1)/(2*b*v0*cos(angulo))^2;

%procedimiento numérico
x0=[0,v0*cos(angulo),0,v0*sin(angulo)];
tspan=[0,2*v0*sin(angulo)/9.8]; %tiempo de vuelo sin rozamiento
%x(1) es x, x(2) es v_x, x(3) es y, x(4) es v_y
fg=@(t,x)[x(2);-b*sqrt(x(2)^2+x(4)^2)*x(2); x(4);-9.8-b*sqrt(x(2)^2+x(4)^2)*x(4)];
opts=odeset('events',@stop_proyectil_v2);
[t,x]=ode45(fg,tspan,x0,opts);

hold on
plot(x(:,1),x(:,3)) %trayectoria
fplot(f,[0,R])
grid on
legend('numérico','ángulo pequeño')
xlabel('x(m)')
ylabel('y(m)');
title('Tiro con rozamiento')
hold off

Comparamos la altura máxima h con el valor obtenido integrando las ecuaciones del movimiento por procedimientos numéricos, max(x(:,3)), el tiempo de vuelo T con t(end), y el rango R con x(end,1), tal como se aprecia en el código más abajo

>> h=v0*sin(angulo)*((beta+1)*log(beta+1)/beta-1)/(2*b*v0*cos(angulo))
h =    0.4806
>> max(x(:,3))
ans =    0.4773
>> T=(exp(-(lambertw(-1,z)+1/(1+beta))/2)-1)/(b*v0*cos(angulo))
T=    0.6246
>> t(end)
ans =    0.6228
>> R 
R =    4.5436
>> x(end,1)
ans =    4.5106

La solución numérica se aproxima a la solución analítica para ángulos pequeños de tiro

Angulos de tiro grande, próximos a 90°

Para ángulos de tiro θ grandes, próximos a 90°, |vy|>>vx, v=|vy|. vy es positivo en la rama ascendente, es negativo en la rama descendente y es nulo en el punto más alto de la trayectoria. Las ecuaciones del movimiento se escriben

d v x dt =mb| v y | v x d v y dt =gb| v y | v y

Debido a que vy cambia de signo de la trayectoria ascendente a la descendente, hay que resolver este sistema de ecuaciones diferenciales en dos etapas

Trayectoria ascendente vy>0

d v x dt =b v y v x d v y dt =gb v y 2

Integramos la segunda ecuación diferencial

v 0y v y d v y g+b v y 2 = 0 t dt arctan( b g v y )arctan( b g v 0y )= bg t v y = g b tan( φ bg t ) φ=arctan( b g v 0y )

Conocido vy, integramos la primera ecuación diferencial para obtener la expresión de vx

v 0x v x d v x v x = 0 t b g b tan( φ bg t )dt ln v x v 0x =ln( cos( φ bg t ) cosφ ) v x = v 0x cosφ cos( φ bg t )

Conocidas las componentes vx y vy de la velocidad en función del tiempo, integramos para obtener la posición x, y, sabiendo que en el instante t=0, el proyectil parte del origen

Vértice de la trayectoria

Denominamos τ al tiempo que tarda la partícula en alcanzar la altura máxima, el vértice de la trayectoria

τ= φ bg

Las coordenadas del vértice de la trayectoria son

y 1 = ln( cosφ ) b x 1 = v 0x cosφ bg ln( tan( π 4 + φ 2 ) )

Las componentes de velocidad v1x y v1y son

v 1x = v 0x cosφ v 1y =0

El vértice de la trayectoria es el punto de partida para la trayectoria descendente

Trayectoria descendente vy<0

Las ecuaciones del movimiento se escriben

d v x dt =b v y v x d v y dt =g+b v y 2

Integramos la segunda ecuación diferencial

0 v y d v y g+b v y 2 = τ t dt 1 gb 0 u du 1 u 2 =(tτ)u= b g v y du 1 u 2 = 1 2 { du 1+u + du 1u }= 1 2 ln( 1+u 1u ) 1 2 gb ln( 1+ b g v y 1 b g v y )=(tτ) v y = g b exp( 2 gb (tτ) )1 exp( 2 gb (tτ) )+1

Conocida vy, integramos la primera ecuación diferencial para obtener la expresión de vx

v 1x v x d v x v x =b τ t g b exp( 2 gb (tτ) )1 exp( 2 gb (tτ) )+1 dt e 2u 1 e 2u +1 du= tanh( u ) ·du=ln( cosh( u ) ) ln( v x v 1x )=ln( cosh( gb (tτ) ) ) v x = v 0x cosφ cosh( gb (tτ) )

Conocidas las componentes vx y vy de la velocidad en función del tiempo, integramos para obtener las expresiones de x e y

Tiempo de vuelo, T que tarda en llegar al suelo, y=0

cosh( gb (tτ) )cosφ=1 e u + e u 2 = 1 cosφ e 2u 2cosφ e u +1=0 T=τ+ 1 gb ln( 1 cosφ +tanφ )

Introduciendo este tiempo T en la expresión de x obtenemos el rango R

Lanzamos un proyectil con velocidad inicial v0=9.8 m/s, con un ángulo de tiro θ=70°. Supondremos que el parámetro b=0.1 m-1

b=0.1; %parámetro
v0=9.8;  %velocidad incial
angulo=70*pi/180; %ángulo de tiro
%solución analítica aproximada 
phi=atan(sqrt(b/9.8)*v0*sin(angulo));
%coordendas del vértice de la trauyectoria
tau=phi/sqrt(9.8*b); %tiempo 
y1=-log(cos(phi))/b;
x1=v0*cos(angulo)*cos(phi)*log(tan(pi/4+phi/2))/sqrt(9.8*b);
f_X=@(t) v0*cos(angulo)*cos(phi)*log(tan(phi/2+pi/4)./
tan(phi/2+pi/4-sqrt(9.8*b)*t/2))/sqrt(9.8*b);
f_Y=@(t) log(cos(phi-sqrt(9.8*b)*t)/cos(phi))/b;
g_X=@(t) 2*v0*cos(angulo)*cos(phi)*
(atan(exp(sqrt(9.8*b)*(t-tau)))-pi/4)/sqrt(9.8*b)+x1;
g_Y=@(t) -log(cosh(sqrt(9.8*b)*(t-tau))*cos(phi))/b;
tFin=tau+log(1/cos(phi)+tan(phi))/sqrt(9.8*b); %tiempo total de vuelo

%procedimiento numérico
x0=[0,v0*cos(angulo),0,v0*sin(angulo)];
tspan=[0,2*v0*sin(angulo)/9.8]; %tiempo de vuelo sin rozamiento
fg=@(t,x)[x(2);-b*sqrt(x(2)^2+x(4)^2)*x(2); x(4);-9.8-b*sqrt(x(2)^2+x(4)^2)*x(4)];
opts=odeset('events',@stop_proyectil_v2);
[t,x]=ode45(fg,tspan,x0,opts);

hold on
plot(x(:,1),x(:,3)) %trayectoria
fplot(f_X,f_Y,[0,tau])
fplot(g_X,g_Y,[tau,tFin])
grid on
legend('numérico','ángulo grande')
xlabel('x(m)')
ylabel('y(m)');
title('Tiro con rozamiento')
hold off

Comparamos la altura máxima y1 con el valor obtenido integrando las ecuaciones del movimiento por procedimientos numéricos, max(x(:,3)), el tiempo de vuelo tFin con t(end), y el rango R con x(end,1), tal como se aprecia en el código más abajo

>> y1
y1 =    3.1173
>> max(x(:,3))
ans =    3.0336
>> tFin
tFin =    1.5965
>> t(end)
ans =    1.5732
>> R=2*v0*cos(angulo)*cos(phi)*(atan(exp(sqrt(9.8*b)*(T-tau)))-pi/4)/sqrt(9.8*b)+x1
R =    3.9180
>> x(end,1)
ans =    3.6386

Vemos que la solución numérica y la aproximación analítica para ángulos próximos a 90° se alejan bastante en comparación con el caso anterior, para ángulos de tiro pequeños

Angulos de tiro próximos a 45°

El módulo de la velocidad en la dirección horizontal es v= 2 v x , y en la dirección vertical v= 2 | v y |

Las ecuaciones del movimiento se escriben

d v x dt = 2 b v x 2 d v y dt =g 2 b| v y | v y

Debido a que vy cambia de signo de la trayectoria ascendente a la descendente, hay que resolver este sistema de ecuaciones diferenciales en dos etapas

Trayectoria ascendente vy>0

d v x dt = 2 b v x 2 d v y dt =g 2 b v y 2

Trayectoria desscendente vy<0

d v x dt = 2 b v x 2 d v y dt =g+ 2 b v y 2

La solución de la primera ecuación diferencial es la misma que hemos obtenido para ángulos de tiro pequeños, cambiando b 2 b . La solución de la segunda ecuación diferencial es la misma que hemos obtenido para ángulos de tiro grandes, haciendo el mismo cambio

Lanzamos un proyectil con velocidad inicial v0=9.8 m/s, con un ángulo de tiro θ=45°. Supondremos que el parámetro b=0.1 m-1

b=0.1*sqrt(2); % (b=0.1), parámetro para la solución analítica
v0=9.8;  %velocidad incial
%solución analítica aproximada 
angulo=45*pi/180; %ángulo de tiro
phi=atan(sqrt(b/9.8)*v0*sin(angulo));
%coordendas del vértice de la trayectoria
tau=phi/sqrt(9.8*b); %tiempo 
y1=-log(cos(phi))/b;
x1=log(1+b*v0*cos(angulo)*tau)/b;
f_X=@(t) log(1+b*v0*cos(angulo)*t)/b;
f_Y=@(t) log(cos(phi-sqrt(9.8*b)*t)/cos(phi))/b;
g_X=@(t) log(1+b*v0*cos(angulo)*t)/b;
g_Y=@(t) -log(cosh(sqrt(9.8*b)*(t-tau))*cos(phi))/b;
tFin=tau+log(1/cos(phi)+tan(phi))/sqrt(9.8*b); %tiempo total de vuelo

%procedimiento numérico
b=0.1; %parámetro para la solución numérica
x0=[0,v0*cos(angulo),0,v0*sin(angulo)];
tspan=[0,2*v0*sin(angulo)/9.8]; %tiempo de vuelo sin rozamiento
fg=@(t,x)[x(2);-b*sqrt(x(2)^2+x(4)^2)*x(2); x(4);-9.8-b*sqrt(x(2)^2+x(4)^2)*x(4)];
opts=odeset('events',@stop_proyectil_v2);
[t,x]=ode45(fg,tspan,x0,opts);

hold on
plot(x(:,1),x(:,3)) %trayectoria
fplot(f_X,f_Y,[0,tau])
fplot(g_X,g_Y,[tau,tFin])
grid on
legend('numérico','ángulo medio')
xlabel('x(m)')
ylabel('y(m)');
title('Tiro con rozamiento')
hold off

Comparamos la altura máxima y1 con el valor obtenido integrando las ecuaciones del movimiento por procedimientos numéricos, max(x(:,3)), el tiempo de vuelo tFin con t(end), y el rango R con x(end,1), tal como se aprecia en el código más abajo

>> y1
y1 =    1.8614
>> max(x(:,3))
ans =    1.7997
>> tFin
tFin =    1.2334
>> t(end)
ans =    1.2078
>> R=log(1+b*v0*cos(angulo)*tFin)/b
R =    5.6032
>> x(end,1)
ans =    5.8162

Movimiento en la dirección tangencial y normal

En este apartado, cambiamos de sistema de referencia, y escribimos las ecuaciones del movimiento en la dirección tangencial y en la dirección normal

m dv dt =mb v 2 mgsinθ m v 2 ρ =mgcosθ

donde ρ es el radio de curvatura de la trayectoria.

En el intervalo de tiempo comprendido entre t y t+dt, la dirección del vector velocidad cambia un ángulo , que es el ángulo entre las tangentes o entre las normales. El móvil se desplaza en este intervalo de tiempo un arco ds=ρ·dθ, tal como se aprecia en la figura.

dv dt = dv ds ds dt =v dv ds = v ρ dv dθ

Hemos de tener en cuenta que la curvatura de la trayectoria es negativa (figura de la derecha). La curva queda a la derecha de la tangente tomada en sentido de las x crecientes. La igualdad anterior se escribe para este caso

dv dt = v ρ dv dθ

Las ecuaciones del movimiento en la dirección tangencial y en la dirección normal se convierten en una única ecuación diferencial de primer orden.

gcosθ v dv dθ =b v 2 +gsinθ 1 v 3 dv dθ 1 v 2 tanθ= b gcosθ

Haciendo el cambo de variable u=1/v2

du dθ +2tanθ·u= 2b gcosθ

Esta ecuación es del tipo lineal (véase Puig Adam P., Curso teórico-práctico de Ecuaciones Diferenciales aplicado a la Física y Técnica. Biblioteca Matemática, 1970. págs. 29-30)

Buscamos una solución de la forma u=w(θ)·z(θ)

w dz dθ +z dw dθ +2tanθ·(w·z)= 2b gcosθ w( dz dθ +2tanθ·z )+z dw dθ = 2b gcosθ

Hacemos que

dz dθ +2tanθ·z=0

La integral se calcula fácilmente

dz z =2tanθ dz z =2 tanθ·dθ lnz=2lncosθ+lnC z= C 1 cos 2 θ

Nos queda ahora que

z dw dθ = 2b gcosθ dw dθ = 2b C 1 g cos 3 θ w= 2b C 1 g dθ cos 3 θ

Integramos por partes

dθ cos 3 θ = tanθ cosθ dθ cos 3 θ + dθ cosθ dθ cos 3 θ = 1 2 tanθ cosθ + 1 2 dθ cosθ  

Resolvemos esta última integral haciendo el cambio de variable t=tan(θ/2)

cosθ= cos 2 (θ/2) sin 2 (θ/2)= 1 1+ t 2 t 2 1+ t 2 = 1 t 2 1+ t 2 dθ= 2 1+ t 2 dt dθ cosθ = 2dt 1 t 2 = dt 1t + dt 1+t = ln 1+t 1t = ln 1+tan(θ/2) 1tan(θ/2) =ln cos(θ/2)+sin(θ/2) cos(θ/2)sin(θ/2) =ln 1+sinθ cosθ dθ cosθ =ln| tanθ+ 1 cosθ |

De este modo,

dθ cos 3 θ = 1 2 tanθ cosθ + 1 2 ln| tanθ+ 1 cosθ |

Finalmente,

w= 2b C 1 g ( 1 2 tanθ cosθ + 1 2 ln| tanθ+ 1 cosθ | )+ C 2 u=w·z= cos 2 θ{ C 2 b g ( tanθ cosθ +ln| tanθ+ 1 cosθ | ) } v 2 = 1 cos 2 θ { C 2 b g ( tanθ cosθ +ln| tanθ+ 1 cosθ | ) } 1

SE calcula la constante de integración C2 a partir de las condiciones iniciales: en el instante t=0, la velocidad de disparo es v0 y hace un ángulo θ0 con la horizontal (véase la figura más abajo)

v 2 = 1 cos 2 θ { C 2 b g ( tanθ cosθ +ln| tanθ+ 1 cosθ | ) } 1 v 0 2 = 1 cos 2 θ 0 { C 2 b g ( tan θ 0 cos θ 0 +ln| tan θ 0 + 1 cos θ 0 | ) } 1

La función que relaciona el módulo de la velocidad v y el ángulo θ, que forma la dirección de la velocidad (tangente a la trayectoria) con la horizontal es

v 2 = 1 cos 2 θ { 1 v 0x 2 b g ( f(θ)f( θ 0 ) ) } 1 f(θ)= tanθ cosθ +ln| tanθ+ 1 cosθ |

Posición del proyectil

dx=ds·cosθ=ρdθ·cosθ

Utilizando la ecuación del movimiento en la dirección normal, y teniendo en cuenta que la trayectoria tiene curvatura negativa

dx= v 2 gcosθ cosθ·dθ= v 2 g dθ x= θ 0 θ v 2 g dθ

Del mismo modo

dy=ds·sinθ=ρdθ·sinθ

dy= v 2 gcosθ sinθ·dθ= v 2 g tanθdθ y= θ 0 θ v 2 g tanθ·dθ

Tiempo de vuelo

ds=v·dt
ρdθ=
v·dt

dt= v 2 vgcosθ dθ= v gcosθ dθ t= θ 0 θ v gcosθ dθ

Calculamos el ángulo θ final que forma la dirección de la velocidad cuando y=0 (véase la figura más arriba).

θ 0 θ v 2 g tanθ·dθ =0

Conocido el ángulo final θf se calcula el alcance x y el tiempo de vuelo t, resolviendo las integrales por procedimientos numéricos

x= θ 0 θ v 2 g dθ t= θ 0 θ v gcosθ dθ

El script, calcula

angIni=45*pi/180;
v0=60.0;
%tiro parabólico
H0=v0^2*sin(angIni)^2/(2*9.8); %altura máxima
R0=v0^2*sin(2*angIni)/9.8; %alcance
T0=2*v0*sin(angIni)/9.8; %tiempo de vuelo

%rozamiento proporcional al cuadrado de la velocidad
b=0.0025;
h=@(x) tan(x)./cos(x)+log(abs(tan(x)+1./cos(x)));
inv_v2=@(x) (cos(x).^2).*(1.0/(v0*cos(angIni))^2-
b*(h(x)-h(angIni))/9.8); %inversa de v^2
f=@(x) -tan(x)./(9.8*inv_v2(x));
g=@(y) integral(f,angIni,y);
angFinal=fzero(g,[-pi*80/180, angIni-eps]);
hMax=integral(f,angIni,0);
f=@(x) -1./(9.8*inv_v2(x));
xMax=integral(f,angIni,angFinal);
f=@(x) -1./(9.8*cos(x).*sqrt(inv_v2(x)));
tVuelo=integral(f,angIni,angFinal);
fprintf('altura máxima %2.1f (%2.1f), alcance %3.1f (%3.1f), 
tiempo de vuelo %1.2f (%1.2f)\n', hMax, H0, xMax,R0, tVuelo, T0);
altura máxima 68.5 (91.8), alcance 223.3 (367.3), 
tiempo de vuelo 7.45 (8.66)

Entre paréntesis se muestran los resultados para el caso del tiro parabólico ideal (sin rozamiento)

Referencias

R. D. H. Warburton, J. Wang, J. Burgdofer. Analytic Approximations of Projectile Motion with Quadratic Air Resistence. J. Service Science & Management, 2010,3:98-105

Erlichson H. Maximum projectile range with drag and lift, with particular application to golf. Am. J. Phys. 51 (4) April 1983, pp. 357-362.

Warburton R. D. H. , Wang J., Analysis of asymptotic projectile motion with air resistance using Lambert W function. Am. J. Phys. 72 (11) November 2004, pp. 1404-1407

Brancazio P. J. Looking into Chapman's homer: The physics of judging a fly ball. Am. J. Phys. 53 (9) September 1985, pp. 849-855.