'''Programma per studiare il moto di un proiettile'''


import numpy as np
import scipy as sp
import matplotlib.pyplot as plt


class air(object):
    def __init__(self, rho, cw, area):
        self.rho=rho
        self.cw=cw
        self.g=9.81
        self.area=area

    def forza_no_attrito(self,y):
        f=np.zeros((4),dtype=np.float64)
        f[0]=y[2]
        f[1]=y[3]
        f[2]=0
        f[3]=-self.g
        return f
        

        return f
        
        

def leap_frog(y0,y1,f,dt):
    y2=y0+2.0*f(y1)*dt 
    return y2

def eulero_esplicito(y0,f,dt):
    y2=y0+f(y0)*dt 
    return y2

def forza_generalizzata(y):
    g=9.81
    f=np.zeros((4),dtype=np.float64)
    f[0]=y[2]
    f[1]=y[3]
    f[2]=0
    f[3]=-g
    return f

def main():

    g=9.81
    
    vini=input("Velocità iniziale (m/s):\t")
    vini=float(vini)

    deg=input("Alzo (gradi):\t")
    deg=float(deg)

    dt=input("Time step (s):\t")
    dt=float(dt)

    cw=input("Coefficiente aereodinamico Cw:\t")
    cw=float(cw)

    raggio=input("Raggio del proiettile (m):\t")
    raggio=float(raggio)

    iprop=input("Propagatore: 1-Eulero esplicito 2-Leap-Frog 3-Runge Kutta: ")
    iprop=int(iprop)

    angle=np.pi*deg/180.

    area=4*np.pi*raggio**2
    aria=air(1.3,cw,area)

    y0=np.zeros((4),dtype=np.float64)
    y0[2]=np.cos(angle)*vini
    y0[3]=np.sin(angle)*vini
   
    it=0
    arrivato=False
    while(not arrivato):
        if(iprop==1 or (iprop==2 and it==0)):
            y1=eulero_esplicito(y0,aria.forza_no_attrito,dt)
        elif(iprop==2):
            y1=leap_frog(ym1,y0,aria.forza_no_attrito,dt)
        elif(iprop==3):
            pass
        if(y1[1]*y0[1]<0):
            arrivato=True
        else:
            ym1=y0
            y0=y1
            it+=1
         

    gittata=(y1[1]*y0[0]-y0[1]*y1[0])/(y1[1]-y0[1])

    gittata_ideale=(2*vini**2/g)*np.sin(angle)*np.cos(angle)

    print("Gittata (m):\t"+str(gittata))
    print("Gittata ideale (m):\t"+str(gittata_ideale))

if __name__ == "__main__":
    main()