@@ -0,0 +1,1996 @@
+import numpy%0Aimport matplotlib.pyplot as plt%0A%0A# use the sign convention and formulae provided by Lanchester (1909)%0A%0A# calculate the radius of curvature of the flight path at any point%0Adef radius_of_curvature(y, yt, C):%0A#%09return 2./(C*(1./y)**1.5 - 2./3./yt)%0A        return yt / (1./3 - C/2.*(yt/y)**1.5)%0A%0A# given the coordinates of a point, the center around which it much revolve,%0A# and the angle of rotation, calculate the new position of the point%0Adef rotate(x, y, xCenter, yCenter, angle):%0A%09dx = x - xCenter%0A%09dy = y - yCenter%0A%09xNew = dx*numpy.cos(angle) - dy*numpy.sin(angle)%0A%09yNew = dx*numpy.sin(angle) + dy*numpy.cos(angle)%0A%09return xCenter + xNew, yCenter + yNew%0A%0A# plot the flight path. The inputs are the trim height and the initial position %0A# and orientation of the aircraft%0Adef plot_flight_path(yt, y0, theta0, case_name='dummy'):%09%0A%09# arrays to store the coordinates of the flight path%0A%09N = 5000%0A%09y = numpy.zeros(N)%0A%09x = numpy.zeros(N)%0A%0A%09# set initial conditions%0A%09y%5B0%5D = y0%0A%09x%5B0%5D = 0.%0A%09theta = theta0%0A%0A%09# calculate the constant C%0A%09#C = numpy.sqrt(y%5B0%5D)*(numpy.cos(theta) - y%5B0%5D/yt/3.)%0A        C = (numpy.cos(theta) - 1./3*y%5B0%5D/yt)*(y%5B0%5D/yt)**.5%0A%0A%09# incremental distance along the flight path%0A        ds = 1%0A%0A%09# obtain the curve coordinates%0A%09for i in xrange(1,N):%0A%09%09normal = numpy.array(%5Bnumpy.cos(theta+numpy.pi/2.), numpy.sin(theta+numpy.pi/2.)%5D)%0A%09%09R = radius_of_curvature(y%5Bi-1%5D, yt, C)%0A%09%09center = numpy.array(%5Bx%5Bi-1%5D+normal%5B0%5D*R, y%5Bi-1%5D+normal%5B1%5D*R%5D)%0A%09%09dtheta = ds/R%0A%09%09x%5Bi%5D, y%5Bi%5D = rotate(x%5Bi-1%5D, y%5Bi-1%5D, center%5B0%5D, center%5B1%5D, dtheta)%0A%09%09theta = theta + dtheta%0A%0A%09# generate a plot%0A%09plt.figure(figsize=(11,8))%0A        plt.plot(x, -y,'k-',label=%22$z_t=%5C %25.1f,%5C%5C,z_1=%5C %25.1f,%5C%5C,%5C%5Ctheta_1=%5C %25.2f$%22 %25 (yt, y%5B0%5D, theta0))%0A%09plt.title(%22Flight path for C = %25.4f%22 %25 C)%0A        plt.xlabel(%22$x$%22)%0A%09plt.ylabel(%22$z$%22)%0A#%09plt.axis(%5B-100, 1000, -300, 50%5D)%0A%09plt.legend()%0A        plt.show()%0A%0A%0A#plot_flight_path(600., 600., -numpy.pi, %221_looped_path%22)%0A#plot_flight_path(64., 16., 0., %222_inflected_path%22)%0A%0A

add phugoid plotter