コード（3次元でのショット軌道シミュレーション）

def search_shot(shot_info, shot_angle, contact_height, net_height, shot_distance):
    ball_radius = 0.033  # ボールの半径 (メートル)
    A = np.pi * (ball_radius**2)  # ボールの断面積 (m^2)
    m = 0.0577  # ボールの質量 (kg)

    sv, spin, spin_axis = shot_info
    shot_v = sv / 3.6

    # サーブの斜め回転に対応するため、回転を物理モデルに組み込む
    w = (spin * 2 * np.pi) / 60
    w_axis = spin_axis  # 回転軸

    net_length = shot_distance[0]
    net_height += ball_radius

    bound_length = shot_distance[1]
    after_bound_length = shot_distance[2]

   # 物理定数
    g = 9.81  # 重力加速度 (m/s^2)
    rho = 1.205  # 空気密度 (kg/m^3)
    e = 0.75  # 反発係数
    myu = 0.72  # 摩擦係数
    micron = 1.822 * 10 ** (-5)  #粘度

    cos0 = (23.77/2 + 6.4) / bound_length
    sin0 = np.sqrt(1 - cos0 ** 2)

    v0 = np.zeros(6)
    v0[0] = -0.3
    v0[2] = contact_height
    v0[3] = shot_v * math.cos(math.radians(shot_angle)) * sin0
    v0[4] = shot_v * math.cos(math.radians(shot_angle)) * cos0
    v0[5] = shot_v * math.sin(math.radians(shot_angle))

    t_min = 0
    t_max = 3
    dt = 0.001
    t = np.arange(t_min, t_max, dt)
    t_span = (t_min, t_max)

    # ショット到達時間
    shot_time = 0
    reach_time = []
    hitting_position = []

    try:
        solved_data = solve_ivp(f, t_span, v0, t_eval=t, args=(ball_radius, w, w_axis, rho, m, g, A, micron))

        flag = True
        not_net_pass = "not checked"

        for i in range(len(solved_data.t)):
            x_point = solved_data.y[0, i]
            y_point = solved_data.y[1, i]
            z_point = solved_data.y[2, i]
            x_velocity = solved_data.y[3, i]
            y_velocity = solved_data.y[4, i]
            z_velocity = solved_data.y[5, i]
            xy_length = np.sqrt(x_point ** 2 + y_point ** 2)
            xy_velocity = np.sqrt(x_velocity ** 2 + y_velocity ** 2)

            if (not_net_pass != "checked") and (y_point >= 23.77/2):
                if (solved_data.y[2, i-1] <= net_height):
                    return [np.nan] * 4
                else:
                    not_net_pass = "checked"

            if z_point <= ball_radius:
                if xy_length <= net_length:
                    return [np.nan] * 4
                elif y_point > 23.77/2 + 6.4:
                    return [np.nan] * 4
                elif x_point < 0 or x_point >8.23/2:
                    return [np.nan] * 4
                elif (x_point < (8.23/2 - 1)) or (y_point < (23.77/2 + 6.4 - 2)):
                    return [np.nan] * 4
                elif (x_point >= (8.23/2 - 1)) or (y_point >= (23.77/2 + 6.4 - 2)):
                    shot_time = solved_data.t[i]
                    y2 = np.zeros(6)
                    y2[0] = x_point
                    y2[1] = y_point
                    y2[2] = z_point
                    y2[3] = x_velocity + myu * (x_velocity / xy_velocity) * (e + 1) * z_velocity
                    y2[4] = y_velocity + myu * (y_velocity / xy_velocity) * (e + 1) * z_velocity
                    y2[5] = - e * z_velocity
                    if xy_velocity - w > 0:
                        w = w - (3 * myu * ball_radius * (1 + e) * z_velocity) / (2 * ball_radius ** 2)
                    elif xy_velocity - w < 0:
                        w = w + (3 * myu * ball_radius * (1 + e) * z_velocity) / (2 * ball_radius ** 2)
                    w_axis = 90
                    break
                else:
                    return [np.nan] * 4

        # print(f"bound: {sv}, {spin_axis}, {spin}, {shot_angle}")

        solved_data = solve_ivp(f, t_span, y2, t_eval=t, args=(ball_radius, w, w_axis, rho, m, g, A, micron))

        check_yline = 23.77
        for i in range(len(solved_data.t)):
            if len(solved_data.t) == i - 1:
                return [np.nan] * 4
            if solved_data.y[1, i] >= check_yline:
                reach_time.append(shot_time + solved_data.t[i])
                hitting_position.append(solved_data.y[0, i])
                hitting_position.append(solved_data.y[1, i])

                return [shot_time] + reach_time + hitting_position
    except Exception as e:
        tb = sys.exc_info()[2]
        print(e.with_traceback(tb))
        return [np.nan] * 4


def f(t, y, ball_radius, w, w_axis, rho, m, g, A, micron):
    u, udot = y[:3], y[3:]
    v = np.sqrt(udot[0]**2 + udot[1]**2 + udot[2]**2)

    v_xy = np.sqrt(udot[0] ** 2 + udot[1] ** 2)
    v_yz = np.sqrt(udot[1] ** 2 + udot[2] ** 2)

    cosXY = udot[1] / v_xy
    sinXY = udot[0] / v_xy


    # レイノルズ数
    Re = (rho * (2 * ball_radius) * v) / micron

    # 抗力係数
    Cd = 24 / Re
    if Re > 3 * 10 ** 5:
        Cd = 0.1
    elif Re > 10 ** 3:
        Cd = 0.44

    # 抗力
    Fd = (1 / 2) * Cd * A * rho * v * udot

    # 揚力
    Cl = 1 / (2.02 + 0.981 * (v / (ball_radius * w)))
    w_vec = [
        w * math.sin(math.radians(w_axis)) * sinXY,
        w * math.sin(math.radians(w_axis)) * cosXY,
        w * math.cos(math.radians(w_axis))
    ]
    omega = [
        [0, w_vec[2], -w_vec[1]],
        [-w_vec[2], 0, w_vec[0]],
        [w_vec[1], -w_vec[0], 0]
    ]
    Fl = (1 / 2) * Cl * A * rho * (v ** 2) * (np.dot(omega, udot)) / np.linalg.norm(np.dot(omega, udot))

    udotdot_x = (1 / m) * (- Fd[0] + Fl[0])
    udotdot_y = (1 / m) * (- Fd[1] + Fl[1])
    udotdot_z = (1 / m) * (- Fd[2] + Fl[2]) - g

    dydt = np.hstack([udot, udotdot_x, udotdot_y, udotdot_z])
    return dydt

columns = ["shot velocity", "spin rate", "spin angle", "shot angle"]
df = pd.DataFrame(columns=columns)
v = [160, 170]
s = [x for x in range(2000, 4001, 50)]
a = [x for x in range(0, 71, 5)]
b = [round(x, 1) for x in np.arange(-5.5, -3.9, 0.1)]

for i, comb in enumerate(itertools.product(v, s, a, b)):
    df.loc[i, columns] = list(comb)

#セットアップ
tqdm.pandas()

contact_height = 2.5
l = 23.77/2+6.4

result_df = pd.DataFrame()

for shot_width in np.arange(0, 1.6, 0.1):
    shot_width = round(shot_width, 1)
    print(f"shot width: {shot_width}")
    bound_length = np.sqrt(l ** 2 + (8.23/2 - shot_width) ** 2)
    shot_length = bound_length + 5
    net_length = bound_length * (((23.77/2)) / l)
    net_height = min(1.07, 0.914 + (0.156 * (abs(((np.sqrt((net_length ** 2 - ((23.77 / 2)) ** 2))) - (8.23 / 2 - shot_width))) / 5.029)))

    result_columns = ["shot time", "reach time", "x", "y"]
    df["shot width"] = shot_width
    df[result_columns] = df[columns].progress_apply(
        lambda x: search_shot([x[0], x[1], x[2]], x[3], contact_height, net_height, [net_length, bound_length, (shot_length - bound_length)]),
        axis=1,
        result_type="expand"
    )

    result_df = pd.concat([result_df, df.dropna()])