import numpy as np v0 ro np pi 81 40 7600 print Galiley str np tan v0

 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
import numpy as np
a, g, v0, ro, r = np.pi / 4, 9.81, 40, 7600, 0.2
print 'Galiley: \n%s' % str(np.tan(a) / (g / (2 * v0 * v0 * np.cos(a) * np.cos(a))))
#second model(Newton)
vx, vy, h = v0 * np.cos(a), v0 * np.sin(a), 0.01
b = 0.15 * 1.15 * np.pi * r **2 / 2.0
m = ro * 4 * np.pi * r **3 / 3
fvx = lambda vx, vy: ((-b * vx * np.sqrt(vx*vx + vy*vy))/m)
fvy = lambda vx, vy: -g - ((b * vy * np.sqrt(vx*vx + vy*vy))/m)
y, x = 0, 0
while (True):
k11, k12 = h * fvx(vx, vy), h * fvy(vx, vy)
k21 = h * fvx(vx + h * k11/2, vy + h * k12/2)
k22 = h * fvy(vx + h * k11/2, vy + h * k12/2)
k31 = h * fvx(vx + h * k21/2, vy + h * k22/2)
k32 = h * fvy(vx + h * k21/2, vy + h * k22/2)
k41 = h * fvx(vx + h * k31, vy + h * k32)
k42 = h * fvy(vx + h * k31, vy + h * k32)
vx += (k11 + 2 * k21 + 2 * k31 + k41)/6
vy += (k12 + 2 * k22 + 2 * k32 + k42)/6
v = np.sqrt(vx*vx + vy*vy)
xprev, yprev = x, y
x += vx * h
y += vy * h
if (y == 0):
print 'Newton:\n%s' % str(x)
break
if (y < 0):
print 'Newton:\n%s' % str((xprev * y - x * yprev)/(y - yprev))
break