В этой статье мы завершим решение задачи внешней баллистики разбором шестого и седьмого случаев. В них мы учтём уменьшение гравитации с высотой, а также кривизну Земли.

Реализуем всё вышенапечатанное на Python:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# Константы
g0 = 9.80665  # м/с² — ускорение свободного падения на уровне моря
RE = 6371000  # м — радиус Земли
rho0 = 1.225  # кг/м³ — плотность воздуха на уровне моря
H = 8000  # м — масштабная высота

# Скорость звука как функция высоты (упрощённая модель)
def sound_speed(h):
    if h < 11000:  # Тропосфера
        return 340 - 0.0065 * h
    else:  # Стратосфера и выше
        return 295

# Параметры снаряда
m = 2000  # кг — масса
S = 0.3  # м² — площадь миделева сечения

# Начальная скорость и угол
v0 = 4000  # м/с
theta0 = 50  # градусов

def cx(M):
    """Коэффициент лобового сопротивления по закону 1943 года"""
    if M < 0.8:
        return 0.15 + 0.25 * M
    elif M < 1.2:
        return 0.4 - 0.3 * (M - 0.8)
    elif M<4:
        return 0.28 + 0.05 * (M - 1.2)
    else:
        return 1

def rho(h):
    """Плотность воздуха с экспоненциальным убыванием"""
    return rho0 * np.exp(-h / H)

def g(h):
    """Ускорение свободного падения с высотой"""
    return g0 * (RE / (RE + h))**2

def equations(t, state):
    """
    Система дифференциальных уравнений движения
    state = [r, phi, vr, vphi]
    r — расстояние от центра Земли
    phi — полярный угол
    vr — радиальная скорость
    vphi — азимутальная скорость
    """
    r, phi, vr, vphi = state
    h = r - RE  # высота над поверхностью

    # Проверка выхода за пределы расчёта
    if h < -100:  # снаряд ушёл под землю
        return [0, 0, 0, 0]
    if h > 1e6:  # высота более 1000 км — сопротивление пренебрежимо мало
        # Упрощённые уравнения без сопротивления
        dvr_dt = -g(h) + vphi**2 / r
        dvphi_dt = -vr * vphi / r
        return [vr, vphi / r, dvr_dt, dvphi_dt]

    # Скорость и число Маха
    v = np.sqrt(vr**2 + vphi**2)
    if v == 0:
        M = 0
    else:
        a = sound_speed(h)
        M = v / a

    # Плотность воздуха и коэффициент сопротивления
    rho_h = rho(h)
    cx_val = cx(M)

    # Силы сопротивления
    F_drag = 0.5 * rho_h * v**2 * S * cx_val
    Fr = -F_drag * vr / v  # радиальная компонента
    Fphi = -F_drag * vphi / v  # азимутальная компонента

    # Ускорения
    dvr_dt = Fr / m - g(h) + vphi**2 / r
    dvphi_dt = Fphi / m - vr * vphi / r

    return [vr, vphi / r, dvr_dt, dvphi_dt]

def event_ground(t, state):
    """Событие — достижение поверхности Земли"""
    r = state[0]
    return r - RE  # ноль при r = RE
event_ground.terminal = True  # останавливаем интегрирование
event_ground.direction = -1  # только при уменьшении r

# Начальные условия
r0 = RE  # начальное расстояние от центра Земли
phi0 = 0  # начальный угол
vr0 = v0 * np.sin(np.radians(theta0))  # радиальная скорость
vphi0 = v0 * np.cos(np.radians(theta0))  # азимутальная скорость

initial_state = [r0, phi0, vr0, vphi0]

# Временной интервал (максимальное время)
t_span = (0, 1000)  # 0–1000 секунд (достаточно для большинства случаев)

# Решение системы ОДУ с событием
solution = solve_ivp(
    equations,
    t_span,
    initial_state,
    method='RK45',
    events=event_ground,
    dense_output=True
)

# Извлечение результатов
t = solution.t
r = solution.y[0]
phi = solution.y[1]
vr = solution.y[2]
vphi = solution.y[3]

# Преобразование в декартовы координаты для визуализации
x = np.sqrt(vr**2 + vphi**2)*t * np.cos(np.radians(phi))

# Расчёт дополнительных параметров
h = r - RE
v_total = np.sqrt(vr**2 + vphi**2)
M_history = [v_total[i] / sound_speed(h[i]) if v_total[i] > 0 else 0 for i in range(len(v_total))]

# Поиск точки падения (из события)
if solution.t_events[0].size > 0:
    fall_time = solution.t_events[0][0]
    # Интерполяция для точного определения координат в момент падения
    state_fall = solution.sol(fall_time)
    r_fall = state_fall[0]
    phi_fall = state_fall[1]
    range_m = (r_fall - RE) * phi_fall  # приближённо для малых углов
    max_height = np.max(h)  # максимальная высота
else:
    fall_time = t[-1]
    range_m = x[-1]
    max_height = np.max(h)

# Визуализация
plt.figure(figsize=(14, 10))

# Траектория
plt.subplot(2, 2, 1)
plt.plot(x/1000,h/1000, label='Траектория', linewidth=2)
plt.xlabel('X, км')
plt.ylabel('Y, км')
plt.title('Траектория полёта')
plt.grid(True, alpha=0.7)
plt.axis('equal')
plt.legend()

# Высота во времени
plt.subplot(2, 2, 2)
plt.plot(t, h/1000, label='Высота', linewidth=2, color='orange')
plt.xlabel('Время, с')
plt.ylabel('Высота, км')
plt.title('Высота полёта')
plt.grid(True, alpha=0.3)
plt.legend()

# Скорость во времени
plt.subplot(2, 2, 3)
plt.plot(t, v_total, label='Скорость', linewidth=2, color='green')
plt.xlabel('Время, с')
plt.ylabel('Скорость, м/с')
plt.title('Скорость полёта')
plt.grid(True, alpha=0.3)
plt.legend()

# Число Маха во времени
plt.subplot(2, 2, 4)
plt.plot(t, M_history, label='Число Маха', linewidth=2, color='red')
plt.xlabel('Время, с')
plt.ylabel('Число Маха')
plt.title('Число Маха полёта')
plt.grid(True, alpha=0.3)
plt.legend()

plt.tight_layout()
plt.show()

Результат моделирования можно увидеть ниже:

Таким образом, исходная задача решена, создана максимально реалистичная математическая модель в общем виде, учитывающая очень многие параметры и факторы, влияющие на движение тела в атмосфере Земли, и предложен алгоритм, который её решает за разумное время с большой точностью.

Литература:

Окунев Б. Н. «Решение основной задачи внешней баллистики при квадратичном законе сопротивления воздуха» (1932).
Окунев Б. Н. «Основная задача внешней баллистики и аналитические методы её решения» (1934).
Дмитриевский А. А., Лысенко Л. Н. «Внешняя баллистика» (2005).
Лысенко А. Н. «Внешняя баллистика» (2024).
Шапиро Я. М. «Внешняя баллистика» (1946).
Беляева С. Д. «Внешняя баллистика с примерами и задачами» (2023).
https://docviewer.yandex.ru/view/1034483224/?*=CBQ5Rag1zM4T3JJp60YBzgS9YU17InVybCI6Imh0dHBzOi8vbGlicmFyeS5rdXpzdHUucnUvZGwucGhwP249MTQwNzQ0JnR5cGU9bnN0dTpjb21tb24iLCJ0aXRsZSI6ImRsLnBocD9uPTE0MDc0NCIsIm5vaWZyYW1lIjp0cnVlLCJ1aWQiOiIxMDM0NDgzMjI0IiwidHMiOjE3NzQ3ODM0MTI0NjQsInl1IjoiMjg4OTYzODY3MTY5NzAzMzU0MCIsInNlcnBQYXJhbXMiOiJ0bT0xNzc0NzgzNDAyJnRsZD1ydSZsYW5nPXJ1Jm5hbWU9ZGwucGhwP249MTQwNzQ0JnR5cGU9bnN0dTpjb21tb24mdGV4dD0lRDAlQjclRDAlQjAlRDAlQjQlRDAlQjAlRDElODclRDAlQjArJUQwJUIyJUQwJUJEJUQwJUI1JUQxJTg4JUQwJUJEJUQwJUI1JUQwJUI5KyVEMCVCMSVEMCVCMCVEMCVCQiVEMCVCQiVEMCVCOCVEMSU4MSVEMSU4MiVEMCVCOCVEMCVCQSVEMCVCOCslRDAlQkQlRDAlQjMlRDElODIlRDElODMmdXJsPWh0dHBzJTNBLy9saWJyYXJ5Lmt1enN0dS5ydS9kbC5waHAlM0ZuJTNEMTQwNzQ0JTI2dHlwZSUzRG5zdHUlM0Fjb21tb24mbHI9NTYmbWltZT1wZGYmbDEwbj1ydSZzaWduPTljZWQzZTc4M2MxNmIwODI0YTE3YjNmZDVmODU0OWMyJmtleW5vPTAifQ%3D%3D&lang=ru
Мандрыка А. П. «Баллистические исследования Леонарда Эйлера» (2017)
https://en.wikipedia.org/wiki/Projectile_motion#Time_of_flight_with_air_resistance

Задача внешней баллистики. Третья часть