How to adjust my basketball-FT(Free throw) simulation to find the necessary torque(s) on joints

I'm working on a simulation of a basketball throw which starts with a angular motion and passes into a projectile motion. Physics set-up of simulation

My goal is to get insights on the amount of torque that is applied on joints(like elbow & shoulder) while throwing a basketball. In my simulation the torques and release angle are the inputs and a trajectory is the output. I want to 'tweek' the input in order to get a ball trajectory that swishes the net (scores). For now I scoped it down to only the elbow joint so basically a catapult as can been seen in the code below, and also the basket isn't in there yet.

But before this expanding I want to run the simulation for multiple 'release angles' and different 'torques on elbow'. As you can see I tried to make a additional list with torques and release angles in line 11 and 17 (commented out) and I wanted to add another loop so the whole simulation would run with the different Angles[re_angl_deg] list and Torques[torq_elb] list as input. Unfortunately this didn't really work. Now my question is: are there other ways to let the simulation run multiple times, with each a different Angle and Torque(like the lists I made).

I hope someone could give me some advise! Regards.TF

import numpy as np
import matplotlib.pyplot as plt
import math

# Define settings.  
endTime = 60 # The time (seconds) that we simulate.  
dt = 0.001 # The time step (seconds) that we use in the discretization.  
w0 = 0 # The initial velocity [m/s].  
tet0 = 15 # The initial position [m].  
torq_elb = 70     #[Nm]
#torq_elb = np.array([50,55,60,65,70,75,80,85,90,95,100])       #[Nm]
I = 0.16        #[kgm^2]  inertia
x_pos = 0       #initial x position
y_pos = 0       # initial y Position
r = 0.3         # lower arm length [m]
re_angl_deg = 50  # release angle degrees
#re_angl_deg = np.array([30,32,34,36,38,40,42,44,45,47,49,50,52,54])
re_angl_rad = math.radians(re_angl_deg)
g = 9.81 # The gravitational acceleration [m/s^2].    
m_ball = 0.6 # The mass of the ball [kg]. 
rho_air = 1.2041 #[kg/m3]
c_drag =  0.3       # drag coefficient
area_b = (4*math.pi*r**2)/2 # surface area of the ball

 
# Set up variables.  
time = np.arange(0, endTime + dt, dt) # A list with all times we want to plot at.  
w = np.zeros(len(time)) # A list for the angular velocity. [rad/s]
w_deg = np.zeros(len(time))   # list for angle of angular velocity 
tet = np.zeros(len(time)) # A list for the distance.  
x_pos = np.zeros(len(time)) # A list for the xpos
y_pos = np.zeros(len(time)) # A list for the ypos.
x_pos_released = np.zeros(len(time)) # A list for the xpos when released
y_pos_released = np.zeros(len(time)) # A list for the ypos when released
v_x = np.zeros(len(time)) # A list for the xpos
v_y = np.zeros(len(time)) # A list for the ypos.
v_t = np.zeros(len(time)) # List for tangent speed
v_ang = np.zeros(len(time))  # list for angle of the tangent speed
air_res = np.zeros(len(time))    # list for total air resistance
air_res_x = np.zeros(len(time)) # list for X vector of air resistance
air_res_y = np.zeros(len(time))     # list for Y vector of air resistance

w[0] = w0  #start w
tet[0] = math.radians(tet0)  
x_pos[0] = 0   # xpos start
y_pos[0] = 1.8  # approximate hight of a person that is throwing
released=False
  
# Run simulation.
for i in range(1, len(time)):  
    if (released==False):
        w[i] = w[i-1] + torq_elb / I *dt  
        tet[i] = tet[i-1] + w[i] * dt  
        x_pos[i] = x_pos[i-1] + r * math.cos(w[i] - w[i-1]) * dt
        y_pos[i] = y_pos[i-1] + r * math.sin(w[i] - w[i-1]) * dt
    if (tet[i]  > re_angl_rad) and (released==False):
        v_release = w[i] * r
        tet_release = tet[i]
        v_x[i-1] = v_release * math.cos(tet_release)
        v_y[i-1] = v_release * math.sin(tet_release)
        x_pos_released[i-1]=x_pos[i]
        y_pos_released[i-1]=y_pos[i]
        released=True
        w[i] = 0       
    if (released==True):
        v_x[i]=v_x[i-1]
        v_y[i]=v_y[i-1]
        v_t[i] = math.sqrt(v_x[i-1]**2+v_y[i-1]**2)
        v_ang[i] = math.atan(v_y[i-1]/v_x[i-1])
        air_res[i] = 0.5*rho_air*v_t[i]**2*c_drag*area_b ##Force airresitance    
        air_res_x[i] = (air_res[i] * math.cos(v_ang[i]))/m_ball
        air_res_y[i] = (air_res[i] * math.sin(v_ang[i]))/m_ball
        v_x[i]= v_x[i-1] - air_res_x[i] * dt
        v_y[i]= v_y[i-1] - (g + air_res_y[i]) * dt 
        x_pos_released[i] = x_pos_released[i-1] + ((v_x[i-1] + v_x[i])/2) * dt
        y_pos_released[i] = y_pos_released[i-1] + ((v_y[i-1] + v_y[i])/2) * dt
          
# Display results.  
plt.plot(time, tet, color='pink',label='Angular Displacement') 
plt.plot(time, w, color='yellow', label='Angular Velocity')  
plt.plot(time, y_pos_released,color='blue',label='Y_position')
plt.plot(time, x_pos_released,color='purple',label='X_position')
plt.plot(time, v_t, color='black',label='Tangent Velocity')
plt.plot(time, air_res, color='cyan',label='Force Air_res')  
plt.xlabel('Time [s]')  
plt.ylabel('YYYY') 
plt.legend()
plt.xlim(0,1)
plt.ylim(0,10) 
plt.show()  

plt.plot(x_pos_released,y_pos_released,':',color='green' ,label='Trajectory ball')
plt.xlabel('Distance[m]')  
plt.ylabel('Height[m]') 
plt.xlim(0,4)
plt.ylim(0,4)
plt.legend() 
plt.show()

I would probably store the values of each simulation in a dictionary, that way you can write that to file and then just call on a specific release angle and torque to display, as opposed to trying to plot each one. But you just need to work out the logic for the nested loop here.

You could then loop through the results to get the graphs or leave it is an input so you could call a specific combination to output. You'll need to pip install choice to implement the input part of my solution.

Code:

import numpy as np
import matplotlib.pyplot as plt
import math

# Define settings.  
endTime = 60 # The time (seconds) that we simulate.  
dt = 0.001 # The time step (seconds) that we use in the discretization.  
w0 = 0 # The initial velocity [m/s].  
tet0 = 15 # The initial position [m].  
#torq_elb = 70     #[Nm]
torq_elb_list = np.array([50,55,60,65,70,75,80,85,90,95,100])       #[Nm]
I = 0.16        #[kgm^2]  inertia
x_pos = 0       #initial x position
y_pos = 0       # initial y Position
r = 0.3         # lower arm length [m]
#re_angl_deg = 50  # release angle degrees
re_angl_deg_list = np.array([30,32,34,36,38,40,42,44,45,47,49,50,52,54])
g = 9.81 # The gravitational acceleration [m/s^2].    
m_ball = 0.6 # The mass of the ball [kg]. 
rho_air = 1.2041 #[kg/m3]
c_drag =  0.3       # drag coefficient
area_b = (4*math.pi*r**2)/2 # surface area of the ball


results = {}
count = 1
tot = len(torq_elb_list) * len(re_angl_deg_list)
for torq_elb in torq_elb_list:
    for re_angl_deg in re_angl_deg_list:
        results[(torq_elb, re_angl_deg)] = {}
 
        re_angl_rad = math.radians(re_angl_deg ) 
 
        # Set up variables.  
        time = np.arange(0, endTime + dt, dt) # A list with all times we want to plot at.  
        w = np.zeros(len(time)) # A list for the angular velocity. [rad/s]
        w_deg = np.zeros(len(time))   # list for angle of angular velocity 
        tet = np.zeros(len(time)) # A list for the distance.  
        x_pos = np.zeros(len(time)) # A list for the xpos
        y_pos = np.zeros(len(time)) # A list for the ypos.
        x_pos_released = np.zeros(len(time)) # A list for the xpos when released
        y_pos_released = np.zeros(len(time)) # A list for the ypos when released
        v_x = np.zeros(len(time)) # A list for the xpos
        v_y = np.zeros(len(time)) # A list for the ypos.
        v_t = np.zeros(len(time)) # List for tangent speed
        v_ang = np.zeros(len(time))  # list for angle of the tangent speed
        air_res = np.zeros(len(time))    # list for total air resistance
        air_res_x = np.zeros(len(time)) # list for X vector of air resistance
        air_res_y = np.zeros(len(time))     # list for Y vector of air resistance
        
        w[0] = w0  #start w
        tet[0] = math.radians(tet0)  
        x_pos[0] = 0   # xpos start
        y_pos[0] = 1.8  # approximate hight of a person that is throwing
        released=False

  
        # Run simulation.
        for i in range(1, len(time)):  
            if (released==False):
                w[i] = w[i-1] + torq_elb / I *dt  
                tet[i] = tet[i-1] + w[i] * dt  
                x_pos[i] = x_pos[i-1] + r * math.cos(w[i] - w[i-1]) * dt
                y_pos[i] = y_pos[i-1] + r * math.sin(w[i] - w[i-1]) * dt
            if (tet[i]  > re_angl_rad) and (released==False):
                v_release = w[i] * r
                tet_release = tet[i]
                v_x[i-1] = v_release * math.cos(tet_release)
                v_y[i-1] = v_release * math.sin(tet_release)
                x_pos_released[i-1]=x_pos[i]
                y_pos_released[i-1]=y_pos[i]
                released=True
                w[i] = 0       
            if (released==True):
                v_x[i]=v_x[i-1]
                v_y[i]=v_y[i-1]
                v_t[i] = math.sqrt(v_x[i-1]**2+v_y[i-1]**2)
                v_ang[i] = math.atan(v_y[i-1]/v_x[i-1])
                air_res[i] = 0.5*rho_air*v_t[i]**2*c_drag*area_b ##Force airresitance    
                air_res_x[i] = (air_res[i] * math.cos(v_ang[i]))/m_ball
                air_res_y[i] = (air_res[i] * math.sin(v_ang[i]))/m_ball
                v_x[i]= v_x[i-1] - air_res_x[i] * dt
                v_y[i]= v_y[i-1] - (g + air_res_y[i]) * dt 
                x_pos_released[i] = x_pos_released[i-1] + ((v_x[i-1] + v_x[i])/2) * dt
                y_pos_released[i] = y_pos_released[i-1] + ((v_y[i-1] + v_y[i])/2) * dt
                
        # Store simulation results
        results[(torq_elb, re_angl_deg)]['time'] = time
        results[(torq_elb, re_angl_deg)]['tet'] = tet
        results[(torq_elb, re_angl_deg)]['w'] = w
        results[(torq_elb, re_angl_deg)]['y_pos_released'] = y_pos_released
        results[(torq_elb, re_angl_deg)]['x_pos_released'] = x_pos_released
        results[(torq_elb, re_angl_deg)]['v_t'] = v_t
        results[(torq_elb, re_angl_deg)]['air_res'] = air_res
        
        print(f'Finished simulation: torque: {torq_elb}    release: {re_angl_deg}\t{count} of {tot}')
        count+=1
          
# Display results for a particular (torq, re_angl_deg).

#pip install choice
import choice

print('Choose torque.')
torque_choice = int(choice.Menu([str(x) for x in torq_elb_list.tolist()]).ask())

print('Choose release degree.')
release_choice = int(choice.Menu([str(x) for x in re_angl_deg_list.tolist()]).ask())


time = results[(torque_choice, release_choice)]['time']
tet = results[(torque_choice, release_choice)]['tet']
w = results[(torque_choice, release_choice)]['w']
y_pos_released = results[(torque_choice, release_choice)]['y_pos_released']
x_pos_released = results[(torque_choice, release_choice)]['x_pos_released']
v_t = results[(torque_choice, release_choice)]['v_t']
air_res = results[(torque_choice, release_choice)]['air_res']

  
plt.plot(time, tet, color='pink',label='Angular Displacement') 
plt.plot(time, w, color='yellow', label='Angular Velocity')  
plt.plot(time, y_pos_released,color='blue',label='Y_position')
plt.plot(time, x_pos_released,color='purple',label='X_position')
plt.plot(time, v_t, color='black',label='Tangent Velocity')
plt.plot(time, air_res, color='cyan',label='Force Air_res')  
plt.xlabel('Time [s]')  
plt.ylabel('YYYY') 
plt.title('Torque: %s - Release: %s°' %(torque_choice,release_choice))
plt.legend()
plt.xlim(0,1)
plt.ylim(0,10) 
plt.show()  

plt.plot(x_pos_released,y_pos_released,':',color='green' ,label='Trajectory ball')
plt.xlabel('Distance[m]')  
plt.ylabel('Height[m]') 
plt.title('Torque: %s - Release: %s°' %(torque_choice,release_choice))
plt.xlim(0,4)
plt.ylim(0,4)
plt.legend() 
plt.show()

Input:

Choose torque.
Make a choice:
 0: 50
 1: 55
 2: 60
 3: 65
 4: 70
 5: 75
 6: 80
 7: 85
 8: 90
 9: 95

Enter number or name; return for next page

? 0

Choose release degree.
Make a choice:
 0: 30
 1: 32
 2: 34
 3: 36
 4: 38
 5: 40
 6: 42
 7: 44
 8: 45
 9: 47

Enter number or name; return for next page

? 4

Output: