Speed-to-fly analysis for different glider classes

polar.py

# Units:
# Distance: metres
# Time: seconds
# All vertical velocities are positive upwards, i.e. climb rates are always
# postive and sink rates are always negative
# To convert km/h into m/s multiply by 3.6

# Notation:
# vx: Horizontal velocity along the course line
# vxmin: Minimum velocity (velocity at minimum sink)
# vxmax: Maximum horizontal velocity (maximum speed)
# vz: Vertical velocity (always negative, as the glider is always sinking)
# cvz: Climb velocity (thermal strength minus glider sink rate)


import math
import matplotlib.pyplot as plt
import numpy
import scipy.optimize


def kph(*args):
    """Convert kilometres per hour to metres per second"""
    return map(lambda x: x / 3.6, args)


class Paraglider(object):

    def __init__(self, label, vxs, vzs):
        self.label = label
        self.vxmin = min(vxs)
        self.vxmax = max(vxs)
        self.vzmin = min(vzs)
        self.polar = numpy.polyfit(vxs, vzs, 2)

    def real_speed(self, vx):
        return max(min(vx, self.vxmax), self.vxmin)

    def sink_rate(self, vx):
        return numpy.polyval(self.polar, vx)

    def course_speed(self, vx, cvz):
        vxreal = self.real_speed(vx)
        vzreal = self.sink_rate(vxreal)
        return vxreal * cvz / (cvz - vzreal)

    def inverse_course_speed(self, vx, cvz):
        if vx < self.vxmin:
            return 0
        if self.vxmax < vx:
            return 0
        return -self.course_speed(vx, cvz)

    def speed_to_fly(self, cvz):
        return scipy.optimize.fmin_powell(lambda vx: self.inverse_course_speed(vx, cvz), (self.vxmin + self.vxmax) / 2, disp=False)

    def inverse_glide(self, vx):
        vz = self.sink_rate(vx)
        return vx / vz

    def best_glide(self):
        vx = scipy.optimize.fmin_powell(self.inverse_glide, (self.vxmin + self.vxmax) / 2, disp=False)
        vz = self.sink_rate(vx)
        return vx / vz

r11 = Paraglider('Mantra R11', kph(30, 38, 50, 70), [-0.9, -1.1, -1.3, -2.8])
enzo = Paraglider('EnZo', kph(30, 38, 50, 60), [-0.9, -1.1, -1.3, -2.0])
m4 = Paraglider('Mantra M4', kph(30, 36, 50, 58), [-0.9, -1.3, -1.4, -2.2])
paragliders = [m4, enzo, r11]
dcvz = 1.5

vxmin = 20
vxmax = 80

plt.figure(1)
for paraglider in paragliders:
    vxs = numpy.arange(30 / 3.6, 70 / 3.6, 0.1 / 3.6)
    vxs = vxs[(paraglider.vxmin <= vxs) & (vxs <= paraglider.vxmax)]
    vzs = numpy.array(map(paraglider.sink_rate, vxs))
    plt.plot(3.6 * vxs, vzs, label=paraglider.label)
plt.axis([vxmin, vxmax, -3.0, 0.0])
plt.grid(True)
plt.legend()
plt.title('Polar curves')
plt.xlabel('Speed (km/h)')
plt.ylabel('Sink rate (m/s)')
plt.savefig('1-polar.png')

plt.figure(2)
cvzs = numpy.arange(0.1, 5.1, 0.1)
for paraglider in paragliders:
    speeds = numpy.array(map(paraglider.speed_to_fly, cvzs))
    plt.plot(cvzs, 3.6 * speeds, label=paraglider.label)
plt.axis([0, 5, 40, 75])
plt.grid(True)
plt.legend(loc='lower right')
plt.title('Speed to fly')
plt.xlabel('Thermal strength (m/s)')
plt.ylabel('Speed to fly (km/h)')
plt.savefig('2-speed-to-fly.png')

plt.figure(3)
cvzs = numpy.arange(0.1, 5.1, 0.1)
for paraglider in paragliders:
    speeds = numpy.array(map(lambda cvz: paraglider.course_speed(paraglider.speed_to_fly(cvz), cvz), cvzs))
    plt.plot(cvzs, 3.6 * speeds, label=paraglider.label)
plt.axis([0, 5, 10, 45])
plt.grid(True)
plt.legend(loc='lower right')
plt.title('Speed around course')
plt.xlabel('Thermal strength (m/s)')
plt.ylabel('Course speed (km/h)')
plt.savefig('3-course-speed.png')

plt.figure(4)
cvzs = numpy.arange(0.1, 5.1, 0.1)
for paraglider in paragliders:
    speeds1 = numpy.array(map(lambda cvz: paraglider.course_speed(paraglider.speed_to_fly(cvz), cvz), cvzs))
    speeds2 = numpy.array(map(lambda cvz: paraglider.course_speed(paraglider.speed_to_fly(cvz * dcvz), cvz * dcvz), cvzs))
    plt.plot(cvzs, 3.6 * (speeds2 - speeds1), label=paraglider.label)
plt.axis([0, 5, 4, 7])
plt.grid(True)
plt.legend(loc='lower right')
plt.title('Absolute speed difference with 50% better climb rates')
plt.xlabel('Thermal strength (m/s)')
plt.ylabel('Absolute speed difference (km/h)')
plt.savefig('4-absolute-speed-difference.png')

plt.figure(5)
cvzs = numpy.arange(0.1, 5.1, 0.1)
for paraglider in paragliders:
    speeds1 = numpy.array(map(lambda cvz: paraglider.course_speed(paraglider.speed_to_fly(cvz), cvz), cvzs))
    speeds2 = numpy.array(map(lambda cvz: paraglider.course_speed(paraglider.speed_to_fly(cvz * dcvz), cvz * dcvz), cvzs))
    plt.plot(cvzs, 60 / (3.6 * (speeds2 - speeds1)), label=paraglider.label)
plt.axis([0, 5, 5, 15])
plt.grid(True)
plt.legend(loc='lower left')
plt.title('Time to catch up 1km if with 50% better climb rates')
plt.xlabel('Thermal strength on direct line (m/s)')
plt.ylabel('Time to catch up 1km with 50% better climb rates (min)')
plt.savefig('5-catch-up-time.png')

def ellipse_area(r):
    return math.pi * r * numpy.sqrt((r ** 2 - 1) / 4) / 4

plt.figure(6)
cvzs = numpy.arange(0.1, 5.1, 0.1)
for paraglider in paragliders:
    speeds1 = numpy.array(map(lambda cvz: paraglider.course_speed(paraglider.speed_to_fly(cvz), cvz), cvzs))
    speeds2 = numpy.array(map(lambda cvz: paraglider.course_speed(paraglider.speed_to_fly(cvz * dcvz), cvz * dcvz), cvzs))
    plt.plot(cvzs, ellipse_area(speeds2 / speeds1), label=paraglider.label)
plt.axis([0, 5, 0.2, 0.4])
plt.grid(True)
plt.legend(loc='upper right')
plt.title('Option area with 50% better climb rates')
plt.xlabel('Thermal strength on direct line (m/s)')
plt.ylabel('Option area with 50% better climb rates')
plt.savefig('6-option-area.png')