Projectile motion sim, index out of bounds error, nested for/while loops

Question

I'm trying to graph the optimum angle of launch as a function of the mass and graph it. for each element of mlist I want it to calculate the distance traveled for that m for each element of thetalist, then take the highest distance for that m and get the theta that gave that answer and make a new list of these so called theta_optimum and plot it vs mlist.

Here is my code:

import matplotlib.pylab as plt
import numpy as np
import math
#############
Die = 10.0
B = 1.6e-4
b = B*Die
g = 9.81
tmax = 1000
dt = 0.01
v = 100.0
mlist = np.arange(1.0, 100.0, 1.0)
thetalist = np.arange(0.1, math.pi / 2.0, 0.1)
theta_op = []

for j in range(len(mlist)):
    for i in range(len(thetalist)):
        vy = math.sin(thetalist[i]) * v
        vx = math.cos(thetalist[i]) * v
        t = 0.0
        x = 0.0
        y = 0.0
        xlist = []
        while y >= 0.0:
            vx1 = -b/mlist[j]*vx
            vx = vx + vx1*dt
            dx = vx * dt
            x = x + dx

            vy1 = -g-b/mlist[j]*vy
            vy = vy + vy1 * dt
            dy = vy * dt
            y = y + dy

            t = t + dt

            xlist.append(x)

        theta_op.append(thetalist[xlist.index(max(xlist))])

plt.plot(mlist, theta_op, color='red')
plt.show()

The error I am getting is:

line 39, in <module>
theta_op.append(thetalist[xlist.index(max(xlist))])
IndexError: index 202 is out of bounds for axis 0 with size 15

My question is how do I fix this error in my code and why has it arisen. I know the physics is right but as I am new to coding I am not sure if the way I have made my lists and loops is correct. I have looked for other questions like this but haven't found any.

Answer 1

It seems like you are severely mixing up your variable scopes.

For every value of mass and theta, you save every value of x during the integration procedure. Assuming that's what you wanted to do, you then wrote the line causing this error:

theta_op.append(thetalist[xlist.index(max(xlist))])

The above returns the index of the maximum value in xlist; you then uses this index to access thetalist, which makes no sense at all. The size of thetalist and xlist are independent, the former depending on the theta step, and the latter on the simulation time step; so using the index of one to access the other may well result in this overflow.

This brings us back to the question of why you want to save all x values in the first place. Why not just save the last x value for each value of theta, and get the optimum theta for each mass, which exactly matches your problem description?

for j in range(len(mlist)):

    best_theta = 0.0
    best_x = 0.0

    for i in range(len(thetalist)):

        # initialisation code as before
        # ...

        while y >= 0.0:
            # integration code as before
            # ...

        # update the best theta for this mass
        if x > best_x:
            best_theta = thetalist[i]
            best_x = x

    # add the best theta to theta_op without having to do any searches   
    theta_op.append(best_theta)

---

EDIT: I tested this code with various theta step values, and all results were within reasonable range of the zero-resistance limit (45 degrees = 0.785398... radians), so you may wish to consider increasing the drag coefficient.