Dummy Gliding Simulator with Python

As we shown in previous post, gliding flight mechanics can be easily modeled and simulated. At that post, we used Scilab for solving dynamic system for a determined range of time and with constant lifting coefficient (CL). This time we are going to improve code using Python, flavoring simulation with PyGame and giving a dummy control to user.

Let's start preparing software platform for building our dummy flight simulator:

• Python is a very powerful and clear programming language. If you have programmed before, you'll have no problem for understanding python codes;


• SciPy is Python library for scientific programming. It is necessary download NumPy and, actual, SciPy;


• PyGame "is a set of Python modules designed for writing games".

All these tools are available for both Linux and Windows users. Once you have installed these required software, the following instructions are valid independently of your operating system.

Designing code

Using Python, which has a different purpose than pure scientific tools as Scilab, we got more flexibility to modularize code and to make it become more reusable for future improvements.

So we are going to design our simulator with object oriented programming (OOP). If you are not familiar with this paradigm, we are going to give a very simple explanation that will hopefully become more clear with the codes written later.

OOP can be understood as a way to crack your code in little modules called classes, which are recipes for keeping together a set of data and a set of functions that change these data. Semantically, when modeling a problem, we bound the entities of the problem in classes, giving to them attributes and methods, i.e. actions that can change the entities behavior.

OK! Let's stop being abstract. See our physics model for gliding. We can identify an entity for Atmosphere stuff. So we can create the class Atmosphere with attributes about the values of temperature and density, for example; and with methods for changing temperature and density according height. Another class could be Aircraft with attributes like position and velocity, and methods to increase or decrease CL (this would be equivalent to give control to make aircraft to go up or to go down).

We call object, an individual of a given class. In other words, class is a recipe to build an object. For example, since we have an Aircraft class, we could create two different objects depending on values of their attributes, different polar drag, for instance.

Our model will contain the following main classes: Atmosphere, Aircraft and Flight.

Atmosphere class

For simplicity we're going to put all our code in just one python module (.py file extension). So, choose your favorite text editor and start putting the following code:

#!/usr/bin/python
# -*- coding: utf-8 -*-

from math import exp

class Atmosphere(object):
    T0 = 288.15
    H0 = 0
    g0 = 9.80665
    L = -6.5e-3
    R = 287.05287
    rho0 = 1.225
    
    def temperature(self,H):
        T = self.T0
        if H<=11000:
            T = self.T0 + self.L*(H-self.H0)
        else:
            T = self.temperature(11000)
        return T

    def density(self,H):
        rho = self.rho0
        if H<=11000:
            rho = self.rho0*(self.temperature(H)/self.T0)\
                  **(-(1+self.g0/(self.R*self.L)))
        else:
            rho = self.density(11000)*exp(-self.g0*\
                  (H-11000)/(self.R*self.temperature(11000)))
        return rho

Two initial lines are just some conventions and they can be even neglected for our purposes. Then we import exponential function (exp) from math library. Finally, we start code Atmosphere class. See that this code is very similar to Scilab code, but this time we have unified all atmospheric stuff in just one code structure, the class. Realize that methods are declared with def keyword and for accessing attribute values we have to use self identifier.

To test the code you can save the module as gliding.py and call it by command line:

$ python -i gliding.py 
>>> a = Atmosphere()
>>> a.density(600)
1.1559768882668726
>>> a.temperature(600)
284.25

In this test you are creating an object a from Atmosphere class and calling its methods density and, after, temperature for calculating these measures at height 600m.

Aircraft class

Now, the entity Aircraft. We should put CL and CD as methods, receiving angle of attack as parameter, but for our gliding flight we're going to consider constant CL or to change it directly. In other words, in real aircrafts, pilot got control over aerodynamic surfaces to change angle of attack. That changes CL. Here, we're going to make user (pilot) change CL value, directly, i.e. CL will be our control. To make it, we put an attribute CL_constant. Below it is the code you can append to previous code.

import math

class Aircraft(object):
    
    CD0 = 0
    taper_ratio = 0
    aspect_ratio = 0
    wing_area = 0
    mass = 0
    k = 0
    CL_constant = 1.0
    
    def __init__(self, CD0, taper_ratio, aspect_ratio, wing_area, mass):
        self.CD0 = CD0
        self.taper_ratio = taper_ratio
        self.aspect_ratio = aspect_ratio
        self.wing_area = wing_area
        self.mass = mass
        self.k = 1/(math.pi*aspect_ratio*taper_ratio);
        
    def CD(self,alfa=0.0):
        return self.CD0+self.k*self.CL(alfa)**2
    
    def CL(self,alfa=0.0):
        return self.CL_constant

Testing the code in python interpreter:

$ python -i gliding.py 
>>> a = Aircraft(CD0=0.0115, taper_ratio=0.94, \
                 aspect_ratio=17, wing_area=16.01, mass=512)
>>> a.CL()
1.0
>>> a.CD()
0.031419266970199665
>>> a.CL_constant = 1.5
>>> a.CL()
1.5
>>> a.CD()
0.056318350682949242

Flight class

This class uses some more specific object concepts as inheritance and composition. One class can inherits behavior from another class so that this subclass got the same methods of your superclass and, possibly, further methods. In our case, we have created a big class called Flight which can be reusable for many different kind of flights with different differential equation systems (ODES) modeling it. In Python you specify superclass for your class in parentheses following class name in class declaration. All previous classes inherits from object class (here you can see the reason for that). FlightSail class below inherits from Flight.

Flight class got methods for solving ODES, using SciPy library. We solve ODES by time slicing (says 1 second slice) and store ODES solution in a vector called state_history containing tuples in the form (t,y) where y is the vector representing one state of flight. In gliding flight y contains velocity, height, range, etc. Below we can see both classes Flight and subclass FlightSail.

import math
from scipy.integrate import ode
from scipy.integrate import odeint
from scipy.optimize import fsolve

class Flight(object):
    """Subclass this class and implement mech_ode_sys method
    """
    
    aircraft = None
    atmosphere = None
    
    def __init__(self,aircraft,atmosphere):
        self.aircraft = aircraft
        self.atmosphere = atmosphere
        
    def mech_ode_sys(self, t, state):
        """Must be implemented by subclass
        """    
    
    def start(self, initial_state, t_0 = 0, dt = 1):
        self.ode_solver = ode(self.mech_ode_sys).set_integrator('vode')
        self.ode_solver.set_initial_value(initial_state, t_0)
        self.dt = dt
        self.state_history = [[0.0,initial_state]]
    
    def next_state(self):
        if self.ode_solver.successful():
            self.ode_solver.integrate(self.ode_solver.t+self.dt)
            self.state_history.append((self.ode_solver.t, self.ode_solver.y,))
            return [self.ode_solver.t, self.ode_solver.y]
        else:
            return None
    
    def solve(self,t):
        """Solve ODE system for a given array of values t
        """
        def func(state,t):
            return self.mech_ode_sys(t, state)
        
        y0 = self.state_history[-1][1]
        y = odeint(func,y0,t)
        states = zip(t,y[:])
        self.state_history += states 
        return states

class FlightSail(Flight):
    
    VELOCITY = 0
    GAMMA = 1
    HEIGHT = 2
    RANGE = 3
    
    def __init__(self,aircraft,atmosphere):
        super(FlightSail,self).__init__(aircraft,atmosphere)
        
    def mech_ode_sys(self, t, state):
        """state is a tuple in the form (Velocity,gamma,Height,Range)
        """
        V = state[VELOCITY]
        gamma = state[GAMMA]
        H = state[HEIGHT]
        R = state[RANGE]
        
        CL = self.aircraft.CL()
        CD = self.aircraft.CD()
        m = self.aircraft.mass
        
        rho = self.atmosphere.rho(H)
        
        A = self.aircraft.wing_area
        
        L = 0.5*rho*V**2*A*CL
        D = 0.5*rho*V**2*A*CD
        
        V_p = (-D-m*self.atmosphere.g_0*math.sin(gamma))/m
        gamma_p = (L-m*self.atmosphere.g_0*math.cos(gamma))/(m*V)
        H_p = V*math.sin(gamma)
        R_p = V*math.cos(gamma)
        return [V_p, gamma_p, H_p, R_p]

Putting it all together

Now, instantiating the same problem from previous post.

import math

atmosphere = Atmosphere()
a = Aircraft(CD0=0.0115, taper_ratio=0.94, \
             aspect_ratio=17, wing_area=16.01, mass=512)
a.CL_constant = 0.75
sail = FlightSail(aircraft=a, atmosphere=atmosphere)
V_0 = 30
gamma_0 = -3*math.pi/180
H_0 = 600
R_0 = 0
state_0 = (V_0, gamma_0, H_0, R_0)
sail.start(initial_state=state_0)
aircraft_state = list(state_0)
MAX_RANGE = 14000

For using PyGame we have put a profile aircraft figure (in our case Super Tucano) and just update its position in a 640 x 480 screen, let CL ranging according up and down keys, and adjusting path angle, rotating the figure (this is not correct, but we are neglecting angle of attack here). PyGame code is shown below:

pygame.init()
SCREEN_X_MAX = 640
SCREEN_Y_MAX = 480
screen = pygame.display.set_mode((SCREEN_X_MAX,SCREEN_Y_MAX))
tucano_img_original = pygame.image.\
                        load('supertucano_lateral_menor.bmp').convert()

pygame.display.update()
tucano_img = None 
clock = pygame.time.Clock()
clock.tick(10)
start = False

while True:
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            sys.exit()
        elif event.type == pygame.KEYDOWN:
            start = True
    
    key = pygame.key.get_pressed()
    if key[pygame.K_UP]:
        sail.aircraft.CL_constant += 0.1
        print "CL increasing to: "+str(sail.aircraft.CL_constant)
    elif key[pygame.K_DOWN]:
        sail.aircraft.CL_constant -= 0.1
        print "CL decreasing to: "+str(sail.aircraft.CL_constant)
    
    if start:
        state = sail.next_state()
        #print state
        pygame.display.set_caption('Sail Flight t = '+str(state[0]))
        aircraft_state = [i for i in state[1]]
        x = aircraft_state[FlightSail.RANGE]/float(MAX_RANGE)\
                         *SCREEN_X_MAX
        y = SCREEN_Y_MAX - float(aircraft_state[FlightSail.HEIGHT])\
                         /H_0*SCREEN_Y_MAX
        tucano_img = pygame.transform.rotate(tucano_img_original,\
                         aircraft_state[FlightSail.GAMMA]/math.pi*180)
        screen.fill((255,255,255))
        screen.blit(tucano_img,(x,y))
        pygame.display.update()
        pygame.time.delay(100)

Figure below shows whole code running (using python gliding.py at command line).



And that is it! Keys UP and DOWN changes CL directly, what is not quiet representative but it is fun though.