# Задача о бросании камешка в вакууме, гл.11 "Задача о бросании камешка: численное решение".
#
# Эта программа рассчитывает и показывает траекторию падения камешка в невязкой среде (в вакууме).
# Для вычисления траектории (интегрирования уравнений движения) используется метод Эйлера 1-го порядка.
# Для отображения траектории используется библиотека matplotlib.
#

from math import pi,sin,cos

#== Параметры задачи ==
Vo = 20             # начальная скорость броска
angle_deg = 86.149  # угол броска
dt = 1e-4           # шаг по времени, 1e-3 = 0.001, 1e-4 = 0.0001

g = 10              # ускорение свободного падения
Xo = 0              # координаты точки броска
Yo = 0              #   --//--
L =  5              # дистанция до окна
h =  5              # высота до окна
hw = 1              # высота самого окна

#== Изменяемые переменные ==
Ax = 0              # ускорение
Ay = -g

phi = angle_deg*pi/180
Vx = Vo*cos(phi)    # скорость
Vy = Vo*sin(phi)

x = Xo              # координаты
y = Yo

Xsav = []           # списки сохранённых координат
Ysav = []           #   - для отрисовки траектории

#===============================
#   Rate(Ax,Ay, Vx,Vy, x,y, dt)
# Расчёт скорости изменения величин (уравнения движения):
#   V' = A      ( dV = A*dt )
#   R' = V      ( dR = V*dt )
#===============================
def Rate(Ax,Ay, Vx,Vy, x,y, dt):
    dAx = 0
    dAy = 0
    dVx = Ax
    dVy = Ay
    dx  = Vx
    dy  = Vy
    return dAx,dAy, dVx, dVy, dx, dy 
#===============================

FMT = """
В окно:      {}
   высота:   {} м
Скорость:    {} м/сек
Угол броска: {}°
"""

saveKth = int((2*Vo/(g*dt))/1000)   # saveKth = сколько точек можно пропускать, 
                                    #     чтобы сохранить ~1 тыс.точек
k = 0
time = 0

if Vx*dt *10**5<L:
    print("горизонтальная скорость слишком мала!")
else:
    while x<L:
        if k%saveKth==0:
            Xsav.append(x)
            Ysav.append(y)
        
        dAx,dAy, dVx,dVy, dx,dy = Rate(Ax,Ay, Vx,Vy, x,y, dt)

        # расчёт величин на следующем шаге по времени
        Ax,Ay, Vx,Vy, x,y = Ax+dAx*dt,Ay+dAy*dt, Vx+dVx*dt,Vy+dVy*dt, x+dx*dt, y+dy*dt
        time+=dt
        k+=1

    if h<=y<=h+hw:
        print(FMT.format("попали!", y, Vx, angle_deg))
    else:
        print(FMT.format("не попали", y, Vx, angle_deg))

#===============================
#  Визуализация траектории
#===============================
import matplotlib.pyplot as PLT
PLT.plot(Xsav, Ysav, label="Траектория тела") # рисуем траекторию тела
PLT.plot([L,L], [h,h+hw], label="окно")       # рисуем окно
PLT.legend()                                  # показывать легенду - описание графиков

PLT.show()                                    # показать окно matplotlib