ProjectMagnus/TestMark.py at master · MarkG98/ProjectMagnus · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
from modsim import *
from matplotlib import pyplot as plt

# declare all units
m = UNITS.meter
s = UNITS.second
kg = UNITS.kilogram
degree = UNITS.degree
radian = UNITS.radian

hero_height = 1.5;

condition = Condition(x = 2 * m,
                      y = (365.76 + hero_height) * m,
                      g = 9.8 * m/s**2,
                      mass = 145e-3 * kg,
                      diameter = 73e-3 * m,
                      rho = 1.2 * kg/m**3,
                      C_d = 0.3,
                      angle = -10 * degree,
                      velocity = 45 * m/s,
                      w = 33.4425156891 * radian/s,
                      duration = 15 * s)

def make_system(condition):
    theta = np.deg2rad(condition.angle)
    vx, vy = pol2cart(theta, condition.velocity)
    init = State(x = condition.x, y = condition.y, vx = vx, vy = vy)
    area = np.pi * (condition.diameter / 2) ** 2
    ts = linspace(0, condition.duration, 101)

    return System(init = init, g = condition.g, mass = condition.mass, area = area, rho = condition.rho, C_d = condition.C_d, ts = ts)

baseballSystem = make_system(condition)

def slope_func(state, t, system):
    x, y, vx, vy = state
    unpack(system)

    a_grav = Vector(0, -g, 0) # acceleration due to gravity
    v_ball = Vector(vx, vy,0) # ball velocity vector

    # creates unit vector in direction of the magus force
    x, y = pol2cart(v_ball.angle - (pi / 2) * radian, 1)
    magnus_direction = Vector(x, y, 0)

    w_vector = Vector(0, 0, - condition.w)

    # calculates acceleration due to the magnus force
    #f_magnus = (pi ** 2) * ((condition.diameter / 2) ** 3) * rho * v_ball.mag * condition.w * magnus_direction.hat()
    f_magnus = ((pi ** 2) * ((condition.diameter / 2) ** 3) * rho) * w_vector.cross(v_ball)
    a_magnus = f_magnus / mass

    # calculates acceleration due to drag
    f_drag = -rho * v_ball.mag * v_ball * C_d * (area / 2)
    a_drag = f_drag / mass

    a_total = a_grav + a_magnus + a_drag

    return v_ball.x, v_ball.y, a_total.x, a_total.y

run_odeint(baseballSystem, slope_func)
xs = baseballSystem.results.x
ys = baseballSystem.results.y

def sweep_func():
  for angle in linspace(-10,-80,2):
    condition.set(angle=angle * degree)
    for velocity in linspace(30,150,2):
      condition.set(velocity=velocity* m/s)
      for w in linspace(20,120,2):
        condition.set(w=w* radian/s)
        baseballSystem = make_system(condition)
        run_odeint(baseballSystem,slope_func)
        xs = baseballSystem.results.x
        ys = baseballSystem.results.y
        plt.xlabel('x (m)')
        plt.ylabel('y (m)')
        plt.title('Ball flight Path')
        plt.xlim([0,70])
        plt.ylim([0,400])
        plt.grid(True)
        plt.axvline(x = 2, color='r')
        plt.plot(xs, ys)
        plt.show()

def interpolate_range(results,y_value):
  ys = results
  descent = ys
  T = interp_inverse(descent)
  t_land = T(y_value)
  return t_land

def error_func(w):
  condition.set(w = w * radian/s)
  baseballSystem = make_system(condition)
  run_odeint(baseballSystem, slope_func)
  X = interpolate(baseballSystem.results.x)
  t_land = interpolate_range(baseballSystem.results.y, 0)
  return X(t_land) + 0.78

def heightAtDoor(system):
    xs = system.results.x
    ys = system.results.y

    t_maxRange = xs.idxmax()

    finalDescent = xs.loc[t_maxRange:]
    T = interp_inverse(finalDescent)
    t_throughDoor = T(2)

    Y = interpolate(ys, kind='cubic')
    return Y(t_throughDoor)


#value = error_func(w = 67, interpolate_range = interpolate_range)
solution = fsolve(error_func, 67)
w_ideal = solution[0]
print(w_ideal)

angle_array = linspace(0, 60, 13)

def angleSweep(system, angle_array):
    for angle in angle_array:
        condition.set(angle = angle * degree)
        solution = fsolve(error_func, 60)
        spinForButt = solution[0]

        condition.set(w = spinForButt * radian/s)
        system = make_system(condition)
        run_odeint(system, slope_func)

        height = heightAtDoor(system)

        print("angle = ", angle, "w = ", condition.w, "height at door = ", height, '\n')

angleSweep(baseballSystem, angle_array)


plt.xlabel('x (m)')
plt.ylabel('y (m)')
plt.title('Ball flight Path')
plt.xlim([0,70])
plt.ylim([0,400])
plt.grid(True)
plt.axvline(x = 2, color='r')
plt.plot(xs, ys)
plt.show()