-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathMain.py
More file actions
129 lines (106 loc) · 3.61 KB
/
Copy pathMain.py
File metadata and controls
129 lines (106 loc) · 3.61 KB
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
from math import sin, cos, radians
import matplotlib.pyplot as plt
name = input('''Choose an option:
m - solveMid-Point
e - solveEuler
''')
theta = radians(float(input("Enter Initial Launch Angle of the Object from 0 Degree to 360 Degree: ")))
DT = 0.01 # The delta time step of the motion
v = float(input("Enter Launch velocity in (m/s): ")) # Launch velocity
G = -9.81 # Gravitational Acceleration
k = float(input("Enter Resistance Force Consant: ")) # Resistance force constant
m = float(input("Enter Mass of the object in(kg): ")) # Mass of the object
t = 0.0
x = 0.0
y = 0.0
vx,vy = v*cos(theta), v*sin(theta) # Initial vertical and horizontal speeds
a_x1 = -(k*vx)/m # Horizontal forces with drag force case.
a_y1 = (G-(k*vy)/m) # Both drag force and gravitational force.
Max_Height = 0.0 # Maximum Height,initially set to zero
xa_1 = list()
ya_1 = list()
xa_2 = list()
ya_2 = list()
xa_3 = list()
ya_3 = list()
def euler(y, x, t, vx, vy, xa_1, ya_1, Max_Height):
while y >= 0:
xa_1.append(x)
ya_1.append(y)
y = y + vy * DT
x = x + vx * DT
vy = vy + a_y1 * DT
vx = vx + a_x1 * DT
t = t + DT
if y > Max_Height:
Max_Height = y
xa_1.append(x)
ya_1.append(y)
return Max_Height
def analitical(y, x, t, vx, vy, xa_2, ya_2, Max_Height):
while y >= 0:
xa_2.append(x)
ya_2.append(y)
vx_2 = vx + G*0.5*t*t*DT
vy_2 = vy + G*0.5*t*t*DT
y += vy_2 * DT
x += vx_2 * DT
vy = vy + a_y1 * DT
vx = vx + a_x1 * DT
t += DT
if y > Max_Height:
Max_Height = y
xa_2.append(x)
ya_2.append(y)
return Max_Height
def mid_point(y, x, t, vx, vy, xa_3, ya_3, Max_Height):
while y >= 0:
xa_3.append(x)
ya_3.append(y)
vx_2 = vx + (a_x1*(DT/2))
vy_2 = vy + (a_y1*(DT/2))
y = y + vy_2 * DT
x = x + vx_2 * DT
vy = vy + a_y1 * DT
vx = vx + a_x1 * DT
t = t + DT
if y > Max_Height:
Max_Height = y
xa_3.append(x)
ya_3.append(y)
return Max_Height
if name == 'm':
h1 = analitical(y, x, t, vx, vy, xa_1, ya_1, Max_Height)
h0 = mid_point(y, x, t, vx, vy, xa_3, ya_3, Max_Height)
if h1 > h0:
Max_Height = h1
else:
Max_Height = h0
elif name == 'e':
h1 = analitical(y, x, t, vx, vy, xa_1, ya_1, Max_Height)
h0 = euler(y, x, t, vx, vy, xa_2, ya_2, Max_Height)
if h1 > h0:
Max_Height = h1
else:
Max_Height = h0
elif name == 'all':
h1 = analitical(y, x, t, vx, vy, xa_1, ya_1, Max_Height)
h2 = mid_point(y, x, t, vx, vy, xa_3, ya_3, Max_Height)
h0 = euler(y, x, t, vx, vy, xa_2, ya_2, Max_Height)
if h0 > h1 and h0 > h2:
Max_Height = h0
elif h0 <= h1 and h0 >= h2:
Max_Height = h1
else:
Max_Height = h2
else:
print('Bad argument')
plt.figure()
plt.plot(xa_1, ya_1, color = 'r') #Analitic
plt.plot(xa_2, ya_2, color = 'b') #Euler
plt.plot(xa_3, ya_3, color = 'g') #Mid-Point
plt.ylim(0, Max_Height + Max_Height/10)
plt.ylabel('Horizontal Axis (m)')
plt.xlabel('Vertical Axis (m)')
plt.title('Trajectory Motion of Object')
plt.show()