diff --git a/.clang-format b/.clang-format
index c0c314a..f4cb2f3 100644
--- a/.clang-format
+++ b/.clang-format
@@ -3,8 +3,8 @@ Language:        Cpp
 # BasedOnStyle:  Google
 AccessModifierOffset: -1
 AlignAfterOpenBracket: Align
-AlignConsecutiveMacros: false
-AlignConsecutiveAssignments: false
+AlignConsecutiveMacros: true
+AlignConsecutiveAssignments: true
 AlignConsecutiveDeclarations: false
 AlignEscapedNewlines: Left
 AlignOperands:   true
diff --git a/OpenGL_Flightsim/CMakeLists.txt b/OpenGL_Flightsim/CMakeLists.txt
index 0392e15..681ea83 100644
--- a/OpenGL_Flightsim/CMakeLists.txt
+++ b/OpenGL_Flightsim/CMakeLists.txt
@@ -11,7 +11,7 @@ add_executable(flightsim
     # Source files
     src/ai.h
     src/clipmap.h
-    src/collisions.h
+    src/collider.h
     src/data.h
     src/flightmodel.h
     src/gfx.cpp
diff --git a/OpenGL_Flightsim/OpenGL_Flightsim.vcxproj b/OpenGL_Flightsim/OpenGL_Flightsim.vcxproj
index bd16b31..ada39bd 100644
--- a/OpenGL_Flightsim/OpenGL_Flightsim.vcxproj
+++ b/OpenGL_Flightsim/OpenGL_Flightsim.vcxproj
@@ -177,7 +177,7 @@ xcopy /y /i /s $(ProjectDir)..\lib\glew-2.2.0\bin\Release\x64\*.dll $(OutDir)</C
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="src\ai.h" />
-    <ClInclude Include="src\collisions.h" />
+    <ClInclude Include="src\collider.h" />
     <ClInclude Include="src\data.h" />
     <ClInclude Include="src\clipmap.h" />
     <ClInclude Include="src\flightmodel.h" />
diff --git a/OpenGL_Flightsim/OpenGL_Flightsim.vcxproj.filters b/OpenGL_Flightsim/OpenGL_Flightsim.vcxproj.filters
index bfdd08e..07cc736 100644
--- a/OpenGL_Flightsim/OpenGL_Flightsim.vcxproj.filters
+++ b/OpenGL_Flightsim/OpenGL_Flightsim.vcxproj.filters
@@ -81,7 +81,7 @@
     <ClInclude Include="src\data.h">
       <Filter>Header Files</Filter>
     </ClInclude>
-    <ClInclude Include="src\collisions.h">
+    <ClInclude Include="src\collider.h">
       <Filter>Header Files</Filter>
     </ClInclude>
     <ClInclude Include="lib\imgui\imgui.h">
diff --git a/OpenGL_Flightsim/src/ai.h b/OpenGL_Flightsim/src/ai.h
index 394fdd5..6f308de 100644
--- a/OpenGL_Flightsim/src/ai.h
+++ b/OpenGL_Flightsim/src/ai.h
@@ -11,28 +11,28 @@
 glm::vec3 get_intercept_point(const glm::vec3& position, const glm::vec3& velocity, const glm::vec3& target_position,
                               const glm::vec3& target_velocity)
 {
-  auto velocity_delta = target_velocity - velocity;
-  auto position_delta = target_position - position;
+  auto velocity_delta    = target_velocity - velocity;
+  auto position_delta    = target_position - position;
   auto time_to_intercept = glm::length(position_delta) / glm::length(velocity_delta);
   return target_position + target_velocity * time_to_intercept;
 }
 
 void fly_towards(Airplane& airplane, const glm::vec3& target)
 {
-  auto& rb = airplane;
+  auto& rb       = airplane;
   auto& joystick = airplane.joystick;
-  auto position = rb.position;
+  auto position  = rb.position;
   auto direction = glm::normalize(rb.inverse_transform_direction(target - rb.position));
-  auto angle = glm::angle(phi::FORWARD, direction);
+  auto angle     = glm::angle(phi::FORWARD, direction);
 
-  float rudder = direction.z;
+  float rudder   = direction.z;
   float elevator = direction.y * 5.0f;
 
-  float m = phi::PI / 4.0f;
-  float agressive_roll = direction.z;
-  float wings_level_roll = rb.right().y;
+  float m                     = phi::PI / 4.0f;
+  float agressive_roll        = direction.z;
+  float wings_level_roll      = rb.right().y;
   float wings_level_influence = phi::inverse_lerp(0.0f, m, glm::clamp(angle, -m, m));
-  float aileron = phi::lerp(wings_level_roll, agressive_roll, wings_level_influence);
+  float aileron               = phi::lerp(wings_level_roll, agressive_roll, wings_level_influence);
 
   joystick = glm::clamp(glm::vec4(aileron, rudder, elevator, 0.0f), glm::vec4(-1.0f), glm::vec4(1.0f));
 }
diff --git a/OpenGL_Flightsim/src/clipmap.h b/OpenGL_Flightsim/src/clipmap.h
index 01f39e1..fb7b54a 100644
--- a/OpenGL_Flightsim/src/clipmap.h
+++ b/OpenGL_Flightsim/src/clipmap.h
@@ -3,7 +3,7 @@
 #include "gfx.h"
 
 constexpr unsigned int primitive_restart = 0xFFFFU;
-const std::string path = "assets/textures/terrain/1/";
+const std::string path                   = "assets/textures/terrain/1/";
 
 struct Seam {
   gfx::gl::VertexBuffer vbo;
@@ -12,7 +12,7 @@ struct Seam {
 
   Seam(int columns, float size)
   {
-    int rows = 1;
+    int rows    = 1;
     index_count = columns;
 
     std::vector<glm::vec3> vertices;
@@ -132,8 +132,8 @@ class Clipmap : public gfx::Object3D
   {
 #if 1
     if (!context.is_shadow_pass) {
-      auto camera_pos = context.camera->get_world_position();
-      float height = camera_pos.y;
+      auto camera_pos    = context.camera->get_world_position();
+      float height       = camera_pos.y;
       auto camera_pos_xy = glm::vec2(camera_pos.x, camera_pos.z);
 
       heightmap.bind(2);
@@ -159,12 +159,12 @@ class Clipmap : public gfx::Object3D
       for (int l = min_level; l <= levels; l++) {
         const int rows = 5, cols = 5;
         // float border = 0.0f;
-        float scale = std::pow(2.0f, l);
-        float next_scale = std::pow(2.0f, l + 2);
+        float scale               = std::pow(2.0f, l);
+        float next_scale          = std::pow(2.0f, l + 2);
         float scaled_segment_size = segment_size * scale;
-        float tile_size = segments * scaled_segment_size;
-        glm::vec2 snapped = glm::floor(camera_pos_xy / next_scale) * next_scale;
-        auto base = calc_base(l, camera_pos_xy);
+        float tile_size           = segments * scaled_segment_size;
+        glm::vec2 snapped         = glm::floor(camera_pos_xy / next_scale) * next_scale;
+        auto base                 = calc_base(l, camera_pos_xy);
 
         shader.uniform("u_Scale", scale);
         shader.uniform("u_SegmentSize", scaled_segment_size);
@@ -183,7 +183,7 @@ class Clipmap : public gfx::Object3D
           center.draw();
         } else {
           auto prev_base = calc_base(l - 1, camera_pos_xy);
-          auto diff = glm::abs(base - prev_base);
+          auto diff      = glm::abs(base - prev_base);
 
           auto l_offset = glm::vec2(tile_size, tile_size);
           if (diff.x == tile_size) {
@@ -283,11 +283,11 @@ class Clipmap : public gfx::Object3D
 
   glm::vec2 calc_base(int level, glm::vec2 camera_pos)
   {
-    float scale = std::pow(2.0f, level);
-    float next_scale = std::pow(2.0f, level + 2);
-    float tile_size = segments * segment_size * scale;
+    float scale       = std::pow(2.0f, level);
+    float next_scale  = std::pow(2.0f, level + 2);
+    float tile_size   = segments * segment_size * scale;
     glm::vec2 snapped = glm::floor(camera_pos / next_scale) * next_scale;
-    glm::vec2 base = snapped - tile_size * 2.0f;
+    glm::vec2 base    = snapped - tile_size * 2.0f;
     return base;
   }
 
diff --git a/OpenGL_Flightsim/src/collisions.h b/OpenGL_Flightsim/src/collider.h
similarity index 79%
rename from OpenGL_Flightsim/src/collisions.h
rename to OpenGL_Flightsim/src/collider.h
index 2026b7a..a461c25 100644
--- a/OpenGL_Flightsim/src/collisions.h
+++ b/OpenGL_Flightsim/src/collider.h
@@ -3,26 +3,13 @@
 */
 #pragma once
 
-#include <glm/gtc/matrix_transform.hpp>
-#include <glm/gtx/euler_angles.hpp>
-#include <glm/gtx/vector_angle.hpp>
-#include <glm/mat4x4.hpp>
 #include <glm/vec3.hpp>
-#include <tuple>
 
 #include "phi.h"
 
 namespace col
 {
 
-constexpr float EPSILON = 1e-8f;
-
-template <typename T>
-constexpr inline T sq(T x)
-{
-  return x * x;
-}
-
 typedef void CollisionCallback(const glm::vec3& point, const glm::vec3& normal);
 
 struct Collider {
@@ -45,13 +32,15 @@ struct AABB : public Collider {
 struct Sphere : public Collider {
   glm::vec3 center;
   float radius;
+  Sphere() : Sphere(glm::vec3(0.0f), 1.0f) {}
   Sphere(const glm::vec3& center, float radius) : center(center), radius(radius) {}
 };
 
-// ray = origin + direction * t
 struct Ray : public Collider {
   glm::vec3 origin, direction;
+  Ray() : Ray(glm::vec3(0.0f), glm::vec3(1.0f, 0.0f, 0.0f)) {}
   Ray(const glm::vec3& origin, const glm::vec3& direction) : origin(origin), direction(direction) {}
+  inline glm::vec3 point_at(float t) const { return origin + direction * t; }
 };
 
 struct Heightmap {
@@ -85,9 +74,9 @@ struct Heightmap {
   float get_height(const glm::vec2& coord) const
   {
     glm::vec2 tmp = glm::clamp(coord / magnification, glm::vec2(-1.0f), glm::vec2(1.0f));
-    auto uv = phi::scale(tmp, glm::vec2(-1.0f), glm::vec2(1.0f), glm::vec2(0.0f), glm::vec2(1.0f));
-    float value = sample(uv).r;
-    float height = scale * value + shift;
+    auto uv       = phi::scale(tmp, glm::vec2(-1.0f), glm::vec2(1.0f), glm::vec2(0.0f), glm::vec2(1.0f));
+    float value   = sample(uv).r;
+    float height  = scale * value + shift;
     return height;
   }
 };
@@ -96,15 +85,13 @@ struct Heightmap {
 bool test_collision(const Ray& r, const Sphere& s, float* t)
 {
   // Christer_Ericson-Real-Time_Collision_Detection.pdf#page=178
-  assert(std::abs(glm::length(r.direction) - 1.0f) < EPSILON);
-
   auto m = r.origin - s.center;
   auto b = glm::dot(m, r.direction);
-  auto c = glm::dot(m, m) - sq(s.radius);
+  auto c = glm::dot(m, m) - phi::sq(s.radius);
 
   if (c > 0.0f && b > 0.0f) return false;
 
-  auto discr = sq(b) - c;
+  auto discr = phi::sq(b) - c;
 
   if (discr < 0.0f) return false;
 
@@ -113,11 +100,19 @@ bool test_collision(const Ray& r, const Sphere& s, float* t)
 }
 
 // test collision between two spheres
-bool test_collision(const Sphere& s0, const Sphere& s1)
+bool test_collision(const Sphere& a, const Sphere& b, phi::CollisionInfo* info)
 {
-  float distance = glm::length(s0.center - s1.center);
-  float radius_sum = s0.radius + s1.radius;
-  return distance < radius_sum;
+  float distance   = glm::length(a.center - b.center);
+  float radius_sum = a.radius + b.radius;
+
+  if (distance < radius_sum) {
+    info->normal      = glm::normalize(b.center - a.center);
+    info->penetration = radius_sum - distance;
+    info->point       = a.center + a.radius * info->normal;
+    return true;
+  } else {
+    return false;
+  }
 }
 
 // test collision between two axis aligned bounding boxes
@@ -129,9 +124,10 @@ bool test_collision(const AABB& a, const AABB& b)
           (a_max.z < b_min.z || a_min.z > b_max.z));
 }
 
-bool test_collision(const Heightmap& heightmap, const glm::vec3& point)
+bool test_collision(const Heightmap& heightmap, const glm::vec3& point, float* height)
 {
-  return point.y <= heightmap.get_height({point.x, point.z});
+  *height = heightmap.get_height({point.x, point.z});
+  return point.y <= *height;
 }
 
 // test collision of two moving spheres
@@ -171,8 +167,8 @@ bool test_moving_collision(const Sphere& s0, const glm::vec3& velocity0, const S
 #else
   // Christer_Ericson-Real-Time_Collision_Detection.pdf#page=226
   float tmp_t = 0.0f;
-  auto v = velocity0 - velocity1;
-  auto vlen = glm::length(v);
+  auto v      = velocity0 - velocity1;
+  auto vlen   = glm::length(v);
 
   Ray ray(s0.center, v / vlen);
   Sphere sphere(s1.center, s0.radius + s1.radius);
diff --git a/OpenGL_Flightsim/src/flightmodel.h b/OpenGL_Flightsim/src/flightmodel.h
index 9283f27..1446d17 100644
--- a/OpenGL_Flightsim/src/flightmodel.h
+++ b/OpenGL_Flightsim/src/flightmodel.h
@@ -27,39 +27,43 @@ float get_air_density(float altitude)
 {
   assert(0.0f <= altitude && altitude <= 11000.0f);
   float temperature = get_air_temperature(altitude);
-  float pressure = 101325.0f * std::pow(1 - 0.0065f * (altitude / 288.15f), 5.25f);
+  float pressure    = 101325.0f * std::pow(1 - 0.0065f * (altitude / 288.15f), 5.25f);
   return 0.00348f * (pressure / temperature);
 }
 
 const float sea_level_air_density = get_air_density(0.0f);
 };  // namespace isa
 
-typedef glm::vec3 AeroData;  // AoA, Cl, Cd
+// AoA, Cl, Cd
+typedef glm::vec3 AeroData;
 
+// aerodynamic data sampler
 struct Airfoil {
   float min_alpha, max_alpha;
   std::vector<AeroData> data;
 
   Airfoil(const std::vector<AeroData>& curve) : data(curve) { min_alpha = curve.front().x, max_alpha = curve.back().x; }
 
-  // lift_coeff, drag_coeff, moment_coeff
+  // lift_coeff, drag_coeff
   std::tuple<float, float> sample(float alpha) const
   {
-    int max_index = static_cast<int>(data.size() - 1);
-    float t = phi::inverse_lerp(min_alpha, max_alpha, alpha) * max_index;
-    float integer = std::floor(t);
+    int max_index    = static_cast<int>(data.size() - 1);
+    float t          = phi::inverse_lerp(min_alpha, max_alpha, alpha) * max_index;
+    float integer    = std::floor(t);
     float fractional = t - integer;
-    int index = static_cast<int>(integer);
-    auto value = (index < max_index) ? phi::lerp(data[index], data[index + 1], fractional) : data[max_index];
+    int index        = static_cast<int>(integer);
+    auto value       = (index < max_index) ? phi::lerp(data[index], data[index + 1], fractional) : data[max_index];
     return {value.y, value.z};
   }
 };
 
+// base engine
 struct Engine : public phi::ForceGenerator {
   float throttle = 0.25f;
   void apply_forces(phi::RigidBody* rigid_body, phi::Seconds dt) override {}
 };
 
+// simple jet-like engine
 struct SimpleEngine : public Engine {
   const float thrust;
   SimpleEngine(float thrust) : thrust(thrust) {}
@@ -70,6 +74,7 @@ struct SimpleEngine : public Engine {
   }
 };
 
+// does not yet implement engine torque
 struct PropellorEngine : public Engine {
   float horsepower, rpm, propellor_diameter;
 
@@ -80,18 +85,18 @@ struct PropellorEngine : public Engine {
 
   void apply_forces(phi::RigidBody* rigid_body, phi::Seconds dt) override
   {
-    float speed = rigid_body->get_speed();
-    float altitude = rigid_body->position.y;
+    float speed        = rigid_body->get_speed();
+    float altitude     = rigid_body->position.y;
     float engine_power = phi::units::watts(horsepower);
 
     const float a = 1.83f, b = -1.32f;  // efficiency curve fit coefficients
-    float turnover_rate = rpm / 60.0f;
+    float turnover_rate           = rpm / 60.0f;
     float propellor_advance_ratio = speed / (turnover_rate * propellor_diameter);
-    float propellor_efficiency = a * propellor_advance_ratio + b * std::pow(propellor_advance_ratio, 3);
+    float propellor_efficiency    = a * propellor_advance_ratio + b * std::pow(propellor_advance_ratio, 3);
     assert(0.0f <= propellor_efficiency && propellor_efficiency <= 1.0f);
 
-    const float c = 0.12f;  // mechanical power loss factor
-    float air_density = isa::get_air_density(altitude);
+    const float c               = 0.12f;  // mechanical power loss factor
+    float air_density           = isa::get_air_density(altitude);
     float power_drop_off_factor = ((air_density / isa::sea_level_air_density) - c) / (1 - c);
     assert(0.0f <= power_drop_off_factor && power_drop_off_factor <= 1.0f);
 
@@ -101,6 +106,7 @@ struct PropellorEngine : public Engine {
   }
 };
 
+// not only a wing, can be any kind of aerodynamic surface
 class Wing : public phi::ForceGenerator
 {
  private:
@@ -112,19 +118,19 @@ class Wing : public phi::ForceGenerator
   const glm::vec3 normal;
   const glm::vec3 center_of_pressure;
 
-  float lift_multiplier = 1.0f;
-  float drag_multiplier = 1.0f;
+  float lift_multiplier   = 1.0f;
+  float drag_multiplier   = 1.0f;
   float efficiency_factor = 1.0f;
 
-  float deflection = 0.0f;
-  float control_input = 0.0f;
-  float min_deflection = -10.0f;
-  float max_deflection = +10.0f;
-  float max_actuator_speed = 90.0f;
+  float deflection          = 0.0f;
+  float control_input       = 0.0f;
+  float min_deflection      = -10.0f;
+  float max_deflection      = +10.0f;
+  float max_actuator_speed  = 90.0f;
   float max_actuator_torque = 6000.0f;
 
  public:
-  float incidence = 0.0f;
+  float incidence         = 0.0f;
   bool is_control_surface = true;
 
   Wing(const Airfoil* airfoil, const glm::vec3& position, float area, float span, const glm::vec3& normal = phi::UP)
@@ -159,7 +165,7 @@ class Wing : public phi::ForceGenerator
   void apply_forces(phi::RigidBody* rigid_body, phi::Seconds dt) override
   {
     glm::vec3 local_velocity = rigid_body->get_point_velocity(center_of_pressure);
-    float speed = glm::length(local_velocity);
+    float speed              = glm::length(local_velocity);
 
     if (speed <= phi::EPSILON) return;
 
@@ -185,8 +191,8 @@ class Wing : public phi::ForceGenerator
     float air_density = isa::get_air_density(0.0f);  // something is not right here, so let's assume sea level
 
     float dynamic_pressure = 0.5f * std::pow(speed, 2) * air_density * area;
-    glm::vec3 lift = lift_direction * lift_coeff * lift_multiplier * dynamic_pressure;
-    glm::vec3 drag = drag_direction * (drag_coeff + induced_drag_coeff) * drag_multiplier * dynamic_pressure;
+    glm::vec3 lift         = lift_direction * lift_coeff * lift_multiplier * dynamic_pressure;
+    glm::vec3 drag         = drag_direction * (drag_coeff + induced_drag_coeff) * drag_multiplier * dynamic_pressure;
 
     // aerodynamic forces are applied at the center of pressure
     rigid_body->add_force_at_point(lift + drag, center_of_pressure);
@@ -209,25 +215,27 @@ class Wing : public phi::ForceGenerator
     deflection = (control_input >= 0.0f ? max_deflection : min_deflection) * std::abs(control_input);
 #endif
 
-    auto axis = glm::normalize(glm::cross(phi::FORWARD, normal));
+    auto axis     = glm::normalize(glm::cross(phi::FORWARD, normal));
     auto rotation = glm::rotate(glm::mat4(1.0f), glm::radians(incidence + deflection), axis);
     return glm::vec3(rotation * glm::vec4(normal, 1.0f));
   }
 };
 
+// simple flightmodel
 struct Airplane : public phi::RigidBody {
   glm::vec4 joystick{};  // roll, yaw, pitch, elevator trim
   Engine* engine;
   std::vector<Wing> surfaces;
+  bool is_landed = false;
 
 #if LOG_FLIGHT
-  float log_timer = 0.0f;
+  float log_timer     = 0.0f;
   float log_intervall = 0.1f;
-  float flight_time = 0.0f;
+  float flight_time   = 0.0f;
   std::ofstream log_file;
 #endif
 
-  Airplane(float mass, Engine* engine, glm::mat3 inertia, std::vector<Wing> elements)
+  Airplane(float mass, const glm::mat3& inertia, std::vector<Wing> elements, Engine* engine)
       : phi::RigidBody({.mass = mass, .inertia = inertia}), surfaces(elements), engine(engine)
   {
 #if LOG_FLIGHT
@@ -243,9 +251,10 @@ struct Airplane : public phi::RigidBody {
 
     surfaces[1].set_deflection_limits(-15.0f, 15.0f);
     surfaces[2].set_deflection_limits(-15.0f, 15.0f);
+    surfaces[5].set_deflection_limits(-5.0f, 5.0f);
   }
 
-  void update_flightmodel(phi::Seconds dt)
+  void update(phi::Seconds dt) override
   {
     float aileron = joystick.x, rudder = joystick.y, elevator = joystick.z, trim = joystick.w;
 
@@ -258,8 +267,8 @@ struct Airplane : public phi::RigidBody {
 #if LOG_FLIGHT
     flight_time += dt;
     if ((log_timer -= dt) <= 0.0f) {
-      log_timer = log_intervall;
-      auto av = glm::degrees(rigid_body.angular_velocity);
+      log_timer  = log_intervall;
+      auto av    = glm::degrees(rigid_body.angular_velocity);
       auto euler = glm::degrees(glm::eulerAngles(glm::normalize(rigid_body.orientation)));
 
       /* clang-format off */
@@ -287,13 +296,26 @@ struct Airplane : public phi::RigidBody {
       wing.apply_forces(this, dt);
     }
 
+    if (is_landed) {
+      // std::cout << "contact\n";
+      // calculate friction with ground
+      const float static_friction_coeff  = 0.2f;
+      const float kinetic_friction_coeff = 0.55f;
+      if (get_speed() > phi::EPSILON) {
+        add_friction(phi::UP, glm::normalize(velocity), kinetic_friction_coeff);
+      }
+    } else {
+      // std::cout << "no contact\n";
+    }
+
     engine->apply_forces(this, dt);
 
-    update(dt);
+    phi::RigidBody::update(dt);
   }
 
   float get_altitude() const { return position.y; }
 
+  // pitch g force
   float get_g() const
   {
     auto velocity = get_body_velocity();
@@ -312,23 +334,26 @@ struct Airplane : public phi::RigidBody {
     return g_force;
   }
 
+  // mach number
   float get_mach() const
   {
-    float temperature = isa::get_air_temperature(get_altitude());
+    float temperature    = isa::get_air_temperature(get_altitude());
     float speed_of_sound = std::sqrt(1.402f * 286.f * temperature);
     return get_speed() / speed_of_sound;
   }
 
+  // angle of attack
   float get_aoa() const
   {
     auto velocity = get_body_velocity();
     return glm::degrees(std::asin(glm::dot(glm::normalize(-velocity), phi::UP)));
   }
 
+  // indicated air speed
   float get_ias() const
   {
     // See: https://aerotoolbox.com/airspeed-conversions/
-    float air_density = isa::get_air_density(get_altitude());
+    float air_density      = isa::get_air_density(get_altitude());
     float dynamic_pressure = 0.5f * std::pow(get_speed(), 2) * air_density;  // bernoulli's equation
     return std::sqrt(2 * dynamic_pressure / isa::sea_level_air_density);
   }
diff --git a/OpenGL_Flightsim/src/gfx.cpp b/OpenGL_Flightsim/src/gfx.cpp
index 3dbc07f..019a82a 100644
--- a/OpenGL_Flightsim/src/gfx.cpp
+++ b/OpenGL_Flightsim/src/gfx.cpp
@@ -27,7 +27,7 @@ Shader::Shader(const std::string& path) : Shader(load_text_file(path + ".vert"),
 Shader::Shader(const std::string& vertShader, const std::string& fragShader)
 {
   // std::cout << "create Shader\n";
-  const char* vertexShaderSource = vertShader.c_str();
+  const char* vertexShaderSource   = vertShader.c_str();
   const char* fragmentShaderSource = fragShader.c_str();
 
   unsigned int vertexShader = glCreateShader(GL_VERTEX_SHADER);
@@ -341,13 +341,13 @@ glm::vec3 Object3D::get_scale() const { return m_scale; }
 
 void Object3D::set_scale(const glm::vec3& scale)
 {
-  m_scale = scale;
+  m_scale     = scale;
   m_dirty_dof = true;
 }
 
 void Object3D::set_position(const glm::vec3& pos)
 {
-  m_position = pos;
+  m_position  = pos;
   m_dirty_dof = true;
 }
 
@@ -366,7 +366,7 @@ void Object3D::set_transform(const glm::vec3& position, const glm::quat& rotatio
 
 void Object3D::set_rotation(const glm::vec3& rot)
 {
-  m_rotation = glm::quat(rot);
+  m_rotation  = glm::quat(rot);
   m_dirty_dof = true;
 }
 
@@ -394,8 +394,8 @@ void Object3D::traverse(const std::function<bool(Object3D*)>& func)
 void Object3D::override_transform(const glm::mat4& matrix)
 {
   m_dirty_transform = true;
-  m_dirty_dof = true;
-  transform = matrix;
+  m_dirty_dof       = true;
+  transform         = matrix;
   // TODO: relcalculate position, rotation etc
 }
 
@@ -464,8 +464,8 @@ Object3D::Type Light::get_type() const { return Object3D::Type::LIGHT; }
 glm::mat4 Light::light_space_matrix()
 {
   float near_plane = 0.1f, far_plane = 10.0f, m = 10.0f;
-  auto wp = get_world_position();
-  glm::mat4 light_view = glm::lookAt(wp, glm::vec3(0.0f), glm::vec3(0.0f, 1.0f, 0.0f));
+  auto wp                    = get_world_position();
+  glm::mat4 light_view       = glm::lookAt(wp, glm::vec3(0.0f), glm::vec3(0.0f, 1.0f, 0.0f));
   glm::mat4 light_projection = glm::ortho(-m, m, -m, m, near_plane, far_plane);
   return light_projection * light_view;
 }
@@ -475,9 +475,9 @@ void Renderer::render(Camera& camera, Object3D& scene)
   scene.update_world_matrix(false);
 
   RenderContext context;
-  context.camera = &camera;
-  context.shadow_map = shadow_map;
-  context.shadow_caster = nullptr;
+  context.camera           = &camera;
+  context.shadow_map       = shadow_map;
+  context.shadow_caster    = nullptr;
   context.background_color = background;
 
   scene.traverse([&context](Object3D* obj) {
@@ -546,7 +546,7 @@ void Mesh::draw_self(RenderContext& context)
 
     for (int i = 0; i < context.lights.size(); i++) {
       auto index = std::to_string(i);
-      auto type = context.lights[i]->type;
+      auto type  = context.lights[i]->type;
 
       shader->uniform("u_Lights[" + index + "].type", type);
       shader->uniform("u_Lights[" + index + "].color", context.lights[i]->rgb);
@@ -648,9 +648,9 @@ void FirstPersonController::move(const Direction& direction)
 void OrbitController::update(Object3D& target, const glm::vec3& center, float dt)
 {
   glm::vec3 front;
-  front.x = cos(glm::radians(m_yaw)) * cos(glm::radians(m_pitch));
-  front.y = sin(glm::radians(m_pitch));
-  front.z = sin(glm::radians(m_yaw)) * cos(glm::radians(m_pitch));
+  front.x     = cos(glm::radians(m_yaw)) * cos(glm::radians(m_pitch));
+  front.y     = sin(glm::radians(m_pitch));
+  front.z     = sin(glm::radians(m_yaw)) * cos(glm::radians(m_pitch));
   auto offset = glm::normalize(front) * radius;
 
   const auto pos = center + offset;
@@ -856,10 +856,10 @@ std::shared_ptr<Geometry> make_plane_geometry(int x_elements, int y_elements, fl
 
   for (int y = 0; y < y_elements; y++) {
     for (int x = 0; x < x_elements; x++) {
-      auto bottom_left = glm::vec3((x + 0) * width, 0.0f, (y + 0) * height);
+      auto bottom_left  = glm::vec3((x + 0) * width, 0.0f, (y + 0) * height);
       auto bottom_right = glm::vec3((x + 1) * width, 0.0f, (y + 0) * height);
-      auto top_left = glm::vec3((x + 0) * width, 0.0f, (y + 1) * height);
-      auto top_right = glm::vec3((x + 1) * width, 0.0f, (y + 1) * height);
+      auto top_left     = glm::vec3((x + 0) * width, 0.0f, (y + 1) * height);
+      auto top_right    = glm::vec3((x + 1) * width, 0.0f, (y + 1) * height);
 
       auto tex_coord = glm::vec2(static_cast<float>(x) / static_cast<float>(x_elements),
                                  static_cast<float>(y) / static_cast<float>(y_elements));
@@ -917,8 +917,8 @@ void Billboard::draw_self(RenderContext& context)
 
   auto camera = context.camera;
 
-  auto view = camera->get_view_matrix();
-  glm::vec3 up = {view[0][1], view[1][1], view[2][1]};
+  auto view       = camera->get_view_matrix();
+  glm::vec3 up    = {view[0][1], view[1][1], view[2][1]};
   glm::vec3 right = {view[0][0], view[1][0], view[2][0]};
 
   shader.bind();
@@ -975,7 +975,7 @@ void Skybox::draw_self(RenderContext& context)
 void SkyboxMaterial::bind()
 {
   gl::Shader* shader = get_shader();
-  int unit = 2;
+  int unit           = 2;
   cubemap->bind(unit);
   shader->bind();
   shader->uniform("u_Skybox", unit);
diff --git a/OpenGL_Flightsim/src/gfx.h b/OpenGL_Flightsim/src/gfx.h
index 962bd3a..d480426 100644
--- a/OpenGL_Flightsim/src/gfx.h
+++ b/OpenGL_Flightsim/src/gfx.h
@@ -108,9 +108,9 @@ struct ElementBufferObject {
 };
 
 struct TextureParams {
-  bool flip_vertically = false;
-  GLint texture_wrap_s = GL_REPEAT;
-  GLint texture_wrap_t = GL_REPEAT;
+  bool flip_vertically     = false;
+  GLint texture_wrap_s     = GL_REPEAT;
+  GLint texture_wrap_t     = GL_REPEAT;
   GLint texture_min_filter = GL_LINEAR_MIPMAP_LINEAR;
   GLint texture_mag_filter = GL_NEAREST;
 };
@@ -157,8 +157,8 @@ struct RenderContext {
 };
 
 #define OBJ3D_TRANSFORM 1U << 0U
-#define OBJ3D_ROTATE 1U << 1U
-#define OBJ3D_SCALE 1U << 2U
+#define OBJ3D_ROTATE    1U << 1U
+#define OBJ3D_SCALE     1U << 2U
 
 class Object3D
 {
@@ -179,7 +179,7 @@ class Object3D
   std::vector<Object3D*> children;
   glm::mat4 transform;
   bool receive_shadow = true;
-  bool visible = true;
+  bool visible        = true;
 
   Object3D& add(Object3D* child);
 
@@ -211,7 +211,7 @@ class Object3D
   void traverse(const std::function<bool(Object3D*)>& func);
 
  protected:
-  bool m_dirty_dof = false;
+  bool m_dirty_dof       = false;
   bool m_dirty_transform = false;
 
   glm::vec3 m_position;
@@ -241,7 +241,7 @@ class Light : public Object3D
 {
  public:
   enum LightType {
-    POINT = 0,
+    POINT       = 0,
     DIRECTIONAL = 1,
   };
 
@@ -420,9 +420,9 @@ class Renderer
     };
 
     auto geometry = std::make_shared<Geometry>(quad_vertices, Geometry::POS_UV);
-    auto texture = std::make_shared<gfx::gl::Texture>(shadow_map->depth_map.id);
+    auto texture  = std::make_shared<gfx::gl::Texture>(shadow_map->depth_map.id);
     auto material = std::make_shared<ScreenMaterial>(texture);
-    screen_quad = std::make_shared<Mesh>(geometry, material);
+    screen_quad   = std::make_shared<Mesh>(geometry, material);
   }
 
   ~Renderer()
@@ -481,7 +481,7 @@ class OrbitController
   void move_mouse(float x, float y);
 
  private:
-  float m_yaw = 0.0f;
+  float m_yaw   = 0.0f;
   float m_pitch = 0.0f;
 };
 };  // namespace gfx
diff --git a/OpenGL_Flightsim/src/main.cpp b/OpenGL_Flightsim/src/main.cpp
index d86cbff..2539292 100644
--- a/OpenGL_Flightsim/src/main.cpp
+++ b/OpenGL_Flightsim/src/main.cpp
@@ -14,7 +14,7 @@
 #include "../lib/imgui/imgui_impl_sdl2.h"
 #include "ai.h"
 #include "clipmap.h"
-#include "collisions.h"
+#include "collider.h"
 #include "flightmodel.h"
 #include "gfx.h"
 #include "phi.h"
@@ -36,16 +36,16 @@ EQ      control yaw
 JK      control thrust
 )";
 
-#define CLIPMAP 1
-#define SKYBOX 1
-#define SMOOTH_CAMERA 1
-#define NPC_AIRCRAFT 0
+#define CLIPMAP            1
+#define SKYBOX             1
+#define SMOOTH_CAMERA      1
+#define NPC_AIRCRAFT       0
 #define SHOW_MASS_ELEMENTS 0
-#define USE_PID 1
+#define USE_PID            1
 
 /* select flightmodel */
-#define FAST_JET 0
-#define CESSNA 1
+#define FAST_JET    0
+#define CESSNA      1
 #define FLIGHTMODEL FAST_JET
 
 #if 0
@@ -73,15 +73,14 @@ struct GameObject {
 
   void update(float dt)
   {
-    airplane.update_flightmodel(dt);
+    airplane.update(dt);
     transform.set_transform(airplane.position, airplane.orientation);
     collider.center = airplane.position;
   }
 };
 
 void get_keyboard_state(Joystick& joystick, phi::Seconds dt);
-void solve_constraints(phi::RigidBody& rigid_body);
-void apply_to_object3d(const phi::RigidBody& rigid_body, gfx::Object3D& object);
+void apply_to_object3d(const phi::RigidBody& rigid_body, gfx::Object3D* object);
 
 int main(void)
 {
@@ -96,11 +95,15 @@ int main(void)
                                         RESOLUTION.y, SDL_WINDOW_OPENGL);
 
   SDL_GLContext context = SDL_GL_CreateContext(window);
-  glewExperimental = GL_TRUE;
+  glewExperimental      = GL_TRUE;
 
   if (GLEW_OK != glewInit()) return -1;
 
   std::cout << glGetString(GL_VERSION) << std::endl;
+  std::cout << glGetString(GL_SHADING_LANGUAGE_VERSION) << std::endl;
+  std::cout << glGetString(GL_VENDOR) << std::endl;
+  std::cout << glGetString(GL_RENDERER) << std::endl;
+
   std::cout << USAGE << std::endl;
 
   IMGUI_CHECKVERSION();
@@ -125,25 +128,25 @@ int main(void)
   // use pid for keyboard control
   PID pitch_control_pid(1.0f, 0.0f, 0.0f);
 
-  int num_joysticks = SDL_NumJoysticks();
+  int num_joysticks     = SDL_NumJoysticks();
   bool joystick_control = num_joysticks > 0;
   if (joystick_control) {
     std::cout << "found " << num_joysticks << " joysticks\n";
     SDL_JoystickEventState(SDL_ENABLE);
     SDL_Joystick* sdl_joystick = SDL_JoystickOpen(0);
-    joystick.num_axis = SDL_JoystickNumAxes(sdl_joystick);
-    joystick.num_hats = SDL_JoystickNumHats(sdl_joystick);
-    joystick.num_buttons = SDL_JoystickNumButtons(sdl_joystick);
+    joystick.num_axis          = SDL_JoystickNumAxes(sdl_joystick);
+    joystick.num_hats          = SDL_JoystickNumHats(sdl_joystick);
+    joystick.num_buttons       = SDL_JoystickNumButtons(sdl_joystick);
     printf("found %d buttons, %d axis\n", joystick.num_buttons, joystick.num_axis);
   }
 
   gfx::Renderer renderer(RESOLUTION.x, RESOLUTION.y);
 
   gfx::gl::TextureParams params = {.flip_vertically = true};
-  auto tex = make_shared<gfx::gl::Texture>("assets/textures/f16_large.jpg", params);
-  auto texture = make_shared<gfx::Phong>(tex);
-  auto obj = gfx::load_obj("assets/models/falcon.obj");
-  auto model = std::make_shared<gfx::Geometry>(obj, gfx::Geometry::POS_NORM_UV);
+  auto tex                      = make_shared<gfx::gl::Texture>("assets/textures/f16_large.jpg", params);
+  auto texture                  = make_shared<gfx::Phong>(tex);
+  auto obj                      = gfx::load_obj("assets/models/falcon.obj");
+  auto model                    = std::make_shared<gfx::Geometry>(obj, gfx::Geometry::POS_NORM_UV);
 
   gfx::Object3D scene;
 
@@ -170,42 +173,48 @@ int main(void)
   Clipmap clipmap;
   scene.add(&clipmap);
 #endif
+#if 1
+  int width, height, channels;
+  const std::string heightmap_path = "assets/textures/terrain/1/heightmap.png";
+  uint8_t* data                    = gfx::gl::Texture::load_image(heightmap_path, &width, &height, &channels, 0);
+  col::Heightmap terrain_collider(data, width, height, channels);
+#endif
 
   std::vector<GameObject*> objects;
 
 #if (FLIGHTMODEL == CESSNA)
-  constexpr float speed = phi::units::meter_per_second(200.0f /* km/h */);
+  constexpr float speed    = phi::units::meter_per_second(200.0f /* km/h */);
   constexpr float altitude = 800.0f;
 
   // airplane mass
   const float mass = 1000.0f;
 
   // engine
-  const float rpm = 2400.0f;
-  const float horsepower = 160.0f;
+  const float rpm           = 2400.0f;
+  const float horsepower    = 160.0f;
   const float prop_diameter = 1.9f;
 
   // main wing
   const float total_wing_area = 16.17f;
   const float total_wing_span = 11.00f;
-  const float main_wing_span = total_wing_span / 2;
-  const float main_wing_area = total_wing_area / 2;
+  const float main_wing_span  = total_wing_span / 2;
+  const float main_wing_area  = total_wing_area / 2;
   const float main_wing_chord = main_wing_area / main_wing_span;
 
   // aileron
-  const float aileron_area = 1.70f;
-  const float aileron_span = 1.26f;
+  const float aileron_area  = 1.70f;
+  const float aileron_span  = 1.26f;
   const auto aileron_offset = glm::vec3(-main_wing_chord * 0.75f, 0.0f, 0.0f);
 
   // horizontal tail
   const float elevator_area = 1.35f;
-  const float h_tail_area = 2.0f + elevator_area;
-  const float h_tail_span = 2.0f;
-  const float h_tail_chord = h_tail_area / h_tail_span;
+  const float h_tail_area   = 2.0f + elevator_area;
+  const float h_tail_span   = 2.0f;
+  const float h_tail_chord  = h_tail_area / h_tail_span;
 
   // vertical tail
-  const float v_tail_area = 2.04f;  // modified
-  const float v_tail_span = 2.04f;
+  const float v_tail_area  = 2.04f;  // modified
+  const float v_tail_span  = 2.04f;
   const float v_tail_chord = v_tail_area / v_tail_span;
 
   const float wing_offset = -0.2f;
@@ -271,10 +280,17 @@ int main(void)
   auto engine = new PropellorEngine(horsepower, rpm, prop_diameter);
 
 #elif (FLIGHTMODEL == FAST_JET)
-  constexpr float altitude = 3000.0f;
+
+  glm::vec3 position = glm::vec3(-7000.0f, 0, 0.0f);
+
+  // constexpr float altitude = 700.0f;
+  float altitude = terrain_collider.get_height(glm::vec2{position.x, position.z}) + 100.0f;
+  position.y     = altitude;
+
   constexpr float speed = phi::units::meter_per_second(500.0f /* km/h */);
+  // constexpr float speed = 0.0f;
 
-  const float mass = 10000.0f;
+  const float mass   = 10000.0f;
   const float thrust = 50000.0f;
 
   const float wing_offset = -1.0f;
@@ -309,10 +325,10 @@ int main(void)
 #endif
 
   GameObject player = {.transform = gfx::Mesh(model, texture),
-                       .airplane = Airplane(mass, engine, inertia, wings),
-                       .collider = col::Sphere({0.0f, 0.0f, 0.0f}, 5.0f)};
+                       .airplane  = Airplane(mass, inertia, wings, engine),
+                       .collider  = col::Sphere({0.0f, 0.0f, 0.0f}, 15.0f)};
 
-  player.airplane.position = glm::vec3(-7000.0f, altitude, 0.0f);
+  player.airplane.position = position;
   player.airplane.velocity = glm::vec3(speed, 0.0f, 0.0f);
   scene.add(&player.transform);
   objects.push_back(&player);
@@ -321,7 +337,7 @@ int main(void)
   auto red_texture = make_shared<gfx::Phong>(glm::vec3(1.0f, 0.0f, 0.0f));
 
   for (int i = 0; i < mass_elements.size(); i++) {
-    auto& mass = mass_elements[i];
+    auto& mass   = mass_elements[i];
     auto element = new gfx::Mesh(gfx::make_cube_geometry(1.0f), red_texture);
     element->set_position(mass.offset);
     element->set_scale(mass.size);
@@ -331,10 +347,10 @@ int main(void)
 
 #if NPC_AIRCRAFT
   GameObject npc = {.transform = gfx::Mesh(model, texture),
-                    .airplane = Airplane(mass, engine, inertia, wings),
-                    .collider = col::Sphere({0.0f, 0.0f, 0.0f}, 5.0f)};
+                    .airplane  = Airplane(mass, inertia, wings, engine),
+                    .collider  = col::Sphere({0.0f, 0.0f, 0.0f}, 15.0f)};
 
-  npc.airplane.position = glm::vec3(-6990.0f, altitude, 10.0f);
+  npc.airplane.position = glm::vec3(-6950.0f, altitude, 10.0f);
   npc.airplane.velocity = glm::vec3(speed, 0.0f, 0.0f);
   scene.add(&npc.transform);
   objects.push_back(&npc);
@@ -348,7 +364,7 @@ int main(void)
 #endif
 
 #if 1
-  float size = 0.1f;
+  float size                = 0.1f;
   float projection_distance = 150.0f;
   glm::vec3 green(0.0f, 1.0f, 0.0f);
   gfx::Billboard cross(make_shared<gfx::gl::Texture>("assets/textures/sprites/cross.png"), green);
@@ -377,14 +393,6 @@ int main(void)
 
   gfx::OrbitController controller(30.0f);
 
-#if 1
-  int width, height, channels;
-  const std::string heightmap_path = "assets/textures/terrain/1/heightmap.png";
-  uint8_t* data = gfx::gl::Texture::load_image(heightmap_path, &width, &height, &channels, 0);
-  col::Heightmap terrain_collider(data, width, height, channels);
-
-#endif
-
   SDL_Event event;
   bool quit = false, paused = false, orbit = false;
   uint64_t last = 0, now = SDL_GetPerformanceCounter();
@@ -397,14 +405,14 @@ int main(void)
   while (!quit) {
     // delta time in seconds
     last = now;
-    now = SDL_GetPerformanceCounter();
-    dt = static_cast<phi::Seconds>((now - last) / static_cast<phi::Seconds>(SDL_GetPerformanceFrequency()));
+    now  = SDL_GetPerformanceCounter();
+    dt   = static_cast<phi::Seconds>((now - last) / static_cast<phi::Seconds>(SDL_GetPerformanceFrequency()));
     flight_time += dt;
     dt = std::min(dt, 0.02f);
 
     if ((timer += dt) >= 1.0f) {
       timer = 0.0f;
-      fps = 1.0f / dt;
+      fps   = 1.0f / dt;
     }
 
     while (SDL_PollEvent(&event) != 0) {
@@ -445,7 +453,7 @@ int main(void)
         }
         case SDL_JOYAXISMOTION: {
           if ((event.jaxis.value < -3200) || (event.jaxis.value > 3200)) {
-            uint8_t axis = event.jaxis.axis;
+            uint8_t axis  = event.jaxis.axis;
             int16_t value = event.jaxis.value;
             switch (axis) {
               case 0:
@@ -495,12 +503,12 @@ int main(void)
 
     if ((hud_timer += dt) > 0.1f) {
       hud_timer = 0.0f;
-      alt = player.airplane.get_altitude();
-      spd = phi::units::kilometer_per_hour(player.airplane.get_speed());
-      ias = phi::units::kilometer_per_hour(player.airplane.get_ias());
-      thr = static_cast<int>(player.airplane.engine->throttle * 100.0f);
-      gee = player.airplane.get_g();
-      aoa = player.airplane.get_aoa();
+      alt       = player.airplane.get_altitude();
+      spd       = phi::units::kilometer_per_hour(player.airplane.get_speed());
+      ias       = phi::units::kilometer_per_hour(player.airplane.get_ias());
+      thr       = player.airplane.engine->throttle;
+      gee       = player.airplane.get_g();
+      aoa       = player.airplane.get_aoa();
     }
 
     ImGui::SetNextWindowPos(ImVec2(10, 10));
@@ -510,16 +518,16 @@ int main(void)
     ImGui::Text("ALT:   %.1f m", alt);
     ImGui::Text("SPD:   %.1f km/h", spd);
     ImGui::Text("IAS:   %.1f km/h", ias);
-    ImGui::Text("THR:   %d %%", thr);
+    ImGui::Text("THR:   %.1f %%", player.airplane.engine->throttle * 100.0f);
     ImGui::Text("G:     %.1f", gee);
     ImGui::Text("AoA:   %.2f", aoa);
     ImGui::Text("Trim:  %.2f", player.airplane.joystick.w);
     ImGui::Text("FPS:   %.1f", fps);
     ImGui::End();
 
-#if 0
+#if 1
     auto angular_velocity = glm::degrees(player.airplane.angular_velocity);
-    auto orientation = glm::degrees(glm::eulerAngles(glm::normalize(player.airplane.orientation)));
+    auto orientation      = glm::degrees(player.airplane.get_euler_angles());
 
     ImVec2 size(140, 140);
     ImGui::SetNextWindowPos(ImVec2(RESOLUTION.x - size.y - 10.0f, RESOLUTION.y - size.y - 10.0f));
@@ -540,10 +548,10 @@ int main(void)
 
     player.airplane.joystick = glm::vec4(joystick.aileron, joystick.rudder, joystick.elevator, joystick.trim);
 #if USE_PID
-    if (!joystick_control) {
-      float max_av = 45.0f;  // deg/s
-      float target_av = max_av * joystick.elevator;
-      float current_av = glm::degrees(player.airplane.angular_velocity.z);
+    {
+      float max_av               = 45.0f;  // deg/s
+      float target_av            = max_av * joystick.elevator;
+      float current_av           = glm::degrees(player.airplane.angular_velocity.z);
       player.airplane.joystick.z = pitch_control_pid.calculate(current_av, target_av, dt);
     }
 #endif
@@ -551,7 +559,7 @@ int main(void)
 
 #if NPC_AIRCRAFT
     target_marker.visible = glm::length(camera.get_world_position() - npc.airplane.position) > 500.0f;
-    fly_towards(npc.airplane, player.airplane.position);
+    // fly_towards(npc.airplane, player.airplane.position);
 #endif
 
     if (!paused) {
@@ -563,17 +571,36 @@ int main(void)
       for (int i = 0; i < objects.size(); i++) {
         auto& a = objects[i];
 #if 1
-        if (col::test_collision(terrain_collider, a->airplane.position)) {
-          printf("[%.1f] terrain collision!\n", flight_time);
+        float terrain_height;
+        float gear_height = 2.0f;
+        if (col::test_collision(terrain_collider, a->airplane.position - glm::vec3(0.0f, gear_height, 0.0f),
+                                &terrain_height)) {
+          // printf("[%.1f] terrain collision!, %f\n", flight_time, terrain_height);
+
+          if (a->airplane.velocity.y < 0.0f) {
+            auto euler              = a->airplane.get_euler_angles();
+            euler.x                 = 0;                        // roll
+            euler.z                 = std::max(0.0f, euler.z);  // pitch, allow pitching up
+            a->airplane.orientation = glm::quat(euler);
+          }
+
+          a->airplane.is_landed  = true;
+          a->airplane.velocity.y = std::max(0.0f, a->airplane.velocity.y);
+          a->airplane.position.y = terrain_height + gear_height;
+
+        } else {
+          a->airplane.is_landed = false;
         }
 #endif
 
         for (int j = i + 1; j < objects.size(); j++) {
           auto& b = objects[j];
-          if (col::test_collision(a->collider, b->collider)) {
-            printf("[%.1f] collision!\n", flight_time);
-            auto collision_normal = glm::normalize(a->airplane.position - b->airplane.position);
-            phi::linear_collision_response(&a->airplane, &b->airplane, collision_normal, 0.5f);
+
+          phi::CollisionInfo collision{};
+
+          if (col::test_collision(a->collider, b->collider, &collision)) {
+            printf("[%.1f] collision!, p = %f\n", flight_time, collision.penetration);
+            phi::RigidBody::impulse_collision_response(&a->airplane, &b->airplane, collision);
           }
         }
       }
@@ -613,7 +640,7 @@ inline float center(float value, float factor, float dt)
 
 void get_keyboard_state(Joystick& joystick, phi::Seconds dt)
 {
-  const glm::vec3 factor = {3.0f, 0.5f, 1.0f};  // roll, yaw, pitch
+  const glm::vec3 factor    = {3.0f, 0.5f, 1.0f};  // roll, yaw, pitch
   const uint8_t* key_states = SDL_GetKeyboardState(NULL);
 
   if (key_states[SDL_SCANCODE_A] || key_states[SDL_SCANCODE_LEFT]) {
@@ -629,7 +656,7 @@ void get_keyboard_state(Joystick& joystick, phi::Seconds dt)
   } else if (key_states[SDL_SCANCODE_S] || key_states[SDL_SCANCODE_DOWN]) {
     joystick.elevator = move(joystick.elevator, -factor.z, dt);
   } else if (joystick.num_axis <= 0) {
-    joystick.elevator = center(joystick.elevator, factor.z, dt);
+    joystick.elevator = center(joystick.elevator, factor.z * 3.0f, dt);
   }
 
   if (key_states[SDL_SCANCODE_E]) {
@@ -656,14 +683,7 @@ void get_keyboard_state(Joystick& joystick, phi::Seconds dt)
   }
 }
 
-void apply_to_object3d(const phi::RigidBody& rigid_body, gfx::Object3D& object3d)
-{
-  object3d.set_transform(rigid_body.position, rigid_body.orientation);
-}
-
-void solve_constraints(phi::RigidBody& rigid_body)
+void apply_to_object3d(const phi::RigidBody& rigid_body, gfx::Object3D* object3d)
 {
-  if (rigid_body.position.y <= 0) {
-    rigid_body.position.y = 0, rigid_body.velocity.y = 0;
-  }
+  object3d->set_transform(rigid_body.position, rigid_body.orientation);
 }
diff --git a/OpenGL_Flightsim/src/phi.h b/OpenGL_Flightsim/src/phi.h
index cce3bcf..4a62cd9 100644
--- a/OpenGL_Flightsim/src/phi.h
+++ b/OpenGL_Flightsim/src/phi.h
@@ -1,3 +1,7 @@
+/*
+    copyright (c) 2023 jakob maier
+    'phi.h' is a simple, header-only rigidbody physics library.
+*/
 #pragma once
 
 #include <glm/gtc/matrix_transform.hpp>
@@ -5,31 +9,37 @@
 #include <glm/gtx/vector_angle.hpp>
 #include <glm/mat4x4.hpp>
 #include <glm/vec3.hpp>
+#include <vector>
+
+// RigidBody::update() can be marked virtual
+#if 0
+#define RB_VIRTUAL_UPDATE
+#else
+#define RB_VIRTUAL_UPDATE virtual
+#endif
 
 namespace phi
 {
 
 typedef float Seconds;
-typedef float Radians;
-typedef float Degrees;
 
 // constants
-constexpr float EPSILON = 1e-8f;
+constexpr float EPSILON       = 1e-8f;
 constexpr float EARTH_GRAVITY = 9.80665f;
-constexpr float PI = 3.141592653589793f;
+constexpr float PI            = 3.141592653589793f;
 
 // directions in body space
-constexpr glm::vec3 UP = {0.0f, 1.0f, 0.0f};
-constexpr glm::vec3 DOWN = {0.0f, -1.0f, 0.0f};
-constexpr glm::vec3 RIGHT = {0.0f, 0.0f, 1.0f};
-constexpr glm::vec3 LEFT = {0.0f, 0.0f, -1.0f};
-constexpr glm::vec3 FORWARD = {1.0f, 0.0f, 0.0f};
-constexpr glm::vec3 BACKWARD = {-1.0f, 0.0f, 0.0f};
-
 constexpr glm::vec3 X_AXIS = {1.0f, 0.0f, 0.0f};
 constexpr glm::vec3 Y_AXIS = {0.0f, 1.0f, 0.0f};
 constexpr glm::vec3 Z_AXIS = {0.0f, 0.0f, 1.0f};
 
+constexpr glm::vec3 FORWARD  = X_AXIS;
+constexpr glm::vec3 UP       = Y_AXIS;
+constexpr glm::vec3 RIGHT    = Z_AXIS;
+constexpr glm::vec3 BACKWARD = -X_AXIS;
+constexpr glm::vec3 DOWN     = -Y_AXIS;
+constexpr glm::vec3 LEFT     = -Z_AXIS;
+
 // utility functions
 template <typename T>
 constexpr inline T sq(T a)
@@ -68,6 +78,7 @@ inline T move_towards(T current, T target, T speed)
 }
 
 // inertia tensor calculations
+// formulas according to https://en.wikipedia.org/wiki/List_of_moments_of_inertia
 namespace inertia
 {
 struct Element {
@@ -80,14 +91,16 @@ struct Element {
   float volume() const { return size.x * size.y * size.z; }
 };
 
-constexpr glm::vec3 cube(const glm::vec3& size, float mass)
+// solid sphere
+constexpr glm::vec3 sphere(float radius, float mass) { return glm::vec3((2.0f / 5.0f) * mass * sq(radius)); }
+
+// cube with side length 'size'
+constexpr glm::vec3 cube(float size, float mass) { return glm::vec3((1.0f / 6.0f) * mass * sq(size)); }
+
+constexpr glm::vec3 cuboid(const glm::vec3& size, float mass)
 {
   const float C = (1.0f / 12.0f) * mass;
-  glm::vec3 I(0.0f);
-  I.x = C * (sq(size.y) + sq(size.z));
-  I.y = C * (sq(size.x) + sq(size.z));
-  I.z = C * (sq(size.x) + sq(size.y));
-  return I;
+  return glm::vec3(sq(size.y) + sq(size.z), sq(size.x) + sq(size.z), sq(size.x) + sq(size.y)) * C;
 }
 
 constexpr glm::vec3 cylinder(float radius, float length, float mass)
@@ -99,21 +112,26 @@ constexpr glm::vec3 cylinder(float radius, float length, float mass)
   return I;
 }
 
-// inertia tensor
-glm::mat3 tensor(const glm::vec3& moment_of_inertia)
+// helper function for the creation of a cuboid mass element
+constexpr Element cube(const glm::vec3& position, const glm::vec3& size, float mass)
 {
-  return {
-      moment_of_inertia.x, 0.0f, 0.0f, 0.0f, moment_of_inertia.y, 0.0f, 0.0f, 0.0f, moment_of_inertia.z,
-  };
+  return {size, position, cuboid(size, mass), position, mass};
 }
 
-constexpr Element cube(const glm::vec3& position, const glm::vec3& size, float mass = 1.0f)
+// inertia tensor from moment of inertia
+constexpr glm::mat3 tensor(const glm::vec3& moment_of_inertia)
 {
-  return {.size = size, .position = position, .inertia = cube(size, mass), .offset = position, .mass = mass};
+  // clang-format off
+  return {
+      moment_of_inertia.x, 0.0f, 0.0f, 
+      0.0f, moment_of_inertia.y, 0.0f, 
+      0.0f, 0.0f, moment_of_inertia.z,
+  };
+  // clang-format on
 }
 
 // distribute mass among elements depending on element volume
-void compute_uniform_mass(std::vector<Element>& elements, float mass)
+void compute_uniform_mass(std::vector<Element>& elements, float total_mass)
 {
   float total_volume = 0.0f;
   for (const auto& element : elements) {
@@ -121,7 +139,7 @@ void compute_uniform_mass(std::vector<Element>& elements, float mass)
   }
 
   for (auto& element : elements) {
-    element.mass = (element.volume() / total_volume) * mass;
+    element.mass = (element.volume() / total_volume) * total_mass;
   }
 }
 
@@ -162,8 +180,15 @@ glm::mat3 tensor(std::vector<Element>& elements, bool precomputed_offset = false
     *cg = center_of_gravity;
   }
 
-  return {Ixx, -Ixy, -Ixz, -Ixy, Iyy, -Iyz, -Ixz, -Iyz, Izz};
+  // clang-format off
+  return {
+      Ixx, -Ixy, -Ixz, 
+      -Ixy, Iyy, -Iyz, 
+      -Ixz, -Iyz, Izz
+  };
+  // clang-format on
 }
+
 };  // namespace inertia
 
 // unit conversions
@@ -180,14 +205,25 @@ constexpr inline float kelvin(float celsius) { return celsius - 273.15f; }
 constexpr inline float watts(float horsepower) { return horsepower * 745.7f; }
 };  // namespace units
 
+struct CollisionInfo {
+  glm::vec3 point;
+  glm::vec3 normal;
+  float penetration;
+};
+
+// default rigid body is a sphere with radius 1 meter and a mass of 100 kg
+constexpr float DEFAULT_RB_MASS            = 100.0f;
+constexpr glm::mat3 DEFAULT_RB_INERTIA     = inertia::tensor(inertia::sphere(1.0f, DEFAULT_RB_MASS));
+constexpr glm::quat DEFAULT_RB_ORIENTATION = glm::quat(1.0f, 0.0f, 0.0f, 0.0f);
+
 struct RigidBodyParams {
-  float mass = 1.0f;
-  glm::mat3 inertia{};
+  float mass        = DEFAULT_RB_MASS;
+  glm::mat3 inertia = DEFAULT_RB_INERTIA;
   glm::vec3 position{};
   glm::vec3 velocity{};
   glm::vec3 angular_velocity{};
-  glm::quat orientation = glm::quat(glm::vec3(0.0f));
-  bool apply_gravity = true;
+  glm::quat orientation = DEFAULT_RB_ORIENTATION;
+  bool apply_gravity    = true;
 };
 
 class RigidBody
@@ -201,13 +237,12 @@ class RigidBody
   glm::vec3 position{};                    // position in world space
   glm::quat orientation{};                 // orientation in world space
   glm::vec3 velocity{};                    // velocity in world space
-  glm::vec3 angular_velocity{};            // angular velocity in object space, x
-                                           // represents rotation around x axis
+  glm::vec3 angular_velocity{};            // angular velocity in object space, (x = roll, y = yaw, z = pitch)
   glm::mat3 inertia{}, inverse_inertia{};  // inertia tensor
   bool apply_gravity = true;
-  bool active = true;
+  bool active        = true;
 
-  RigidBody() : RigidBody({.mass = 1.0f, .inertia = inertia::tensor(inertia::cube(glm::vec3(1.0f), 1.0f))}) {}
+  RigidBody() : RigidBody({DEFAULT_RB_MASS, DEFAULT_RB_INERTIA}) {}
 
   RigidBody(const RigidBodyParams& params)
       : mass(params.mass),
@@ -221,15 +256,24 @@ class RigidBody
   {
   }
 
-  // get velocity of point in body space
+  // get angular velocity of relative point in body space
+  inline glm::vec3 get_point_angular_velocity(const glm::vec3& point) const
+  {
+    return glm::cross(angular_velocity, point);
+  }
+
+  // get velocity of relative point in body space
   inline glm::vec3 get_point_velocity(const glm::vec3& point) const
   {
-    return inverse_transform_direction(velocity) + glm::cross(angular_velocity, point);
+    return inverse_transform_direction(velocity) + get_point_angular_velocity(point);
   }
 
   // get velocity in body space
   inline glm::vec3 get_body_velocity() const { return inverse_transform_direction(velocity); }
 
+  // self explantory
+  inline float get_inverse_mass() const { return 1.0f / mass; };
+
   // force and point vectors are in body space
   inline void add_force_at_point(const glm::vec3& force, const glm::vec3& point)
   {
@@ -249,11 +293,23 @@ class RigidBody
   // set inertia tensor
   inline void set_inertia(const glm::mat3& tensor) { inertia = tensor, inverse_inertia = glm::inverse(tensor); }
 
+  // set moment of inertia
+  inline void set_inertia(const glm::vec3& moment) { set_inertia(inertia::tensor(moment)); }
+
   // linear impulse in world space
-  inline void add_impulse(const glm::vec3& impulse) { velocity += impulse / mass; }
+  inline void add_linear_impulse(const glm::vec3& impulse) { velocity += impulse / mass; }
 
   // linear impulse in body space
-  inline void add_relative_impulse(const glm::vec3& impulse) { velocity += transform_direction(impulse) / mass; }
+  inline void add_relative_linear_impulse(const glm::vec3& impulse) { velocity += transform_direction(impulse) / mass; }
+
+  // angular impulse in world space
+  inline void add_angular_impulse(const glm::vec3& impulse)
+  {
+    angular_velocity += inverse_transform_direction(impulse) * inverse_inertia;
+  }
+
+  // angular impulse in body space
+  inline void add_relative_angular_impulse(const glm::vec3& impulse) { angular_velocity += impulse * inverse_inertia; }
 
   // force vector in world space
   inline void add_force(const glm::vec3& force) { m_force += force; }
@@ -267,25 +323,39 @@ class RigidBody
   // torque vector in body space
   inline void add_relative_torque(const glm::vec3& torque) { m_torque += torque; }
 
+  // add terrain friction
+  glm::vec3 add_friction(const glm::vec3& normal, const glm::vec3& sliding_direction, float friction_coeff)
+  {
+    // https://en.wikipedia.org/wiki/Normal_force
+    // https://en.wikipedia.org/wiki/Friction
+    float weight      = mass * EARTH_GRAVITY;
+    auto normal_force = weight * std::max(glm::dot(normal, UP), 0.0f);
+    return -sliding_direction * normal_force * friction_coeff;
+  }
+
   // get speed
   inline float get_speed() const { return glm::length(velocity); }
 
+  // get euler angles in radians
+  inline glm::vec3 get_euler_angles() const { return glm::eulerAngles(orientation); }
+
   // get torque in body space
   inline glm::vec3 get_torque() const { return m_torque; }
 
   // get torque in world space
   inline glm::vec3 get_force() const { return m_force; }
 
-  // get forward direction in world space
+  // get rigidbody forward direction in world space
   inline glm::vec3 forward() const { return transform_direction(phi::FORWARD); }
 
-  // get up direction in world space
+  // get rigidbody up direction in world space
   inline glm::vec3 up() const { return transform_direction(phi::UP); }
 
-  // get right direction in world space
+  // get rigidbody right direction in world space
   inline glm::vec3 right() const { return transform_direction(phi::RIGHT); }
 
-  void update(Seconds dt)
+  // integrate RigidBody
+  RB_VIRTUAL_UPDATE void update(Seconds dt)
   {
     if (!active) return;
 
@@ -303,26 +373,86 @@ class RigidBody
     // reset accumulators
     m_force = glm::vec3(0.0f), m_torque = glm::vec3(0.0f);
   }
-};
 
-struct ForceGenerator {
-  virtual void apply_forces(phi::RigidBody* rigid_body, phi::Seconds dt) = 0;
-};
+  // restitution_coeff:  0 = perfectly inelastic, 1 = perfectly elastic
+  // impulse collision response without angular effects
+  static void linear_impulse_collision_response(RigidBody* a, RigidBody* b, const CollisionInfo& collision,
+                                                float restitution_coeff = 0.66f)
+  {
+    float total_inverse_mass = a->get_inverse_mass() + b->get_inverse_mass();
 
-// no angular effects, collision_normal goes from a to b
-// restitution_coeff:  0 = perfectly inelastic, 1 = perfectly elastic
-void linear_collision_response(phi::RigidBody* a, phi::RigidBody* b, const glm::vec3& collision_normal,
-                               float restitution_coeff = 0.5f)
-{
-  assert(0.0f <= restitution_coeff && restitution_coeff <= 1.0f);
+    // move objects so they are no longer colliding. heavier object gets moved less
+    a->position -= collision.normal * collision.penetration * (a->get_inverse_mass() / total_inverse_mass);
+    b->position += collision.normal * collision.penetration * (b->get_inverse_mass() / total_inverse_mass);
 
-  auto relative_velocity = a->velocity - b->velocity;
+    auto relative_velocity = b->velocity - a->velocity;
 
-  float impulse = (-(1 + restitution_coeff) * glm::dot(relative_velocity, collision_normal)) /
-                  (glm::dot(collision_normal, collision_normal) * (1 / a->mass + 1 / b->mass));
+    // force is highest in head on collision
+    float impulse_force = glm::dot(relative_velocity, collision.normal);
 
-  a->add_impulse(impulse * +collision_normal);
-  b->add_impulse(impulse * -collision_normal);
-}
+    // magnitude of impulse
+    float j = (-(1 + restitution_coeff) * impulse_force) / (total_inverse_mass);
+
+    auto impulse = j * collision.normal;
+
+    // apply linear impulse
+    a->add_linear_impulse(-impulse);
+    b->add_linear_impulse(+impulse);
+  }
+
+  // impulse collision response
+  static void impulse_collision_response(RigidBody* a, RigidBody* b, const CollisionInfo& collision,
+                                         float restitution_coeff = 0.66f)
+  {
+    float total_inverse_mass = a->get_inverse_mass() + b->get_inverse_mass();
+
+    // move objects so they are no longer colliding. heavier object gets moved less
+    a->position -= collision.normal * collision.penetration * (a->get_inverse_mass() / total_inverse_mass);
+    b->position += collision.normal * collision.penetration * (b->get_inverse_mass() / total_inverse_mass);
+
+    // location of collision point relative to rigidbody
+    auto a_relative = collision.point - a->position;
+    auto b_relative = collision.point - b->position;
+
+    // get velocity at this point (body space)
+    auto a_velocity = a->get_point_velocity(a_relative);
+    auto b_velocity = b->get_point_velocity(b_relative);
+
+    // relative velocity in world space
+    auto relative_velocity = b->transform_direction(b_velocity) - a->transform_direction(a_velocity);
+
+    // force is highest in a head on collision
+    float impulse_force = glm::dot(relative_velocity, collision.normal);
+
+    auto a_inertia         = glm::cross(a->inertia * glm::cross(a_relative, collision.normal), a_relative);
+    auto b_inertia         = glm::cross(b->inertia * glm::cross(b_relative, collision.normal), b_relative);
+    float angular_effect_1 = glm::dot(a_inertia + b_inertia, collision.normal);
+
+    float angular_effect_2 =
+        glm::dot(collision.normal, glm::cross((glm::cross(a_relative, collision.normal) / a->inertia), a_relative)) +
+        glm::dot(collision.normal, glm::cross((glm::cross(b_relative, collision.normal) / b->inertia), b_relative));
+
+    // TODO: find correct implementation
+    printf("a_1 = %f, a_2 = %f\n", angular_effect_1, angular_effect_2);
+    assert(std::abs(angular_effect_1 - angular_effect_2) < phi::EPSILON);
+
+    // magnitude of impulse
+    float j = (-(1 + restitution_coeff) * impulse_force) / (total_inverse_mass + angular_effect_1);
+
+    auto impulse = j * collision.normal;
+
+    // apply linear impulse
+    a->add_linear_impulse(-impulse);
+    b->add_linear_impulse(+impulse);
+
+    // apply angular impulse at position
+    a->add_angular_impulse(glm::cross(a_relative, -impulse));
+    b->add_angular_impulse(glm::cross(b_relative, +impulse));
+  }
+};
+
+struct ForceGenerator {
+  virtual void apply_forces(phi::RigidBody* rigid_body, phi::Seconds dt) = 0;
+};
 
 };  // namespace phi
diff --git a/OpenGL_Flightsim/src/pid.h b/OpenGL_Flightsim/src/pid.h
index fc78b31..81621fa 100644
--- a/OpenGL_Flightsim/src/pid.h
+++ b/OpenGL_Flightsim/src/pid.h
@@ -6,8 +6,8 @@
 
 class PID
 {
-  float integral = 0.0f;
-  bool initialized = false;
+  float integral       = 0.0f;
+  bool initialized     = false;
   float previous_value = 0.0f, previous_error = 0.0f;
 
   const bool use_value;
@@ -25,22 +25,22 @@ class PID
   float calculate(float current_value, float target_value, float dt)
   {
     float error = target_value - current_value;
-    float P = error * proportional_gain;
+    float P     = error * proportional_gain;
 
     if (!initialized) {
       previous_error = error;
       previous_value = current_value;
-      initialized = true;
+      initialized    = true;
     }
 
     integral += error * dt;
     float I = integral * integral_gain;
 
     float error_rate_of_change = (error - previous_error) / dt;
-    previous_error = error;
+    previous_error             = error;
 
     float value_rate_of_change = (current_value - previous_value) / dt;
-    previous_value = current_value;
+    previous_value             = current_value;
 
     float D = (use_value ? value_rate_of_change : error_rate_of_change) * derivative_gain;