Simple Physics-based Flight Simulation with C++

A picture of the author, Jakob Maier
Jakob Maier
Jan 8, 2023



In this blog post I aim to show you how to create a reasonably realistic flight simulator in C++ from scratch, so no physics engine needed.

A little disclaimer: I am by no means a physics expert, so take everything you read here with a grain of salt. What I show here is what I have learnt from reading various books on the subject, so please do not use this post to create a flight simulator to be used for pilot training. 😉

Simulating a Generic Rigid Body

Before we can get into programming a flight model we start with creating a generic rigid body. Rigid bodies are used to simulate real-time physics for objects that do not deform under forces (they are 'rigid').

We would like to have an object that we can apply abitrary forces to and that then reacts based on its current velocity and angular velocity, as well as it's mass and inertia.

One important concept we need to understand is 'world coordinates' vs 'body coordinates'. World coordinates are exactly what they sound like, the represent the position of our rigidbody in 3D space. I work with OpenGL, so for me the Y-coordinate is up, X points to the right and Z is the direction out from the screen. Any rigidbody will have both a position as well as an orientation. Body coordinates on the other hand are in reference to the rigidbody itself. Think of it like this: if you are sitting in the airplane, to your left there is always the left wing and if you look up you will see the 'roof' of the airplane, even if the airplane is actually maneuvering in some other way. Body coordinates are very useful and we will use it to calculate our lift and drag forces as well as the resulting torque.

class RigidBody {
        glm::vec3 m_force{}; // world space
        glm::vec3 m_torque{}; // body space

        float mass = 1.0f;                                  // kg
        glm::vec3 position{};                               // world space
        glm::quat orientation{};                            // world space
        glm::vec3 velocity{};                               // world space, meter/second
        glm::vec3 angular_velocity{};                       // body space, radians/second
        glm::mat3 inertia{}, inverse_inertia{};             // inertia tensor, body space
        bool apply_gravity = true;

        // transform direction from body space to world space
        inline glm::vec3 transform_direction(const glm::vec3& direction) const
        { return orientation * direction; }

        // transform direction from world space to body space
        inline glm::vec3 inverse_transform_direction(const glm::vec3& direction) const
        { return glm::inverse(orientation) * direction;}

        // get velocity and angular velocity in body space
        inline glm::vec3 get_point_velocity(const glm::vec3& point) const
        { return inverse_transform_direction(velocity) + glm::cross(angular_velocity, point); }

        // force and point vectors are in body space
        inline void add_force_at_point(const glm::vec3& force, const glm::vec3& point)
        { m_force += transform_direction(force), m_torque += glm::cross(point, force); }

        // force vector in body space
        inline void add_relative_force(const glm::vec3& force)
        { m_force += transform_direction(force); }

        // integrate using the euler method
        virtual void update(float dt)
            // integrate position
            glm::vec3 acceleration = m_force / mass;
            if (apply_gravity) acceleration.y -= 9.81f;
            velocity += acceleration * dt;
            position += velocity * dt;
            // integrate orientation
            angular_velocity += inverse_inertia *
                (m_torque - glm::cross(angular_velocity, inertia * angular_velocity)) * dt;
            orientation += (orientation * glm::quat(0.0f, angular_velocity)) * (0.5f * dt);
            orientation = glm::normalize(orientation);

            // reset accumulators
            m_force = glm::vec3(0.0f), m_torque = glm::vec3(0.0f);

Creating the Flightmodel

We simulate our aircraft by splitting it up in a number of aerodynamic surfaces, which we will call wings for simplicity, but they can also be control surfaces and the fuselage. The other part of our aircraft that produces a force is of course the engine.

struct Airplane : public phi::RigidBody {
    Engine engine;
    std::vector<Wing> elements;
    void update(float dt) override

        for (auto& wing : elements)


Wings create lift and drag depending on the speed of the air over them as well as the angle between the wing and the air flow. This angle is called 'angle of attack' or 'alpha'.

The lift and drag force are calculated using these formulas:



Total force

Where v is the velocity, A is the wing area, ρ is the air pressure and Cl/Cd are the coefficients of lift and drag.

Coefficient of lift and coefficient of drag depend on the shape (airfoil) of the wing and the angle of attack. Data on standard airfoils is available online, I got the data from As you can see, we sample these coefficients in Wing from airfoil.

Airplanes manouver by deflecting their control surfaces, which then produces differential lift, creating a torque on the aircraft. We do this by changing the wing normal depending on the deflection of the control surface.

  Wing(const glm::vec3& position, float span, float chord, const Airfoil* airfoil, const glm::vec3& normal,
       float flap_ratio = 0.25f)
      : airfoil(airfoil),
        area(span * chord),
        aspect_ratio(phi::sq(span) / area),

  // controls how much the wing is deflected
  void set_control_input(float input) { control_input = glm::clamp(input, -1.0f, 1.0f); }

  // compute and apply aerodynamic forces
  void apply_forces(phi::RigidBody* rigid_body, phi::Seconds dt)
    glm::vec3 local_velocity = rigid_body->get_point_velocity(center_of_pressure);
    float speed = glm::length(local_velocity);

    if (speed <= 1.0f) return;

    // drag acts in the opposite direction of velocity
    glm::vec3 drag_direction = glm::normalize(-local_velocity);

    // lift is always perpendicular to drag
    glm::vec3 lift_direction = glm::normalize(glm::cross(glm::cross(drag_direction, normal), drag_direction));

    // angle between chord line and air flow
    float angle_of_attack = glm::degrees(std::asin(glm::dot(drag_direction, normal)));

    // sample aerodynamic coefficients
    auto [lift_coeff, drag_coeff] = airfoil->sample(angle_of_attack);

    if (flap_ratio > 0.0f) {
      float deflection_ratio = control_input;

      // lift coefficient changes based on flap deflection
      float delta_lift_coeff = sqrt(flap_ratio) * airfoil->cl_max * deflection_ratio;
      lift_coeff += delta_lift_coeff;

    // induced drag, increases with lift
    float induced_drag_coeff = phi::sq(lift_coeff) / (phi::PI * aspect_ratio * efficiency_factor);
    drag_coeff += induced_drag_coeff;

    // air density depends on altitude
    float air_density = isa::get_air_density(rigid_body->position.y);

    float dynamic_pressure = 0.5f * phi::sq(speed) * air_density * area;

    glm::vec3 lift = lift_direction * lift_coeff * dynamic_pressure;
    glm::vec3 drag = drag_direction * drag_coeff * dynamic_pressure;

    // aerodynamic forces are applied at the center of pressure
    rigid_body->add_force_at_point(lift + drag, center_of_pressure);

To get our aerodynamic coefficients we simply sample from the array of data points we got from airfoiltools:

// NACA 0012: alpha, Cl, Cd
std::vector<glm::vec3> NACA_0012_data = {
    {-18.500f, -1.2258f, 0.10236f},
    {18.500f, 1.2284f, 0.10229f}

struct Airfoil {
    const float min_alpha, max_alpha;
    std::vector<glm::vec3> data;

    Airfoil(const std::vector<glm::vec3> &curve) : data(curve)
        min_alpha = curve[0].x, max_alpha = curve[curve.size() - 1].x;

    // get Cl and Cd
    std::tuple<float, float> sample(float alpha) const
        int i = static_cast<int>(scale(alpha, min_alpha, max_alpha, 0, data.size() - 1));
        return { data[i].y, data[i].z };

Our airplane of course needs an engine, which we simulate by simply applying a force in the forward direction directly to the center of gravity of the aircraft:

struct Engine {
    float throttle = 1.0f, thrust;
    Engine(float thrust) : thrust(thrust) {}
    void apply_force(phi::RigidBody *rigid_body)
        // thrust is applied to the center of gravity and does not produce torque
        rigid_body->add_relative_force(phi::FORWARD * (throttle * thrust));

Putting it all together

One of the hardest parts of doing simulations like this is plugging in the correct values to make it look reasonable. Messing up the wing positions and area does not make a huge difference, but a wrong inertia tensor quickly makes the simulation spin out of control. So here is a setup that ended up working for me:

const float mass = 10000.0f;
const float thrust = 50000.0f;

const float wing_offset = -1.0f;
const float tail_offset = -6.6f;

const Airfoil NACA_0012(NACA_0012_data);
const Airfoil NACA_2412(NACA_2412_data);

std::vector<Wing> wings = {
      Wing({wing_offset, 0.0f, -2.7f}, 6.96f, 2.50f, &NACA_2412),             // left wing
      Wing({wing_offset - 1.5f, 0.0f, -2.0f}, 3.80f, 1.26f, &NACA_0012),      // left aileron
      Wing({wing_offset - 1.5f, 0.0f, 2.0f}, 3.80f, 1.26f, &NACA_0012),       // right aileron
      Wing({wing_offset, 0.0f, +2.7f}, 6.96f, 2.50f, &NACA_2412),             // right wing
      Wing({tail_offset, -0.1f, 0.0f}, 6.54f, 2.70f, &NACA_0012),             // elevator
      Wing({tail_offset, 0.0f, 0.0f}, 5.31f, 3.10f, &NACA_0012, phi::RIGHT),  // rudder

glm::mat3 inertia = {
    48531.0f,-1320.0f, 0.0f,
   -1320.0f, 256608.0f, 0.0f,
    0.0f, 0.0f, 211333.0f

Airplane airplane(mass, thrust, inertia, wings);

airplane.position = { 0.0f, 2000.0f, 0.0f };
airplane.velocity = { phi::units::meter_per_second(600.0f), 0.0f, 0.0f };

As you can see, the results of our simple flightmodel look quite realistic:

Nevermind the ugly graphics, this is about the flightmodel. Generating pretty terrain is a topic for a future blog post...

The code can be found here. If you have any suggestions, questions or problems, don't hesitate to write me an email or open an issue on github.


↑ back to top

© 2023 Jakob Maier
kontakt & impressum