Proyectil que desliza sobre un plano inclinado

Consideremos un proyectil de masa m, que es lanzado desde el origen O con velocidad v0 haciendo un ángulo θ0 con el eje X. En un momento dado, la posición del proyectil es (x,y), su vector velocidad v es tangente a la trayectoria y forma un ángulo θ con el eje X.

Las fuerzas que actúan sobre el proyectil son:

Casos particulares

Los casos particulares más sencillos se producen cuando el ángulo de tiro θ0=π/2 y cuando θ0=-π/2.

θ0=π/2

a=g( sinα+μcosα ) v= v 0 +at y= v 0 t+ 1 2 a t 2

El proyectil se detiene, v=0, en la posición

x=0 y= 1 2 v 0 2 g( sinα+μcosα )

θ0=-π/2

a=g( sinα+μcosα ) v= v 0 +at y= v 0 t+ 1 2 a t 2

El proyectil se detiene, v=0, si μcosα>sinα, μ>tanα, en la posición

x=0 y= 1 2 v 0 2 g( sinα+μcosα )

Si no se cumple esta condición, el proyectil continuará su movimiento descendente a lo largo del eje Y del plano inclinado

Movimiento general

Escribimos las ecuaciones del movimiento en la dirección tangencial y en la dirección normal

m dv dt =mgsinαsinθμmgcosα m v 2 ρ =mgsinαcosθ

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.

gsinαcosθ v dv dθ =gsinαsinθμgcosα 1 v dv dθ =tanθ+ f cosθ

Sea f=μ/tanα

Velocidad del proyectil

Integramos la ecuación diferencial, teniendo en cuenta que en el cuerpo se lanza del origen con velocidad v0 haciendo un ángulo θ0 con la horizontal

v 0 v dv v = θ 0 θ ( tanθ+ f cosθ ) dθ ln( v v 0 )=ln( cosθ cos θ 0 )+fln( (1+sinθ)cos θ 0 cosθ(1+sin θ 0 ) ) v= v 0 ( cos θ 0 cosθ ) f+1 ( 1+sinθ 1+sin θ 0 ) f

Se ha utilizado el resultado de la integral de la inversa del coseno

Tiempo de vuelo

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

dt= v gsinαcosθ dθ= t= v 0 gsinα ( cos θ 0 ) f+1 ( 1+sin θ 0 ) f θ 0 θ ( 1+sinθ cosθ ) f dθ cos 2 θ

Haciendo los cambios de variable

u=tanθdu= dθ cos 2 θ ( 1+ u 2 +u ) f du u=sinhz ( coshz+sinhz ) f coshz·dz= 1 2 ( e (f+1)z + e (f1)z ) dz= e fz 2 ( e z f+1 + e z f1 ) e z = 1+sinθ cosθ e z = 1sinθ cosθ

llegamos al siguiente resultado

t= v 0 gsinα ( cos θ 0 ) f+1 ( 1+sin θ 0 ) f { ( 1+sinθ ) f ( cosθ ) f ( fsinθ ) ( f 2 1 )cosθ } | θ 0 θ = 1 f 2 1 v 0 gsinα { ( fsin θ 0 ) ( cos θ 0 cosθ ) f+1 ( 1+sinθ 1+sin θ 0 ) f ( fsinθ ) }

Abscisa x

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

dx= v 2 gsinαcosθ cosθ·dθ= v 2 gsinα dθ x= v 0 2 gsinα ( cos θ 0 ) 2f+2 ( 1+sin θ 0 ) 2f θ 0 θ ( 1+sinθ cosθ ) 2f dθ cos 2 θ

Haciendo los cambios de variable

u=tanθdu= dθ cos 2 θ ( 1+ u 2 +u ) 2f du u=sinhz ( coshz+sinhz ) 2f coshz·dz= 1 2 ( e (2f+1)z + e (2f1)z ) dt= e 2fz 2 ( e z 2f+1 + e z 2f1 ) e z = 1+sinθ cosθ e z = 1sinθ cosθ

Llegamos al siguiente resultado

x= v 0 2 gsinα ( cos θ 0 ) 2f+2 ( 1+sin θ 0 ) 2f { ( 1+sinθ ) 2f ( cosθ ) 2f ( 2fsinθ ) ( 4 f 2 1 )cosθ } | θ 0 θ = 1 4 f 2 1 v 0 2 gsinα { ( 2fsin θ 0 )cos θ 0 ( cos θ 0 cosθ ) 2f+2 ( 1+sinθ 1+sin θ 0 ) 2f cosθ( 2fsinθ ) }

Ordenada y

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

dy= v 2 gsinαcosθ sinθ·dθ= v 2 gsinα tanθdθ y= v 0 2 gsinα ( cos θ 0 ) 2f+2 ( 1+sin θ 0 ) 2f θ 0 θ ( 1+sinθ cosθ ) 2f tanθ cos 2 θ dθ

Haciendo los cambios de variable

u=tanθdu= dθ cos 2 θ ( 1+ u 2 +u ) 2f u·du u=sinhz ( coshz+sinhz ) 2f sinhz·coshz·dz= 1 4 ( e (2f+2)z e (2f2)z ) dt= e 2fz 8 ( e 2z f+1 e 2z f1 ) e z = 1+sinθ cosθ e z = 1sinθ cosθ

llegamos al siguiente resultado

y= v 0 2 gsinα ( cos θ 0 ) 2f+2 ( 1+sin θ 0 ) 2f { ( 1+sinθ cosθ ) 2f ( 2fsinθ1 sin 2 θ ) 4( f 2 1 ) cos 2 θ } | θ 0 θ = 1 4( f 2 1 ) v 0 2 gsinα { ( 2fsin θ 0 1 sin 2 θ 0 ) ( cos θ 0 cosθ ) 2f+2 ( 1+sinθ 1+sin θ 0 ) 2f ( 2fsinθ1 sin 2 θ ) }

Trayectorias cuando el proyectil se detiene, f>1

Cuando f>1 el proyectil se detiene v=0, la dirección de la velocidad hace un ángulo θ=-π/2.

En la expresión de la velocidad v en función del ángulo θ

v= v 0 ( cos θ 0 cosθ ) f+1 ( 1+sinθ 1+sin θ 0 ) f

obtenemos un cociente 0/0 cuando θ=-π/2. Resolvemos esta indeterminación, del siguiente modo

( 1+sinθ ) f ( cosθ ) f+1 = ( ( cos θ 2 +sin θ 2 ) 2 ) f ( ( cos θ 2 +sin θ 2 )( cos θ 2 sin θ 2 ) ) f+1 = ( cos θ 2 +sin θ 2 ) f1 ( cos θ 2 sin θ 2 ) f+1 lim θ π 2 ( 1+sinθ ) f ( cosθ ) f+1 = lim θ π 2 ( cos θ 2 +sin θ 2 ) f1 ( cos θ 2 sin θ 2 ) f+1 =0f>1

Otro resultado notable, el límite tiende a 1/2 cuando f=1

Teniendo en cuenta este resultado para f>1, el instante t1 en el que el proyectil se detiene es

t 1 = 1 f 2 1 v 0 gsinα ( fsin θ 0 )

Su posición (x1, y1) es

x 1 = 1 4 f 2 1 v 0 2 gsinα ( 2fsin θ 0 )cos θ 0 y 1 = 1 4( f 2 1 ) v 0 2 gsinα ( 2fsin θ 0 1 sin 2 θ 0 )

Representamos la velocidad del proyectil en función del tiempo hasta que se detiene. El ángulo de disparo es θ0=-30°, el ángulo del plano inclinado es α=30°, la velocidad de disparo es v0=1, y el parámetro que mide el rozamiento f=2

f=2;
alfa=pi/6;
ang=-pi/6;
v0=1;
v=@(th) v0*((cos(ang)./cos(th)).^(f+1)).*((1+sin(th))/(1+sin(ang))).^f;
t=@(th) v0*((f-sin(ang))-((cos(ang)./cos(th)).^(f+1)).*
(((1+sin(th))/(1+sin(ang))).^f).*(f-sin(th)))/((f^2-1)*9.8*sin(alfa));
fplot(t,v,[-pi/2,ang])
grid on
xlabel('t')
ylabel('v')
title('velocidad')

Calculamos el tiempo t1 que tarda en detenerse

>> t1 = v0*(f - sin(ang))/((f^2 - 1)*9.8*sin(alfa))
t1 =    0.1701

Casos particulares

Resultados que hemos obtenido al principio de esta página

Dibujamos la curva que une los puntos (x1, y1) para todos los ángulos de tiro -π/2≤θ0≤π/2 para un determinado valor de f>1. Se ha fijado la velocidad de disparo v 0 2 =gsinα

hold on
for f=[2,3,4,6,10]
    x1=@(th) (2*f-sin(th)).*cos(th)/(4*f^2-1);
    y1=@(th) (2*f*sin(th)-1-sin(th).^2)/(4*(f^2-1));
    fplot(x1,y1,[-pi/2,pi/2], 'displayName',num2str(f));
end
axis equal
hold off
grid on
legend('-DynamicLegend','location','northeast')
xlabel('x')
ylabel('y')
title('Envolventes')

El proyectil recorre más distancia cuando el ángulo de tiro apunta hacia abajo θ0=-π/2 que cuando apunta hacia arriba θ0=π/2. A medida que se incrementa el rozamiento f, las distancias hasta que se detiene el proyectil se hacen cada vez más pequeñas

Cuando f se hace grande, (x1, y1) tienden a

x 1 = 1 4 f 2 v 0 2 gsinα 2fcos θ 0 = 1 2f v 0 2 gsinα cos θ 0 y 1 = 1 4 f 2 v 0 2 gsinα 2fsin θ 0 = 1 2f v 0 2 gsinα sin θ 0 x 1 2 + y 1 2 = 1 4 f 2 v 0 4 g 2 sin 2 α

Los puntos (x1, y1) están situados en una semicircunferencia de radio r (la más interior de color verde, con la etiqueta de f=10)

r= 1 2f v 0 2 gsinα

Dibujamos las trayectorias de los proyectiles lanzados con la misma velocidad de disparo v0, con ángulos de tiro espaciados 15°=π/12. El valor del parámetro f=2. Se dibuja la curva que une los puntos (x1, y1) en los que se detiene el proyectil

f=2;
hold on
for ang=-pi/2+pi/12:pi/12:pi/2-pi/12
    x=@(th) ((2*f-sin(ang))*cos(ang)-((cos(ang)./cos(th)).^(2*f+2)).
*(((1+sin(th))/(1+sin(ang))).^(2*f)).*cos(th).*(2*f-sin(th)))/(4*f^2-1);
    y=@(th) ((2*f*sin(ang)-1-sin(ang)^2)-((cos(ang)./cos(th)).^(2*f+2)).
*(((1+sin(th))/(1+sin(ang))).^(2*f)).*(2*f*sin(th)-1-sin(th).^2))/(4*(f^2-1));
    fplot(x,y,[-pi/2,ang])
end
x1=@(th) (2*f-sin(th)).*cos(th)/(4*f^2-1);
y1=@(th) (2*f*sin(th)-1-sin(th).^2)/(4*(f^2-1));
fplot(x1,y1,[-pi/2,pi/2]);
axis equal
xlim([0,0.9])
grid on
hold off
xlabel('x')
ylabel('y')
title('trayectorias')

Utilizamos las herramientas Zoom in y Pan de la ventana gráfica, para observar el final de la trayectoria de un proyectil, cómo se curva hasta que la tangente a la trayectoria alcanza el ángulo límite θ=-π/2

Trayectorias cuando el proyectil no se detiene, f<1

Velocidad del proyectil

Representamos la velocidad del proyectil en función del tiempo, para un ángulo de tiro de θ=45°, la velocidad de disparo es v0=1, el ángulo del plano inclinado α=30°, el parámetro f=0.7

f=0.7;
alfa=pi/6;
ang=pi/4;
v=@(th) ((cos(ang)./cos(th)).^(f+1)).*((1+sin(th))/(1+sin(ang))).^f;
t=@(th) ((f-sin(ang))-((cos(ang)./cos(th)).^(f+1)).*
(((1+sin(th))/(1+sin(ang))).^f).*(f-sin(th)))/((f^2-1)*9.8*sin(alfa));
fplot(t,v,[-pi/2+0.001,ang])
grid on
xlabel('t')
ylabel('v')
title('velocidad')

Vemos que el módulo de la velocidad disminuye, hasta que alcanza la altura máxima y luego vuelve a aumentar

Trayectorias de los proyectiles

Examinamos las trayectorias de los proyectiles disparados desde el origen con ángulos espaciados 15°, para f=0.7

f=0.7;
hold on
for ang=-pi/2+pi/12:pi/12:pi/2-pi/12
    x=@(th) ((2*f-sin(ang))*cos(ang)-((cos(ang)./cos(th)).^(2*f+2)).
*(((1+sin(th))/(1+sin(ang))).^(2*f)).*cos(th).*(2*f-sin(th)))/(4*f^2-1);
    y=@(th) ((2*f*sin(ang)-1-sin(ang)^2)-((cos(ang)./cos(th)).^(2*f+2)).
*(((1+sin(th))/(1+sin(ang))).^(2*f)).*(2*f*sin(th)-1-sin(th).^2))/(4*(f^2-1));
    fplot(x,y,[-pi/2+0.0001,ang])
end
axis equal
ylim([-0.8,0.3])
xlim([0,0.9])
grid on
hold off
xlabel('x')
ylabel('y')
title('trayectorias')

Observamos que los proyectiles disparados con un ángulo próximo a 90° se curvan rápidamente y descienden siguiendo una trayectoria que es casi una línea recta paralela al eje Y. Observamos que el alcance de los proyectiles está limitado

Cuando θ tiende a -π/2, la ordenada y tiende a infinito. Sin embargo, la abscisa x tiene un comportamiento diferente

( 1+sinθ ) 2f ( cosθ ) 2f+1 = ( cos θ 2 +sin θ 2 ) 2f ( ( cos θ 2 +sin θ 2 )( cos θ 2 sin θ 2 ) ) 2f+1 = ( cos θ 2 +sin θ 2 ) 2f1 ( cos θ 2 sin θ 2 ) 2f+1 lim θ π 2 ( 1+sinθ ) 2f ( cosθ ) 2f+1 = lim θ π 2 ( cos θ 2 +sin θ 2 ) 2f1 ( cos θ 2 sin θ 2 ) 2f+1 =00.5<f1

Otro resultado notable, el límite tiende a 1/2 cuando f=0.5

Llevando este resultado a la expresión de x (aquí abajo) obtenemos x

x= 1 4 f 2 1 v 0 2 gsinα { ( 2fsin θ 0 )cos θ 0 ( cos θ 0 cosθ ) 2f+2 ( 1+sinθ 1+sin θ 0 ) 2f cosθ( 2fsinθ ) } x = 1 4 f 2 1 v 0 2 gsinα ( 2fsin θ 0 )cos θ 0

Representamos la trayectoria de un proyectil disparado con ángulo de tiro θ0=π/4, y dibujamos mediante una línea a trazos su asíntota vertical x

f=0.7;
ang=pi/4;
x=@(th) ((2*f-sin(ang))*cos(ang)-((cos(ang)./cos(th)).^(2*f+2)).
*(((1+sin(th))/(1+sin(ang))).^(2*f)).*cos(th).*(2*f-sin(th)))/(4*f^2-1);
y=@(th) ((2*f*sin(ang)-1-sin(ang)^2)-((cos(ang)./cos(th)).^(2*f+2)).
*(((1+sin(th))/(1+sin(ang))).^(2*f)).*(2*f*sin(th)-1-sin(th).^2))/(4*(f^2-1));
fplot(x,y,[-pi/2+0.0001,ang])
x_inf=(2*f-sin(ang))*cos(ang)/(4*f^2-1);
line([x_inf,x_inf],[-0.8,0.4],'linestyle','--','color','k')
axis equal
ylim([-0.8,0.4])
xlim([0,1.3])
grid on
xlabel('x')
ylabel('y')
title('trayectorias')

La altura máxima de un proyectil disparado con ángulo de tiro, θ0>0, se obtiene para θ=0

>> y(0)
ans =    0.1787

Cambiamos el parámetro f=0.3, y volvemos a dibujar las trayectorias

f=0.3;
hold on
for ang=-pi/2+pi/12:pi/12:pi/2-pi/12
    x=@(th) ((2*f-sin(ang))*cos(ang)-((cos(ang)./cos(th)).^(2*f+2)).
*(((1+sin(th))/(1+sin(ang))).^(2*f)).*cos(th).*(2*f-sin(th)))/(4*f^2-1);
    y=@(th) ((2*f*sin(ang)-1-sin(ang)^2)-((cos(ang)./cos(th)).^(2*f+2)).
*(((1+sin(th))/(1+sin(ang))).^(2*f)).*(2*f*sin(th)-1-sin(th).^2))/(4*(f^2-1));
    fplot(x,y,[-pi/2+0.0001,ang])
end
axis equal
ylim([-0.8,0.4])
xlim([0,1.3])
grid on
hold off
xlabel('x')
ylabel('y')
title('trayectorias')

Casos especiales f=1 y f=0.5

Examinamos primeros el caso f=1, la velocidad del proyectil es

v= v 0 ( cos θ 0 cosθ ) 2 ( 1+sinθ 1+sin θ 0 )

Hemos demostrado que el límite cuando θ tiende a -π/2

lim θ π 2 ( 1+sinθ ) ( cosθ ) 2 = lim θ π 2 1 ( cos θ 2 sin θ 2 ) 2 = 1 2

La velocidad tiende al valor constante

v = v 0 2 cos 2 θ 0 1+sin θ 0

Para f=1, tenemos que volver a calcular la expresión del tiempo t y la expresión de la ordenada y, debido a que en el denominador contienen el término (f2-1) que se hace cero

Tiempo t

dt= v gsinαcosθ dθ= t= v 0 gsinα ( cos θ 0 ) 2 ( 1+sin θ 0 ) θ 0 θ ( 1+sinθ cosθ ) dθ cos 2 θ

Haciendo los cambios de variable

u=tanθdu= dθ cos 2 θ ( 1+ u 2 +u ) du u=sinhz ( coshz+sinhz ) coshz·dz= 1 2 ( e 2z +1 ) dz= 1 2 ( e 2z 2 +z ) e z = 1+sinθ cosθ

llegamos al siguiente resultado

t= v 0 gsinα ( cos θ 0 ) 2 ( 1+sin θ 0 ) 1 2 { 1 2 ( 1+sinθ ) 2 ( cosθ ) 2 +ln| 1+sinθ cosθ | } | θ 0 θ = 1 2 v 0 gsinα ( cos θ 0 ) 2 ( 1+sin θ 0 ) { 1 2 ( 1+sin θ 0 ) 2 ( cos θ 0 ) 2 +ln| 1+sin θ 0 cos θ 0 | 1 2 ( 1+sinθ ) 2 ( cosθ ) 2 ln| 1+sinθ cosθ | }

Representamos la velocidad del proyectil en función del tiempo, para un ángulo de tiro θ0=30°, la velocidad de disparo es v0=1, el ángulo del plano inclinado α=30°, el parámetro f=1

f=1;
alfa=pi/6;
ang=pi/6;
v=@(th) ((cos(ang)./cos(th)).^(f+1)).*((1+sin(th))/(1+sin(ang))).^f;
t=@(th) cos(ang)^2*(((1+sin(ang))/cos(ang))^2/2+log((1+sin(ang))/cos(ang))
-((1+sin(th))./cos(th)).^2/2-log((1+sin(th))./cos(th)))
/(2*9.8*sin(alfa)*(1+sin(ang)));
fplot(t,v,[-pi/2+0.001,ang])
grid on
ylim([0,1])
xlabel('t')
ylabel('v')
title('velocidad')

Observamos que la velocidad tiende hacia un valor límite constante

0.5*cos(ang)^2/(1+sin(ang))
ans =    0.2500

Ordenada y

dy= v 2 gsinαcosθ sinθ·dθ= v 2 gsinα tanθdθ y= v 0 2 gsinα ( cos θ 0 ) 4 ( 1+sin θ 0 ) 2 θ 0 θ ( 1+sinθ cosθ ) 2 tanθ cos 2 θ dθ

Haciendo los cambios de variable

u=tanθdu= dθ cos 2 θ ( 1+ u 2 +u ) 2 u·du u=sinhz ( coshz+sinhz ) 2 sinhz·coshz·dz= 1 4 ( e 4z 1 ) dt= 1 4 ( e 4z 4 z ) e z = 1+sinθ cosθ

Llegamos al siguiente resultado

y= 1 4 v 0 2 gsinα ( cos θ 0 ) 4 ( 1+sin θ 0 ) 2 { 1 4 ( 1+sinθ cosθ ) 4 ln| 1+sinθ cosθ | } | θ 0 θ = 1 4 v 0 2 gsinα ( cos θ 0 ) 4 ( 1+sin θ 0 ) 2 { 1 4 ( 1+sin θ 0 cos θ 0 ) 4 ln| 1+sin θ 0 cos θ 0 | 1 4 ( 1+sinθ cosθ ) 4 +ln| 1+sinθ cosθ | }

La expresión de la abscisa x sigue siendo la misma que hemos deducido para los otros casos

x= 1 4 f 2 1 v 0 2 gsinα { ( 2fsin θ 0 )cos θ 0 ( cos θ 0 cosθ ) 2f+2 ( 1+sinθ 1+sin θ 0 ) 2f cosθ( 2fsinθ ) }

Dibujamos las trayectorias de los proyectiles lanzados con la misma velocidad de disparo v0, con ángulos de tiro espaciados 15°=π/12. El valor del parámetro f=1.

f=1;
hold on
for ang=-pi/2+pi/12:pi/12:pi/2-pi/12
    x=@(th) ((2*f-sin(ang))*cos(ang)-((cos(ang)./cos(th)).^(2*f+2)).
*(((1+sin(th))/(1+sin(ang))).^(2*f)).*cos(th).*(2*f-sin(th)))/(4*f^2-1);
    y=@(th) cos(ang)^4*(((1+sin(ang))/cos(ang))^4/4-log((1+sin(ang))/cos(ang))-
((1+sin(th))./cos(th)).^4/4+log((1+sin(th))./cos(th)))/(4*(1+sin(ang))^2);
    fplot(x,y,[-pi/2,ang])
end
axis equal
ylim([-0.8,0.3])
xlim([0,0.7])
grid on
hold off
xlabel('x')
ylabel('y')
title('trayectorias')

La abscisa x está limitada cuando θ tiende a -π/2, situación que ya observamos para 0.5<f≤1

f=0.5

Examinamos ahora el caso f=0.5. El denominador de la abscisa x contiene el término (4f2-1), que se hace cero cuando f=0.5. Es preciso volver a calcular la abscisa x

dx= v 2 gsinαcosθ cosθ·dθ= v 2 gsinα dθ x= v 0 2 gsinα ( cos θ 0 ) 3 ( 1+sin θ 0 ) θ 0 θ ( 1+sinθ cosθ ) dθ cos 2 θ

Haciendo los cambios de variable

u=tanθdu= dθ cos 2 θ ( 1+ u 2 +u ) du u=sinhz ( coshz+sinhz ) coshz·dz= 1 2 ( e 2z +1 ) dt= 1 2 ( e 2z 2 +z ) e z = 1+sinθ cosθ

Llegamos al resultado

x= v 0 2 gsinα cos 3 θ 0 ( 1+sin θ 0 ) 1 2 { 1 2 ( 1+sinθ cosθ ) 2 +ln| 1+sinθ cosθ | } | θ 0 θ = 1 2 v 0 2 gsinα cos 3 θ 0 ( 1+sin θ 0 ) { 1 2 ( 1+sin θ 0 cos θ 0 ) 2 +ln| 1+sin θ 0 cos θ 0 | 1 2 ( 1+sinθ cosθ ) 2 ln| 1+sinθ cosθ | }

Las expresiones de la ordenada y y del tiempo t no han cambiado

Dibujamos las trayectorias de los proyectiles lanzados con la misma velocidad de disparo v0, con ángulos de tiro espaciados 15°=π/12. El valor del parámetro f=0.5.

f=0.5;
hold on
for ang=-pi/2+pi/12:pi/12:pi/2-pi/12
     x=@(th) cos(ang)^3*(((1+sin(ang))/cos(ang))^2/2+log((1+sin(ang))/cos(ang))-
((1+sin(th))./cos(th)).^2/2-log((1+sin(th))./cos(th)))/(2*(1+sin(ang)));
    y=@(th) ((2*f*sin(ang)-1-sin(ang)^2)-((cos(ang)./cos(th)).^(2*f+2)).
*(((1+sin(th))/(1+sin(ang))).^(2*f)).*(2*f*sin(th)-1-sin(th).^2))/(4*(f^2-1));
    fplot(x,y,[-pi/2,ang])
end
axis equal
ylim([-0.8,0.3])
xlim([0,1.1])
grid on
hold off
xlabel('x')
ylabel('y')
title('trayectorias')

Referencias

Xiaosun Wang. Trajectory of a projectile on a frictional inclined plane. Am. J. Phys. 82 (8) August 2014, pp. 764-768