__author__ = 'Lenovo' import numpy as np a = np.pi / 4 g = 9.8 v0 = 30 l1 = np.sin(2*a) * v0 ** 2 / g print 'flight length in Galiley model:' + str(l1) ro = 7600 r = 0.2 u = v0 * np.cos(a) w = 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 x = 0 y = 0 f = False xprev = 0 yprev = 0 def deriv_u_t(u, w): return (-b) * u * np.sqrt(u**2 + w**2) / m def deriv_w_t(u, w): return -g - b * w * np.sqrt(u**2 + w**2) / m while (True): k1u = h * deriv_u_t(u, w) k1w = h * deriv_w_t(u, w) k2u = h * deriv_u_t(u + h * k1u/2, w + h * k1w/2) k2w = h * deriv_w_t(u + h * k1u/2, w + h * k1w/2) k3u = h * deriv_u_t(u + h * k2u/2, w + h * k2w/2) k3w = h * deriv_w_t(u + h * k2u/2, w + h * k2w/2) k4u = h * deriv_u_t(u + h * k3u, w + h * k3w) k4w = h * deriv_w_t(u + h * k3u, w + h * k3w) u += (k1u + 2 * k2u + 2 * k3u + k4u)/6 w += (k1w + 2 * k2w + 2 * k3w + k4w)/6 v = np.sqrt(u**2 + w**2) xprev = x yprev = y x += u * h y += w * h if (y == 0): f = True break if (y < 0): break print 'flight length in Newton model:' if (f): print x else: print (xprev * y - x * yprev) / (y - yprev)