from math import pi,sin,cos,sqrt

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

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

# параметры бросаемого предмета (шарик для пинг-понга)
radii = 0.02        # радиус [м]
mass= 0.0027        # масса  [кг]
Cf = 0.47           # коэффициент формы для сферы
air_viscosity = 1.83e-5 # коэф.вязк.возд.
air_density = 1.293 # плотность возд.(при н.у.)

k1 = 0.002555162    #6.89894E-06/mass = 6*pi*radii*air_viscosity
k2 = 0.141420539    #0.000381835/mass = Cf*pi*radii**2/2*air_density

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

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

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

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

#===============================
#   Next(Ax,Ay, Vx,Vy, x,y, dt)
# Расчёт скорости изменения величин (уравнения движения):
#   V' = A      ( dV = A*dt )
#   R' = V      ( dR = V*dt )
#===============================
def Next(Ax,Ay, Vx,Vy, x,y, dt):
    V2= Vx**2+Vy**2
    V = sqrt(V2)
    Ax1 =  0 -Vx*(k1+k2*V)
    Ay1 = -g -Vy*(k1+k2*V)
    Vx1 = Vx+dt*Ax
    Vy1 = Vy+dt*Ay
    x1  = x+dt*Vx
    y1  = y+dt*Vy
    return Ax1,Ay1, Vx1,Vy1, x1,y1 
#===============================

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 and y>=-1:
        if k%saveKth==0:
            Xsav.append(x)
            Ysav.append(y)
        
        # расчёт величин на следующем шаге по времени
        Ax,Ay, Vx,Vy, x,y = Next(Ax,Ay, Vx,Vy, x,y, 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