Fixing projectile motion formulas in python

The code outputs a graph that is nowhere near what you'd expect of a projectile motion type of graph. In addition, if you change the step number to around 2, the graph doesn't output much of anything.

import numpy as np
import matplotlib.pyplot as plt

grav = 9.8
airDen = 1.225      # kg/m^3
areaProj = 0.8643   # m^2
v0 = 198.1          # m/s
angle = 60.0   
dragCoef = 0.62     #Approximate value
mass = 15.0         #kg
step = 0.02    #time step in seconds 

t = [0]
vx = [v0*np.cos(angle*np.pi/180)]
vy = [v0*np.sin(angle*np.pi/180)]
x = [0]
y = [0]

dragForce = 0.5*dragCoef*areaProj*(v0**2)*airDen


ax = [-(dragForce*np.cos(angle/180*np.pi))/mass]
ay = [-grav-(dragForce*np.sin(angle/180*np.pi)/mass)]


counter = 0
while(y[counter] >= 0):
    t.append(t[counter]+step)

    vx.append(vx[counter]+step*ax[counter])
    vy.append(vy[counter]+step*ay[counter])

    x.append(x[counter]+step*vx[counter])
    y.append(y[counter]+step*vy[counter])

    vel = np.sqrt(vx[counter+1]**2 + vy[counter+1]**2)
    dragForce = 0.5*dragCoef*areaProj*(vel**2)*airDen

    ax.append(-(dragForce*np.cos(angle/180*np.pi))/mass)
    ay.append(-grav-(dragForce*np.sin(angle/180*np.pi)/mass))

    counter=counter+1

plt.plot(x,y)
plt.ylabel("y (m)")
plt.xlabel("x (m)")

Graph-example

Solution

Looks like you're not updating your angle, so instead of reaching zero vx it accelerates backwards since it thinks you're still going at 60 degrees initial angle. Adding an updated angle based on the current velocity vector results in the following:

import numpy as np
import matplotlib.pyplot as plt

grav = 9.8
airDen = 1.225      # kg/m^3
areaProj = 0.8643   # m^2
v0 = 198.1          # m/s
angle = 60.0   
dragCoef = 0.62     #Approximate value
mass = 15.0         #kg
step = 0.02    #time step in seconds 

t = [0]
vx = [v0*np.cos(angle*np.pi/180)]
vy = [v0*np.sin(angle*np.pi/180)]
x = [0]
y = [0]

dragForce = 0.5*dragCoef*areaProj*(v0**2)*airDen


ax = [-(dragForce*np.cos(angle/180*np.pi))/mass]
ay = [-grav-(dragForce*np.sin(angle/180*np.pi)/mass)]


counter = 0
while(y[counter] >= 0):
    t.append(t[counter]+step)

    vx.append(vx[counter]+step*ax[counter])
    vy.append(vy[counter]+step*ay[counter])

    x.append(x[counter]+step*vx[counter])
    y.append(y[counter]+step*vy[counter])

    vel = np.sqrt(vx[counter+1]**2 + vy[counter+1]**2)
    dragForce = 0.5*dragCoef*areaProj*(vel**2)*airDen

    angle = np.arctan(vy[counter]/vx[counter]) * (180 / np.pi)

    ax.append(-(dragForce*np.cos(angle/180*np.pi))/mass)
    ay.append(-grav-(dragForce*np.sin(angle/180*np.pi)/mass))

    counter=counter+1


plt.plot(x,y)
plt.ylabel("y (m)")
plt.xlabel("x (m)")