# pallon lento, mukana ilmanvastus riippuvuuksineen 
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# Pallon lennon parametrit (otettu  baseballista)
g       = 9.80665         # Painovoiman putouskiihtyvyys (m/s^2)
rho = 1.225         # Ilman tiheys (kg/m^3), ICAO, 0m +15C, 0m korkeus
m       = 0.145        # Pallon massa (kg)
Cd      = 0.3          # Pallon muotokerroin
A = 0.004185 # pallon poikkinleikkauksen pinta-ala m²
k1 = (Cd*A)/(2*m) #ilmanvastuskertoimen vakio-osa

# Lähtöpaikka, lähtönopeuden komponentit
x0, y0 = 0.0, 0.0 # paikan alkuarvot
v0      = 300        # Lähtönopeus (m/s)
alfa = np.radians(30) # lähtökulma
vx0, vy0 = v0 * np.cos(alfa), v0 * np.sin(alfa) # nopeuskomponentit x,y-suunnassa

def model(U,t):
    x,y,vx,vy = U # mallin tilamuuttujat, joiden derivaatat palautetaan
    v = np.sqrt(vx**2+vy**2) # nopeus
    k = rho*(288.15/(288.15-0.0065*y))**(-4.25578597)*k1 # ilmanvastuskerroin, mukana lämpötila
    dx = vx # x' eli x:n derivaatta, vaakanopeus
    dy = vy # y' eli y:n derivaatta, pystynopeus
    dvx = -k*v*vx # kiihtyvyys x-suunnassa, vx/dt
    dvy = -g - k*v*vy # kiihtyvyys y-suunnassa, vy/dt
    # mukana ilmanpaineen lampötilariippuvuus
    # palauttaa derivaatat vx,vy,vx/dt,vy/dt
    return[dx,dy,dvx,dvy]              
U0 = [0,0,vx0,vy0] # paikan ja lähtönopeuden alkuarvot
ts = np.linspace(0,11,1000) # laskennan aikaväli 0-11s, 1000 näytettä
Us = odeint(model, U0, ts)
# Pallon lentorata graafisesti
x = Us[:,0]
y = Us[:,1]
vx = Us[:,2]
vy = Us[:,3]
maksimikorkeus = np.amax(y)
print("maksimikorkeus = ",maksimikorkeus)
maksimikohta = np.where(y>=maksimikorkeus)
print("lentoaika maksimikorkeuteen= ",ts[maksimikohta][0])
print("vaakanopeus maksimikorkeudessa = ",vx[maksimikohta][0])
maahantulokohta = np.where(y<0)
print("kantama = ",x[maahantulokohta][0], " m")
print("lentoaika = ",ts[maahantulokohta][0]," s")
vxmaa = Us[:,2][maahantulokohta][0]
vymaa = Us[:,3][maahantulokohta][0]
print("tulokulma = ",np.degrees(np.arctan(abs(vymaa/vxmaa))))
fig, ax = plt.subplots()
ax.plot(x, y, "r-", label="Pallon lentorata")
ax.set_title(r"rho(0m,15C)=1.225, m=0.145 kg,r=0.0365 m, alfa=30,v0=300 m/s")
ax.set_aspect("equal")
ax.grid(b=True)
ax.legend()
ax.set_xlabel("$x$ (m)")
ax.set_ylabel("$y$ (m)")
plt.show()

#plt.plot(x,y)
#plt.show()