Trayectoria de un proyectil disparado desde una altura h sobre la superficie de la Tierra.

Se dispara un proyectil de masa m desde una distancia r0=R+h del centro de la Tierra, con velocidad v0 haciendo un ángulo φ con el radio vector. El momento angular y la energía del proyectil son, respectivamente

L=m r 0 v 0 sinφ E= 1 2 m v 0 2 GMm r 0

La ecuación de la trayectoria en coordenadas polares es

r= d 1+εcos( θ θ 0 ) ε= 1+ 2 L 2 E G 2 M 2 m 3 d= L 2 GM m 2

La semidistancia focal, c=ε·a

El semieje menor b de la elipse b= a 2 c 2

Velocidad del proyectil en el punto de impacto

Como la energía es constante en todos los puntos de la trayectoria, la velocidad v con la que impacta el proyectil en la superficie de la Tierra es independiente de la masa m del proyectil y del ángulo de disparo. Se obtiene poniendo r=R (radio de la Tierra) en la ecuación de la energía, y despejando la incógnita v

E= 1 2 m v 0 2 GMm r 0 = 1 2 m v 2 GMm R

Tiempo de vuelo

Para calcular el tiempo de vuelo, vamos a utilizar el mismo procedimiento que empleamos para deducir la fórmula del periodo de un planeta, a partir de la ley de las áreas. El momento angular en coordenadas polares se escribe

L=m r 2 dθ dt

Integrando

θ π 1 2 r 2 dθ = 0 t L 2m dt

El primer miembro, es el área barrida por el radio vector cuando se mueve desde la posición angular θ, a la posición θ=π. Despejando t se obtiene.

t= 2m·Area L

Vamos ahora a estudiar los distintos casos que se pueden presentar

El ángulo de disparo es φ=0º.

El momento angular L=0, por lo que la trayectoria es una línea recta que pasa por el centro de fuerzas. El proyectil asciende y luego cae hacia la Tierra a lo largo de la dirección radial.

La máxima altura que alcanza, se calcula poniendo v=0 en la ecuación de la energía y se despejando la incógnita r.

E= 1 2 m v 0 2 GMm r 0 = GMm r

En la página titulada Movimiento de dos cuerpos bajo la fuerza de atracción mutua calculamos el tiempo que tarda el proyectil en impactar sobre la superficie de la Tierra

Ejemplo

Lanzamos un proyectil desde la altura h=6000 km con velocidad inicial v0= 4500 m/s en la dirección radial r0=6.0·106+6.37·106 m

El ángulo de disparo es φ=180º.

El momento angular L=0, por lo que la trayectoria es una línea recta que pasa por el centro de fuerzas. El proyectil desciende a lo largo de la dirección radial hasta que llega a la superficie de la Tierra con la misma velocidad que hemos calculado en el apartado anterior.

Ejemplo

Lanzamos un proyectil desde la posición r0=6.0·106+6.37·106 m con velocidad inicial v0= 4500 m/s en la dirección radial y sentido hacia el centro de la Tierra

El ángulo de disparo es φ=90º.

Alcance máximo

El alcance máximo se produce cuando el perigeo es R, y el apogeo es r0=h+R.

Como el momento angular y la energía son constantes en todos los puntos de la trayectoria y en particular, en el perigeo y en el apogeo, tenemos que

m r 0 v 0 ·sin90=mRv·sin90 1 2 m v 0 2 GMm r 0 = 1 2 m v 2 GMm r

Los datos son r0 y R y las incógnitas v y v0.  La velocidad de disparo es

v 0 = 2GMR r 0 (R+ r 0 )

Ejemplo: Sea h=6000 km o bien, la distancia a lo largo de la dirección radial es r0=12.37·106 m

Posición del punto de impacto

Como vemos en la figura, el proyectil sale de la posición θ=π, e impacta en la posición θ=π-α cuando r=R.

Poniendo r=R en la ecuación de la trayectoria, despejamos el ángulo θ.

cosθ= 1 ε ( d R 1 )

Ejemplo:

Continuando con los mismos datos de los casos anteriores:

Obtenemos los valores del momento angular y de la energía del proyectil

L=5.57·1010 m kgm2/s
E
=-22.12·106 m J

Conocida la energía y el momento angular, se determina la ecuación de la trayectoria, el valor del parámetro d y la excentricidad ε

ε=0.372
d
=7.77·106 m

Con estos datos, poniendo r=6.37·106 m en la ecuación de la trayectoria obtenemos el ángulo θ=0.934 rad.

La distancia angular entre el punto de impacto y la posición de disparo es

α=π-0.934=2.20 rad

Denominado alcance a la longitud del arco s de circunferencia de la Tierra que corresponde a esta distancia angular, s=R·α=14.03·106 m

Tiempo de vuelo empleando la constancia del momento angular

El área sombreada es el área barrida por el radio vector entre las posiciones angulares θ y π. En otras palabras, es la porción de elipse comprendida entre x y a menos el área del triángulo de base R·cosθ y altura R·sinθ, siendo x=-c-R·cosθ

Sabiendo que la ecuación de la elipse es

x 2 a 2 + y 2 b 2 =1

donde a es el semieje mayor de la elipse, b el semieje menor, y c la semidistancia focal.

El área de la porción de elipse comprendida entre x y a es

x a b 1 x 2 a 2 dx = z 1 z 2 ab cos 2 z·dz = ab 2 [ z+ 1 2 sin2z ] z 1 z 2 = ab 2 ( π 2 z 1 1 2 sin2 z 1 )

Para integrar, se ha hecho el cambio de variable x=a·sin z. Los nuevos límites de integración son:

El área sombreada vale, por tanto

Area= ab 2 ( π 2 z 1 1 2 sin2 z 1 ) R 2 sin2θ 4

Para calcular el área necesitamos los siguientes datos

A continuación, obtenemos z1 que es a función del ángulo θ=0.934 rad de la posición de impacto. Después de hacer algunas operaciones obtenemos el valor del área barrida por el radio vector A=1.022·1014.

Finalmente, el tiempo de vuelo t es

t= 2mA L =3671.5s

R=6.37e6; %radio Tierra
M=5.98e24; %masa de la Tierra
G=6.67e-11; %constante G

r0=R+6e6; %6000 km de altura
v0=4500; %velocidad de disparo
ang=pi/2;  %ángulo de tiro

L=r0*v0*sin(ang); %momento angular
d=L^2/(G*M);
E=v0^2/2-G*M/r0; %energía
ex=sqrt(1+2*L^2*E/(G*M)^2); %excentricidad
a=d/(1-ex^2);  %semieje mayor de la elipse
b=a*sqrt(1-ex^2); %semieje menor de la elipse
angulo=acos((d-R)/(R*ex));
z1=asin((-R*cos(angulo)-ex*a)/a);
area=a*b*(pi/2-z1-sin(2*z1)/2)/2-R^2*sin(2*angulo)/4;
fprintf('El tiempo de vuelo (s) es %4.1f\n',2*area/L)
El tiempo de vuelo (s) es 3671.5

Cálculo del tiempo de vuelo empleando la ecuación de Kepler

Dada la posición angular θ, determinamos el ángulo E mediante las relaciones

sinE= 1 ε 2 sinθ 1+εcosθ cosE= ε+cosθ 1+εcosθ

Dependiendo del signo del seno y del coseno, se determina el cuadrante y se calcula el ángulo E. Sustituyendo en la ecuación de Kepler, obtenemos el tiempo t

t= E 0 εsin E 0 2π P P 2 = 4 π 2 GM a 3

Creamos un script para calcular el tiempo de vuelo, T=t1-t0

R=6.37e6; %radio Tierra
M=5.98e24; %masa de la Tierra
G=6.67e-11; %constante G


r0=R+6e6; %6000 km de altura
v0=4500; %velocidad de disparo
ang=pi/2;  %ángulo de tiro

L=r0*v0*sin(ang); %momento angular
d=L^2/(G*M);
E=v0^2/2-G*M/r0; %energía
ex=sqrt(1+2*L^2*E/(G*M)^2); %excentricidad
a=d/(1-ex^2);  %semieje mayor de la elipse

angulo=acos((d-R)/(R*ex));
th=angulo; %posición final
c_E=(ex+cos(th))/(1+ex*cos(th));
s_E=sqrt(1-ex^2)*sin(th)/(1+ex*cos(th));
E=atan2(s_E,c_E);
if E<0
    E=2*pi+E;
end
t0=(E-ex*s_E)*a^(3/2)/sqrt(G*M);

th=pi; 
c_E=(ex+cos(th))/(1+ex*cos(th));
s_E=sqrt(1-ex^2)*sin(th)/(1+ex*cos(th));
E=atan2(s_E,c_E);
if E<0
    E=2*pi+E;
end
t1=(E-ex*s_E)*a^(3/2)/sqrt(G*M);
fprintf('El tiempo de vuelo (s) es %4.1f\n',(t1-t0))
El tiempo de vuelo (s) es 3671.5

El ángulo de disparo es φ<90º

Trayectoria

Supongamos que el ángulo de disparo es φ distinto de 0º, 90º, ó 180º.

Como vemos en la figura, la trayectoria que sigue el proyectil es una elipse, pero que está girada un ángulo β. Este ángulo se calcula poniendo r=r0 en la ecuación de la trayectoria y despejando el ángulo θ

cosβ= 1 ε ( d r 0 1 )

Continuando con los datos de los casos anteriores

La energía del proyectil no cambia, pero cambia el momento angular

L=2.78·1010 m kgm2/s
E
=-22.12·106 m J

Conocida la energía y el momento angular se determina la ecuación de la trayectoria, el valor del parámetro d y la excentricidad ε

ε=0.886
d
=1.94·106 m

Con estos datos calculamos el ángulo girado por el eje mayor de la elipse β=2.83 rad.

Posición del punto de impacto

Como vemos en la figura, calculamos el ángulo de impacto poniendo en la ecuación de la elipse r=R, lo que nos da el ángulo θ señalado en la figura, del mismo modo que en el caso anterior

cosθ= 1 ε ( d R 1 )

Relacionamos los ángulos θ, α y β. para calcular la distancia angular α entre el punto de impacto y la posición de disparo.

α=2π-θ-β

Ejemplo:

con los datos anteriores θ=2.47, y β=2.83 rad, la distancia angular α=0.981 rad (56.2º)

Tiempo de vuelo empleando la constancia del momento angular

El tiempo de vuelo es proporcional a la suma de las áreas sombreadas de elipse

Las áreas se calculan como en el caso anterior. En primer lugar, necesitamos los valores de los parámetros de la elipse:

Calculamos el área de la porción de elipse por encima del eje mayor, que es el área barrida por el radio vector desde la posición angular θ=2.47  hasta θ=π. Necesitamos conocer previamente, z1, que es a su vez función del ángulo θ de la posición de impacto.

-R·cosθ-c=a·sin z1

Area= ab 2 ( π 2 z 1 1 2 sin2 z 1 ) R 2 sin2θ 4

El resultado es A1=5.1786·1013

Calculamos el área por debajo del eje mayor barrida por el radio vector desde la posición angular β=2.83 rad  hasta β =π.

Necesitamos conocer previamente, z1, que es a su vez función del ángulo β=2.83 rad que reemplaza al ángulo θ en la fórmula del área y r0 reemplaza a R

-r0·cosβ-c=a·sin z1

Area= ab 2 ( π 2 z 1 1 2 sin2 z 1 ) r 0 2 sin2β 4

El resultado es A2=3.6620·1013

El tiempo de vuelo es

t= 2m( A 1 + A 2 ) L =6352.7s

R=6.37e6; %radio Tierra
M=5.98e24; %masa de la Tierra
G=6.67e-11; %constante G

r0=R+6e6; %6000 km de altura
v0=4500; %velocidad de disparo
ang=pi/6;  %ángulo de tiro

L=r0*v0*sin(ang); %momento angular
d=L^2/(G*M);
E=v0^2/2-G*M/r0; %energía
ex=sqrt(1+2*L^2*E/(G*M)^2); %excentricidad
a=d/(1-ex^2);  %semieje mayor de la elipse
b=a*sqrt(1-ex^2); %semieje menor de la elipse
angulo=acos((d-R)/(R*ex));
z1=asin((-R*cos(angulo)-ex*a)/a);
area1=a*b*(pi/2-z1-sin(2*z1)/2)/2-R^2*sin(2*angulo)/4;
beta=acos((d-r0)/(r0*ex));
z1=asin((-r0*cos(beta)-ex*a)/a);
area2=a*b*(pi/2-z1-sin(2*z1)/2)/2-r0^2*sin(2*beta)/4;
fprintf('El tiempo de vuelo (s) es %4.1f\n',2*(area1+area2)/L)
El tiempo de vuelo (s) es 6351.5

Tiempo de vuelo empleando la ecuación de Kepler

tal como apreciamos en la figura más arriba

R=6.37e6; %radio Tierra
M=5.98e24; %masa de la Tierra
G=6.67e-11; %constante G

r0=R+6e6; %6000 km de altura
v0=4500; %velocidad de disparo
ang=pi/6;  %ángulo de tiro

L=r0*v0*sin(ang); %momento angular
d=L^2/(G*M);
E=v0^2/2-G*M/r0; %energía
ex=sqrt(1+2*L^2*E/(G*M)^2); %excentricidad
a=d/(1-ex^2);  %semieje mayor de la elipse

angulo=acos((d-R)/(R*ex));
th=angulo; %posición inicial
c_E=(ex+cos(th))/(1+ex*cos(th));
s_E=sqrt(1-ex^2)*sin(th)/(1+ex*cos(th));
E=atan2(s_E,c_E);
if E<0
    E=2*pi+E;
end
t0=(E-ex*s_E)*a^(3/2)/sqrt(G*M);

beta=acos((d-r0)/(r0*ex));
th=2*pi-beta; %posición final
c_E=(ex+cos(th))/(1+ex*cos(th));
s_E=sqrt(1-ex^2)*sin(th)/(1+ex*cos(th));
E=atan2(s_E,c_E);
if E<0
    E=2*pi+E;
end
t1=(E-ex*s_E)*a^(3/2)/sqrt(G*M);
fprintf('El tiempo de vuelo (s) es %4.1f\n',(t1-t0))
El tiempo de vuelo (s) es 6351.5

El ángulo de disparo es φ>90º.

Trayectoria

Los proyectiles disparados con ángulos φ y 180-φ tienen la misma energía y el mismo momento angular, la trayectoria es una elipse con los mismos valores del parámetro d, y de la excentricidad ε, pero su orientación es distinta.

Si el ángulo de disparo es 150º, la energía y el momento angular son los mismos que cuando se dispara el proyectil con 30º

ε=0.886
d
=1.94·106 m

Como vemos en la figura, la trayectoria que sigue el proyectil es una elipse, pero que está girada un ángulo β. Este ángulo se calcula poniendo r=r0 en la ecuación de la trayectoria

cosβ= 1 ε ( d r 0 1 )

Con estos datos calculamos el ángulo β=2.83 rad (color rojo) que gira el eje mayor de la elipse que es la solución que tomamos en el caso anterior, pero también es solución el ángulo β=-2.83=3.45 rad (color azul)

Posición del punto de impacto

En el apartado anterior, calculamos el ángulo de impacto poniendo en la ecuación de la elipse r=R, lo que nos daba el ángulo θ=2.47 rad

cosθ= 1 ε ( d R 1 )

Relacionamos los ángulos θ, α y β. para calcular el ángulo de impacto α.

α+θ+β-π =π o bien,

α=2π-β-θ=0.36 rad (20.4º)

que es la misma relación que obtuvimos en el caso anterior.

Tiempo de vuelo empleando la constancia del momento angular

El área barrida por el radio vector desde la posición inicial de salida a la de impacto es la diferencia de dos áreas

Estas dos áreas coinciden con las áreas A1 y A2 calculadas en el caso anterior

El tiempo de vuelo es

t= 2m( A 1 A 2 ) L =1089.7s

Tiempo de vuelo empleando la ecuación de Kepler

R=6.37e6; %radio Tierra
M=5.98e24; %masa de la Tierra
G=6.67e-11; %constante G

r0=R+6e6; %6000 km de altura
v0=4500; %velocidad de disparo
ang=5*pi/6;  %ángulo de tiro

L=r0*v0*sin(ang); %momento angular
d=L^2/(G*M);
E=v0^2/2-G*M/r0; %energía
ex=sqrt(1+2*L^2*E/(G*M)^2); %excentricidad
a=d/(1-ex^2);  %semieje mayor de la elipse

angulo=acos((d-R)/(R*ex));
th=angulo; %posición inicial
c_E=(ex+cos(th))/(1+ex*cos(th));
s_E=sqrt(1-ex^2)*sin(th)/(1+ex*cos(th));
E=atan2(s_E,c_E);
if E<0
    E=2*pi+E;
end
t0=(E-ex*s_E)*a^(3/2)/sqrt(G*M);

beta=2*pi-acos((d-r0)/(r0*ex));
th=2*pi-beta; %posición final
c_E=(ex+cos(th))/(1+ex*cos(th));
s_E=sqrt(1-ex^2)*sin(th)/(1+ex*cos(th));
E=atan2(s_E,c_E);
if E<0
    E=2*pi+E;
end
t1=(E-ex*s_E)*a^(3/2)/sqrt(G*M);
fprintf('El tiempo de vuelo (s) es %4.1f\n',(t1-t0))
El tiempo de vuelo (s) es 1088.5

Cálculo con MATLAB

Creamos un script de MATLAB que se aplica a los tres casos analizados en esta página, calcula el alcance, el tiempo de vuelo y representa la trayectoria, un arco de elipse

R=6.37e6; %radio de la Tierra
M=5.98e24; %masa de la Tierra
G=6.67e-11; %constante G

r0=12.37e6; %(6000 km de altura) posición de disparo
v0=4500; %velocidad de disparo
ang=30;  %ángulo de tiro en grados
L=r0*v0*sin(ang*pi/180); %momento angular
E=v0^2/2-G*M/r0; %energía
ex=sqrt(1+2*L^2*E/(G*M)^2); %excentricidad
d=L^2/(G*M);
%ángulo de giro
beta=acos((d/r0-1.0)/ex);
if ang==90
   beta=pi;
   if v0>sqrt(G*M/r0)
        beta=0.0;
   end  
end    
    
angulo=acos((d-R)/(R*ex));
alcance=(2*pi-angulo-beta); %alcance en radianes
if ang>90
    alcance=(-angulo+beta); %alcance
end
fprintf('Alcance (km) es %2.1f\n',alcance*R/1000)

a=d/(1-ex^2);  %semieje mayor de la elipse
b=a*sqrt(1-ex^2); %semieje menor
z1=asin((-ex*a-R*cos(angulo))/a);
area1=a*b*(pi/2-z1-sin(2*z1)/2)/2-R^2*sin(2*angulo)/4;
z1=asin((-ex*a-r0*cos(beta))/a);
area2=a*b*(pi/2-z1-sin(2*z1)/2)/2-r0^2*sin(2*beta)/4;
if ang<90
    P=2*(area1+area2)/L;
elseif ang>90
    P=2*(area1-area2)/L;
else
    P=2*area1/L; %tiempo de vuelo
end
fprintf('El tiempo de vuelo (min) es %4.1f\n',P/60)

hold on
ax=(1:360)*pi/180;
fill(cos(ax),sin(ax),'c') %la Tierra
%trayectoria
r=@(x) (d/R)./(1+ex*cos(2*pi-x-beta));
if ang>90
    r=@(x) (d/R)./(1+ex*cos(-x+beta));
end
fplot(@(x) r(x).*cos(x), @(x) r(x).*sin(x),[0,2*pi],'lineStyle','--')
fplot(@(x) r(x).*cos(x), @(x) r(x).*sin(x),[0,alcance])
line([0,r0/R],[0,0],'lineStyle','--')
plot(r0/R,0,'ro','markersize',3,'markeredgecolor','r','markerfacecolor','r')
line([0,cos(alcance)],[0,sin(alcance)],'lineStyle','--')
plot(cos(alcance),sin(alcance),'ro','markersize',3,'markeredgecolor',
'r','markerfacecolor','r')
%quiver(r0/R,0,v0*cos(ang*pi/180)/1000,v0*sin(ang*pi/180)/1000);
hold off
axis equal
xlabel('x')
ylabel('y')
title('Disparo de un proyectil')

Alcance (km) es 6246.1
El tiempo de vuelo (min) es 105.9

Actividades

Se introduce

Se pulsa el botón titulado Nuevo

Se excluyen los ángulos 0 y 180º ya que su análisis es más simple y dan lugar, a errores por desbordamiento en la rutina principal de cálculo.

Se observa el movimiento del proyectil y se proporcionan los datos de la distancia angular entre el punto de impacto sobre la superficie de la Tierra y el lugar del lanzamiento, así como el tiempo de vuelo empleado por el proyectil.

Si la velocidad es grande, puede ocurrir que el proyectil se ponga en órbita alrededor de la Tierra.

Como ejercicios se sugiere,



Alcance máximo

En este apartado vamos a describir el problema del disparo de un proyectil con otra ecuación de la trayectoria cuya deducción resumimos a continuación:

Partimos de la ecuación del movimiento del proyectil en coordenadas polares

m( d 2 r d t 2 r ( dθ dt ) 2 )=G Mm r 2 m( r d 2 θ d t 2 +2 dθ dt dr dt )=0

La segunda ecuación nos indica que el momento angular L (expresado en coordenadas polares es constante)

d dt ( m r 2 dθ dt )=0L=m r 2 dθ dt =cte

La ecuación del movimiento se convierte en

d 2 u d θ 2 +u= GM m 2 L 2 u= 1 r

La solución de la ecuación diferencial es

u=Asinθ+Bcosθ+ GM m 2 L 2

Donde A, B y el momento angular L se determinan a partir de las condiciones iniciales.

El momento angular es L=mr0v0sinφ

La posición inicial en coodenadas polares es θ=0, r=r0

Las componentes de la velocidad en coordendas polares son:

Teniendo en cuenta la relación

dr dt = dr dθ dθ dt θ=0{ ( dr dt ) 0 = v 0 cosφ r 0 ( dθ dt ) 0 = v 0 sinφ ( dr dθ ) 0 = r 0 cosφ sinφ

Derivando con respecto a θ la ecuación de la trayectoria

1 r 2 dr dθ =AcosθBsinθ

Para la posición inicial θ=0.

1 r 0 2 ( dr dθ ) 0 =A A= cosφ r 0 sinφ

Expresamos la ecuación de la trayectoria en términos del ángulo de disparo φ y de un parámetro p relacionado con la velocidad de disparo v0 en la posición de disparo a una distancia r0 del centro de la Tierra

1 r = cosφ r 0 sinφ sinθ+( 1 r 0 GM r 0 2 v 0 2 sin 2 φ )cosθ+ GM r 0 2 v 0 2 sin 2 φ r 0 r = sin( φθ ) sinφ + 1cosθ p sin 2 φ p= r 0 v 0 2 GM

Alcance

El alcance angular del proyectil, α, es el ángulo θ que se obtiene para r=R

r 0 R = sin( φα ) sinφ + 1cosα p sin 2 φ

Despejamos el alcance α

psinφcosφsinα=(p sin 2 φ1)cosα+1 r 0 R p sin 2 φ

Elevamos al cuadrado obteniendo una ecuación de segundo grado en cosα

cosα= b± b 2 4ac 2a { a= ( p sin 2 φ1 ) 2 + p 2 sin 2 φ cos 2 φ c= ( 1 r 0 R p sin 2 φ ) 2 p 2 sin 2 φ cos 2 φ b=2( p sin 2 φ1 )( 1 r 0 R p sin 2 φ )

Dibujamos la trayectoria de un proyectil disparado con un ángulo φ=π/6 (30°) desde una distancia al centro de la Tierra r0=1.2R, con una velocidad v0 o parámetro p=0.6,

R=1; %radio de la Tierra
phi=pi/6;  %ángulo de tiro

r0=1.2*R; %distancia al centro
p=0.6; %parámetro
a=(p*sin(phi)^2-1)^2+(p*sin(2*phi))^2/4;
b=2*(p*sin(phi)^2-1)*(1-r0*p*sin(phi)^2/R);
c=(1-r0*p*sin(phi)^2/R)^2-(p*sin(2*phi))^2/4;
%alcance
alfa=real(acos((-b-sqrt(b^2-4*a*c))/(2*a)));
if phi>pi/2
    alfa=real(acos((-b+sqrt(b^2-4*a*c))/(2*a)));
end
hold on
ang=(1:360)*pi/180; 
fill(cos(ang),sin(ang),'c') %la Tierra
%trayectoria
r=@(x) r0./(sin(phi-x)/sin(phi)+(1-cos(x))/(p*sin(phi)^2));
fplot(@(x) r(x).*cos(x)/R, @(x) r(x).*sin(x)/R,[0,real(alfa)])
line([0,r0/R],[0,0],'lineStyle','--')
plot(r0/R,0,'ro','markersize',3,'markeredgecolor','r','markerfacecolor','r')
line([0,cos(alfa)],[0,sin(alfa)],'lineStyle','--')
plot(cos(alfa),sin(alfa),'ro','markersize',3,'markeredgecolor',
'r','markerfacecolor','r')
quiver(r0,0,sqrt(p/r0)*cos(phi),sqrt(p/r0)*sin(phi));
hold off
axis equal
xlabel('x')
ylabel('y')
title('Disparo de un proyectil')

El alcance en grados es

>> alfa*180/pi
ans =   39.6910

Alcance máximo

El alcance α máximo se obtiene para el ángulo φ tal que dα/dφ=0

{ 1 sin 2 φ ( 1 p sinα dα dφ +( 1 dα dφ )cos( φα )sinφ+sin( φα )cosφ ) 2 cosφ sin 3 φ ( 1cosα p +sin( φα )sinφ )=0 } ( 1 p sinαcos( φα )sinφ ) dα dφ =2 cosφ sinφ ( 1cosα p )sinα dα dφ = 2( 1cosα ) / tanφpsinα sinαpsinφcos( φα )

Igulando el numerador a cero

2( 1cos α m ) tan φ m psin α m =0 tan φ m = 2 p tan( α m 2 )

Sustituimos tan(αm/2) en la ecuación del alcance

r 0 R = sin φ m cos α m cos φ m sin α m sin φ m + 1cos α m p sin 2 φ m

Teniendo en cuenta las relaciones trigonométricas

sinθ= tanθ 1+ tan 2 θ cosθ= 1 1+ tan 2 θ sinθ= 2tan( θ 2 ) 1+ tan 2 ( θ 2 ) cosθ= 1 tan 2 ( θ 2 ) 1+ tan 2 ( θ 2 )

Se convierte en la expresión

r 0 R = ( 1 p 2 )+ tan 2 ( α m 2 )( 2 p 1 ) 1+ tan 2 ( α m 2 ) tan 2 ( α m 2 )= r 0 R 1+ p 2 2 p 1 r 0 R

El ángulo óptimo de disparo φm para el cual el alcance es máximo αm es

tan φ m = 2 p tan( α m 2 )= 2 p r 0 R 1+ p 2 2 p 1 r 0 R

Dados la distancia al centro r0 y el parámetro p (relacionado con la velocidad de disparo v0), calculamos el ángulo de tiro óptimo φm para el cual el alcance es máximo αm. Representamos la trayectoria para este ángulo de tiro y para otros dos ángulos de tiro, mayor y menor en 10°, comprobando que el alcance de estos dos últimos es más pequeño

R=1; %radio de la Tierra
r0=1.2*R; %distancia al centro
p=0.8;
%ángulo de tiro para el que el alcance sea máximo
phi_m=atan(2*sqrt((r0/R-1+p/2)/(2/p-1-r0/R))/p);
if phi_m<0
    phi_m=phi_m+pi;
end
hold on
ang=(1:360)*pi/180; %tierra
fill(cos(ang),sin(ang),'c') %la Tierra
for phi=[phi_m-10*pi/180,phi_m,phi_m+10*pi/180]
    a=(p*sin(phi)^2-1)^2+(p*sin(2*phi))^2/4;
    b=2*(p*sin(phi)^2-1)*(1-r0*p*sin(phi)^2/R);
    c=(1-r0*p*sin(phi)^2/R)^2-(p*sin(2*phi))^2/4;
%alcance
    alfa=real(acos((-b-sqrt(b^2-4*a*c))/(2*a)));
    if phi>pi/2
        alfa=real(acos((-b+sqrt(b^2-4*a*c))/(2*a)));
    end
    disp(alfa*180/pi);
%trayectoria
    r=@(x) r0./(sin(phi-x)/sin(phi)+(1-cos(x))/(p*sin(phi)^2));
    fplot(@(x) r(x).*cos(x)/R, @(x) r(x).*sin(x)/R,[0,real(alfa)])
    line([0,cos(alfa)],[0,sin(alfa)],'lineStyle','--')
    plot(cos(alfa),sin(alfa),'ro','markersize',3,'markeredgecolor',
'r','markerfacecolor','r')
end
line([0,r0/R],[0,0],'lineStyle','--')
plot(r0/R,0,'ro','markersize',3,'markeredgecolor','r','markerfacecolor','r')
hold off
axis equal
xlabel('x')
ylabel('y')
title('Disparo de un proyectil')

El ángulo de tiro óptimo es φm=74.2°, los alcances para los ángulos de tiro: 64.2°,74.2°y 84.2° son, respectivamente

103.6312
109.4712
98.2124
>> phi_m*180/pi
ans =   74.2068

El alcance máximo posible es αm=π, como tan(π/2)=∞, el denominador se ha de hacer cero. Cuando el alcance es αm=π, el ángulo de tiro φm=π/2

2 p 1 r 0 R =0 p= 2 1+ r 0 R

Sea el ángulo de tiro φ=π/2, asignamos al parámetro p el valor deducido para que el alcance del proyectil sea el máximo posible α

R=1; %radio de la Tierra
phi=pi/2;  %ángulo de tiro

r0=1.2*R; %distancia al centro
p=2/(1+r0/R); %parámetro
a=(p*sin(phi)^2-1)^2+(p*sin(2*phi))^2/4;
b=2*(p*sin(phi)^2-1)*(1-r0*p*sin(phi)^2/R);
c=(1-r0*p*sin(phi)^2/R)^2-(p*sin(2*phi))^2/4;
%alcance
alfa=real(acos((-b-sqrt(b^2-4*a*c))/(2*a)));
if phi>pi/2
    alfa=real(acos((-b+sqrt(b^2-4*a*c))/(2*a)));
end
hold on
ang=(1:360)*pi/180;
fill(cos(ang),sin(ang),'c') %la Tierra
%trayectoria
r=@(x) r0./(sin(phi-x)/sin(phi)+(1-cos(x))/(p*sin(phi)^2));
fplot(@(x) r(x).*cos(x)/R, @(x) r(x).*sin(x)/R,[0,real(alfa)])
line([0,r0/R],[0,0],'lineStyle','--')
plot(r0/R,0,'ro','markersize',3,'markeredgecolor','r','markerfacecolor','r')
line([0,cos(alfa)],[0,sin(alfa)],'lineStyle','--')
plot(cos(alfa),sin(alfa),'ro','markersize',3,'markeredgecolor',
'r','markerfacecolor','r')
quiver(r0,0,sqrt(p/r0)*cos(phi),sqrt(p/r0)*sin(phi));
hold off
axis equal
xlabel('x')
ylabel('y')
title('Disparo de un proyectil')

El alcance en grados y el valor del parámetro p son

>> alfa*180/pi
ans =   180.0000
>> p
p =    0.9091

Referencias

Leon Blitzer. Maximun Range of a Projectile in Vacuum on a Spherical Earth. Am. J. Phys. 25 1957, pp. 21-24