import numpy as np

# object
A    = .8               # surface area m^2
C    = 1.4              # drag coefficient
m_prop = 65 # mass of object in kg

m = 90                  # number of angles
n = 1000                # number of time steps to simulate
rho = 1.225             # kg/m^3
g = -9.8                # m/s^2
start = 0.0             # seconds
end = 10                # seconds
initial_velocity = 1500 * ((m_prop / 65) ** .45)  # meters/second
initial_height = 5


# parameters of simulation
angles = np.linspace( 1,m,m )
radians = np.radians( angles )

#initalize variables
t  = np.linspace( start,end,n+1 )
vx = np.zeros( ( m,n+1 ) )
vy = np.zeros( ( m,n+1 ) )
x  = np.zeros( ( m,n+1 ) )
y  = np.zeros( ( m,n+1 ) )

for i in range(m):
    y[ i,0 ]  = initial_height
    vx[ i,0 ] = initial_velocity * np.cos( radians[ i ] )
    vy[ i,0 ] = initial_velocity * np.sin( radians[ i ] )

for i in range( m ):
    for j in range( 1,n+1 ):
        if y[ i,j-1 ] <= 0:
            x[ i,j ]  = x[ i,j-1 ]
            y[ i,j ]  = 0
            vx[ i,j ] = 0
            vy[ i,j ] = 0
        else:
            v  = np.sqrt( vx[ i,j-1 ] ** 2 + vy[ i,j-1 ] ** 2 )
            ax =   - ( 0.5*rho*C*A/m_prop ) * v**2 * ( vx[ i,j-1 ] / v )
            ay = g - ( 0.5*rho*C*A/m_prop ) * v**2 * ( vy[ i,j-1 ] / v )
            # https://www.grc.nasa.gov/www/k-12/airplane/flteqs.html
            dt        = t[ j ] - t[ j-1 ]
            vx[ i,j ] = vx[ i,j-1 ] + ax*dt
            vy[ i,j ] = vy[ i,j-1 ] + ay*dt
            x[ i,j ]  =  x[ i,j-1 ] + vx[ i,j ]*dt
            y[ i,j ]  =  y[ i,j-1 ] + vy[ i,j ]*dt
            if y[ i,j ] <= 0:
                y[ i,j ]  = 0
                x[ i,j ]  = x[ i,j-1 ]
                vy[ i,j ] = 0
                vx[ i,j ] = 0

best_d = 0
best_a = 0
for i in range( m ):
    if x[ i,-1 ] > best_d:
        best_d = x[ i,-1 ]
        best_a = angles[ i ]

def dist( angle,m_prop ):
    pass