__author__ = 'vladimirtsvetkov' import numpy as np a = np.pi / 4 g = 10 v0 = 40 ro = 7600 r = 0.2 #first model(simple) x = np.tan(a) / (g / (2 * v0 * v0 * np.cos(a) * np.cos(a))) print 'simple model:' print x #second model(complex) vx = v0 * np.cos(a) vy = v0 * np.sin(a) h = 0.01 b = 0.16 * 1.15 * np.pi * r * r / 2.0 m = ro * 4 * np.pi * r * r * r /3 def fvx(vx, vy): return ((-b * vx * np.sqrt(vx*vx + vy*vy))/m) def fvy(vx, vy): return -g - ((b * vy * np.sqrt(vx*vx + vy*vy))/m) y = 0 x = 0 flag = False xprev = 0 yprev = 0 xmax = 120 while (True): k11 = h * fvx(vx, vy) k12 = 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 = x yprev = y x += vx * h y += vy * h if (x == xmax): flag = True break if (x > xmax): break if (flag): print 'complex model1:' print x else: print 'complex model2:' print x #print (xprev * y - x * yprev)/(y - yprev)