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

from math import pi,sin,cos

#===============================
#   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):
    dAx = 0
    dAy = 0
    dVx = dt*Ax
    dVy = dt*Ay
    dx  = dt*Vx
    dy  = dt*Vy
    return Ax+dAx,Ay+dAy, Vx+dVx, Vy+dVy, x+dx, y+dy 
#===============================

#===============================
# Моделирование движения
#===============================
#   Solve(angle0)
# Для заданного угла броска angle0 рассчитывается
# горизонтальная скорость камня при попадании в стену.
# Если камень не попал в окно, то к скорости прибавляется "штраф".
#
def Solve(angle0):
    #== Параметры задачи ==
    Vo = 20             # начальная скорость броска
    angle_deg = angle0  # угол броска
    dt = 1e-3           # шаг по времени, 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

    k = 0
    time = 0

    if Vx*dt *10**5<L:
        return Vx+L/dt
    else:
        while x<L:
            Ax,Ay, Vx,Vy, x,y = Next(Ax,Ay, Vx,Vy, x,y, dt)
            time+=dt
            k+=1

        if y<h:
            return Vx+(h-y)*0.33
        elif y>h+hw:
            return Vx+(y-(h+hw))*0.33
        else:
            return Vx

#===============================
# Поиск оптимального угла броска
#   - Используем пакет SciPy.optimize
#===============================
#
import numpy as NP
from scipy.optimize import minimize, Bounds

x0 = NP.array( [60.0] )     # задаём начальные значения
            # вызываем функцию минимизации
res = minimize(Solve, x0, method='nelder-mead', options={'xatol': 1e-12, 'disp': False})

#---
#   Другой вариант вызова функции - здесь задаются границы области поиска:
#res = minimize(Solve, x0=80.0, method='nelder-mead', options = {"bounds": Bounds(lb=10, ub=89) } )
#---

    # анализ и печать ответа
if res.success:
    print("[+] Solution found")
    print("optimal angle = {}°, Vx = {} m/s".format(res.x, res.fun))
else:
    print("(!) Failed find a solution")
print("\n\n===\nExhaustive info:\n",res)