import math
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

g = 9.81
μ = 0.12
s = 0.2
R = 0.1085
m = 7
v0 = 8
Θ = 0
ωx = -7.76
ωy = -28.98
ωz = 0
Ix = 0.033
Iy = 0.033
Iz = 0.033

import math

def func(X, t, M, g, R, m, Ix, Iy, Iz) :
    if math.sqrt((M*g/math.sqrt((X[2]-R*X[5])**2+(X[3]+R*X[4])**2)*math.sqrt
                  ((X[2]-R*X[5])**2))**2+M*g/math.sqrt ((X[2]-R*X[5])**2+(X[3]+R*X[4])**2)
                 *math.sqrt((X[3]+R*X[4])**2)**2) < s*m*g:
        dXdt = [X[2], X[3],
                -M*g/math.sqrt((X[2]-R*X[5])**2+(X[3]+R*X[4])**2)*math.sqrt((X[2]-R*X[5])**2),
                -M*g/math.sqrt((X[2]-R*X[5])**2+(X[3]+R*X[4])**2)*math.sqrt((X[3]+R*X[4])**2),
                (-R*M*m*g/math.sqrt((X[2]-R*X[5])**2+(X[3]+R*X[4])**2)
                 *math.sqrt((X[3]+R*X[4])**2)+(Iy - Iz)*X[5]*X[6])*X[7]/Ix,
                (R*M*m*g/math.sqrt((X[2]-R*X[5])**2+(X[3]+R*X[4])**2)
                 *math.sqrt((X[2]-R*X[5])**2)+(Iz - Ix)*X[6]*X[4])*X[8]/Iy,
                (Ix-Iy)*X[4]*X[5]*X[9]/Iz,
                X[6]*X[8]-X[5]*X[9],
                X[4]*X[9]-X[6]*X[7],
                X[5]*X[7]-X[4]*X[8]]
    else:
        dXdt =[X[2], X[3],
               R*(R*μ*m*g/math.sqrt((X[2]-R*X[5])**2+(X[3]+R*X[4])**2)
                  *math.sqrt ((X[2]-R*X[5])**2)+(Iz-Ix)*X[6]*X[4])*X[8]/Iy,
               -R*(-R*M*n*g/math.sqrt((X[2]-R*X[5])**2+(X[3]+R*X[4])**2)
                   *math.sqrt((X[3]+R*X[4])**2)+(Iy-Iz)*X[5]*X[6])*X[7]/Ix,
               (-R*M*m*g/math.sqrt((X[2]-R*X[5])**2+(X[3]+R*X[4])**2)
                *math.sqrt ((X[3]+R*X[4])**2)+(Iy-Iz)*X[5]*X[6])*X[7]/Ix,
               (R*M*n*g/math.sqrt((X[2]-R*X[5])**2+(X[3]+R*X[4])**2)
                *math.sqrt ((X[2]-R*X[5])**2)+(Iz-Ix)*X[6]*X[4])*X[8]/Iy,
               (Ix-Iy) *X[4]*X[5]*X[9]/Iz,
               X[6]*X[8]-X[5]*X[9],X[4]*X[9]-X[6]*X[7],X[5]*X[7]-X[4]*X[8]]
    return dXdt

X0 = [0,0,v0*math.cos(math.radians(0)),v0*math.sin(math.radians(0)), ωx, ωy, ωz,1,1,1]
t = np.arange(0, 10, 0.01)

X = odeint(func, X0, t, args=(μ, g, R, m, Ix, Iy, Iz))

def plot_bowling_pins():
    pin_coordinates = np.array([
        [18.29, 0],
        [18.554, 0.1524],
        [18.554, -0.1524],
        [18.818, 0.30481],
        [18.818, -0.3048],
        [18.818, 0],
        [19.082, 0.1524],
        [19.082, -0.1524],
        [19.082, 0.4572],
        [19.082, -0.45721],
    ])

    plt.figure(figsize=(20.29, 1.066))
    plt.scatter(pin_coordinates[:, 0], pin_coordinates[:, 1], marker='o', color=
               'black', s=30)

plot_bowling_pins()

plt.plot(X[:,0], X[:,1])
plt.xlim(0, 20.29)
plt.ylim (-0.533, 0.533)
plt.xlabel('x', fontsize=10)
plt.ylabel('y', fontsize=10)

plt.title('Bowling Ball Motion')
plt.show()