Projectile Motion (With Air Resistance) In Python

Simulating a projectile motion

By MikeRun — Own work, using, CC BY-SA 4.0, https://commons.wikimedia.org/w/index.php?curid=120407200

I tried to write a python code that simulates a projectile motion and returns the time of flight — considering air resistance too.

It includes all the parameters and variables: acceleration due to gravity, launch angle, launch height, mass of projectile and drag constant.

Changing conditions lead to different projectile paths.

I have provided the full code in the end. First, let me explain the code in detail.

The air resistance model I am using is

equation of air resistance

The equation f = kv² is a simplified model of air resistance, and is most applicable at high speeds where air resistance becomes significant.

The constant “K” is the drag coefficient, which depends on shape, size of the object and the properties of air.

This function calculates the acceleration components after considering air resistance.

def calculating_accel(t, state, k, m):
 x, y, vx, vy = state
 v = np.sqrt(vx**2 + vy**2)
 F_air_x = -k * v * vx
 F_air_y = -k * v * vy
 ax = F_air_x / m
 ay = (F_air_y - m * g) / m
 return [vx, vy, ax, ay]

Defining the Runge-Kutta Step Function:

The Runge-Kutta method is a numerical integration technique, and here it is used to approximate the displacement of an object from its velocity.

def runge_kutta_step(t, state, dt, k, m):
    k1 = calculating_Acclel(t, state, k, m)
    k2 = calculating_Acclel(t + 0.5*dt, [s + 0.5*dt*k for s, k in zip(state, k1)], k, m)
    k3 = calculating_Acclel(t + 0.5*dt, [s + 0.5*dt*k for s, k in zip(state, k2)], k, m)
    k4 = calculating_Acclel(t + dt, [s + dt*k for s, k in zip(state, k3)], k, m)
  
    return [s + dt/6 * (k1_i + 2*k2_i + 2*k3_i + k4_i) for s, k1_i, k2_i, k3_i, k4_i in zip(state, k1, k2, k3, k4)]

The fourth-order Runge-Kutta method (RK4) is popular because it strikes a good balance between accuracy and computational efficiency.

1. At the beginning of each time step, the main function calls the calculating_Accel function to obtain the current acceleration (ax and ay) at time ‘t’.

The t variable helps in calculating the intermediate states at each time step. (t is not a constant here)

2. It then estimates how the object’s position and velocity will change over time step based on the current acceleration. For this, four intermediate estimates (k1, k2, k3, k4) are calculated to refine the update for each variable (position and velocity) using weighted averages.

Projectile Motion with Air Resistance Function:

This function simulates the projectile motion of the object with air resistance considered.

def projectile_motion_with_air_resistance(v0, theta, h0, k, m):
    if k==0 and theta==90:
        vx0=0
        vy0=v0
    else:
        theta = np.radians(theta)
        vx0 = v0 * np.cos(theta)
        vy0 = v0 * np.sin(theta)
    state = [0, h0, vx0, vy0]
    dt = 0.01
    times = [0]
# This is the initial height of projectile. You can change it
    x_values = [0]
    y_values = [h0]
    
    while state[1] >= 0:
        t = times[-1]
        state = runge_kutta_step(t, state, dt, k, m)
        times.append(t + dt)
        x_values.append(state[0])
        y_values.append(state[1])
    
    time_of_flight = times[-1]
    
    return time_of_flight, times, x_values, y_values

• It initializes the state with [x, y, vx, vy] and sets up the time step (dt) and lists to store time, horizontal position (x_values), and vertical position (y_values).
• The simulation runs inside a while loop until the vertical position (y) of the object becomes less than zero (indicating it has hit the ground).
• The runge_kutta_step function updates the object’s state (position and velocity) to its new values. The new position and velocity are then used to calculate the acceleration at the next time step, and the process continues until object hits the ground.
• The time, horizontal position (x_values), and vertical position (y_values) are stored in lists as the simulation progresses.
• t is used implicitly in the loop to update the state of the object at each time step. The loop keeps running until the object hits the ground (state[1] >= 0).

Setting parameters and initial conditions:

g=9.81
v0 = 40   # Initial velocity (m/s)
launch_angle = 30  
initial_height = 20  
k = 1.225 # Air resistance constant 
#depending on the projectile and environment)
mass = 1  #mass of projectile

The code plots the trajectory of the projectile and flight time:

time_of_flight, times, x_values, y_values = \
projectile_motion_with_air_resistance(v0, launch_angle,initial_height, k, mass)

plt.plot(x_values, y_values)
plt.xlabel('Horizontal Distance (m)')
plt.ylabel('Vertical Distance (m)')
plt.title('Projectile Motion with Air Resistance')
plt.ylim(0,max(y_values)+7)
plt.plot([0,0],[0,initial_height],'k',linewidth=3)
plt.grid(True)
plt.show()

print(f'Time of flight : {time_of_flight:.2f} seconds')

The last line, tells the system to format the numerical value as a floating-point number and display it with two decimal places.

Here is the complete code. I ran it on Jupyter Notebook.

GitHub - hamxa708/Projectile_air_resistance: This code in python simulates a projectile motion with…

I would really appreciate any feedback.

Disclaimer:

This simulation is a basic model and may not accurately represent real-world projectile motion in all scenarios. Depending on your specific case, you may need to consider more complex models.