import numpy as np

def forza_generalizzata(y,m,g, Catt):
    f=np.zeros(4)
    f[0]=y[2]
    f[1]=y[3]
    vel=np.sqrt(y[2]**2+y[3]**2)
    u_x=y[2]/vel
    u_y=y[3]/vel
    f[2]=-u_x*Catt*vel**2/m
    f[3]=-u_y*Catt*vel**2/m-g
    return f


def main():
    g=9.81
    rho_aria=1.22
    rho_gran=2700.0
    R=0.1
    Cd=0.47


    A=np.pi*R**2
    Catt=0.5*rho_aria*Cd*A
    m=rho_gran*(4/3)*np.pi*R**3


    theta=float(input('Alzo (gradi) :'))
    theta*=np.pi/180
    v0=float(input('Velocità iniziale (m/s)'))
    dt=float(input('Time step (s): '))
    
    it=0

    y0=np.zeros(4)
    y1=np.zeros(4)

    y0[2]=np.cos(theta)*v0
    y0[3]=np.sin(theta)*v0

    found=False

    it=0
    while(not found):
        
        y1=y0+forza_generalizzata(y0,m,g, Catt)*dt
        if(y0[1]*y1[1]<0):
            gitt=y0[0]-(y1[0]-y0[0])/(y1[1]-y0[1])*y0[1]
            found=True
        it+=1
        y0=y1
       


    print('Gittata (m):' , gitt)
    gitt_ide=2*v0**2/g*np.sin(theta)*np.cos(theta)
    print('Gittata ideale (m):', gitt_ide)
    











if __name__ == "__main__":
    main()