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