Pacejka 96 for tyre friction, combined friction using the Beckman method

.gitignore

@@ -16,3 +16,4 @@ linux
 atlas.png
 pdbs/
 game_web/
+.venv

game/game.odin

@@ -393,9 +393,9 @@ update_runtime_world :: proc(runtime_world: ^Runtime_World, dt: f32) {
 		front_wheels := turn_wheels
 		back_wheels := drive_wheels
 
-		DRIVE_IMPULSE :: 2000
+		DRIVE_IMPULSE :: 4000
 		BRAKE_IMPULSE :: 500
-		TURN_ANGLE :: -f32(20) * math.RAD_PER_DEG
+		TURN_ANGLE :: -f32(30) * math.RAD_PER_DEG
 		// 68% front, 32% rear
 		BRAKE_BIAS :: f32(0.68)
 

game/physics/pacejka.odin

@@ -0,0 +1,22 @@
+package physics
+
+import "core:math"
+
+Pacejka96_Params :: [11]f32
+
+DEFAULT_PACEJKA96_PARAMS :: Pacejka96_Params{1.6, 0, 1688, 0, 229, 0, 0, 0, -10, 0, 0}
+
+// X is slip ratio percentage [-100, 100] or slip angle in degrees, where positive angle is turning left
+// Output is the friction coefficient
+pacejka_96 :: proc(b: Pacejka96_Params, x: f32, f_z: f32, s_v: f32 = 0) -> f32 {
+	C := b[0]
+	D := (b[1] * f_z + b[2]) * f_z
+	B := ((b[3] * f_z * f_z + b[4] * f_z) * math.exp(-b[5] * f_z)) / (C * D)
+	E := -(1 - (b[6] * f_z * f_z + b[7] * f_z + b[8]))
+	s_h := b[9] * f_z + b[10]
+
+	X := x + s_h
+	y := D * math.sin(C * math.atan(B * X - E * (B * X - math.atan(B * X))))
+	Y := y + s_v
+	return Y / (f_z * 1000)
+}

game/physics/simulation.odin

@@ -4,6 +4,7 @@ import "bvh"
 import "collision"
 import "core:container/bit_array"
 import "core:fmt"
+import "core:log"
 import "core:math"
 import lg "core:math/linalg"
 import "core:math/rand"
@@ -549,24 +550,67 @@ pgs_solve_suspension :: proc(sim_state: ^Sim_State, config: Solver_Config, dt: f
 						apply_velocity_correction(body, incremental_impulse * dir, wheel_world_pos)
 					}
 
+					body_right := body_local_to_world_vec(body, Vec3{1, 0, 0})
 					right := wheel_get_right_vec(body, v)
 
-					contact_patch_linear_vel :=
-						body_vel_at_contact_patch + (v.radius * v.w * forward)
+					// Positive means spinning forward
+					wheel_spin_vel := -v.radius * v.w
+					ground_vel := lg.dot(body_vel_at_contact_patch, forward)
+					// contact_patch_linear_vel :=
+					// 	body_vel_at_contact_patch + (v.radius * v.w * forward)
+
+					slip_ratio :=
+						ground_vel == 0 ? 0 : clamp(wheel_spin_vel / ground_vel - 1, -1, 1)
+					slip_angle :=
+						lg.angle_between(forward, body_vel_at_contact_patch) * math.DEG_PER_RAD
+
+					OPTIMAL_SLIP_RATIO :: f32(0.075)
+					OPTIMAL_SLIP_ANGLE :: f32(8)
+					MAX_SLIP_LEN :: f32(2.0)
+
+					slip_vec := Vec2 {
+						slip_angle / OPTIMAL_SLIP_ANGLE / MAX_SLIP_LEN,
+						slip_ratio / OPTIMAL_SLIP_RATIO / MAX_SLIP_LEN,
+					}
+
+					slip_len := lg.length(slip_vec)
+					slip_len = slip_len == 0 ? 0 : min(slip_len, 1) / slip_len
+					slip_vec *= slip_len
+
+					log.debugf("slip_vec: %v", slip_vec)
+
+					long_friction :=
+						abs(
+							pacejka_96(
+								DEFAULT_PACEJKA96_PARAMS,
+								slip_ratio * 100,
+								max(abs(v.spring_impulse), 0.001) * inv_dt * 0.001,
+							),
+						) *
+						abs(slip_vec.y)
+					lat_friction :=
+						abs(
+							pacejka_96(
+								DEFAULT_PACEJKA96_PARAMS,
+								slip_angle,
+								max(abs(v.spring_impulse), 0.001) * inv_dt * 0.001,
+							),
+						) *
+						abs(slip_vec.x)
 
 					// Longitudinal friction
 					if true {
-						vel_long := lg.dot(contact_patch_linear_vel, forward)
+						// Wheel linear velocity relative to ground
+						relative_vel := ground_vel - wheel_spin_vel
 
-						friction := f32(1)
-						friction_clamp := abs(v.spring_impulse) * friction
+						friction_clamp := abs(v.spring_impulse) * long_friction
 
 						w_body := get_body_inverse_mass(body, forward, v.hit_point)
 
 						w_long := w_body + v.inv_inertia
 						inv_w_long := 1.0 / w_long
 
-						incremental_impulse := -inv_w_long * vel_long
+						incremental_impulse := -inv_w_long * relative_vel
 						new_total_impulse := clamp(
 							v.longitudinal_impulse + incremental_impulse,
 							-friction_clamp,
@@ -583,9 +627,7 @@ pgs_solve_suspension :: proc(sim_state: ^Sim_State, config: Solver_Config, dt: f
 						vel_contact := body_vel_at_contact_patch
 
 						lateral_vel := lg.dot(right, vel_contact)
-
-						friction := f32(1)
-						friction_clamp := -v.spring_impulse * friction
+						friction_clamp := -v.spring_impulse * lat_friction
 
 						incremental_impulse := -inv_w_normal * lateral_vel
 						new_total_impulse := clamp(
@@ -596,7 +638,7 @@ pgs_solve_suspension :: proc(sim_state: ^Sim_State, config: Solver_Config, dt: f
 						applied_impulse := new_total_impulse - v.lateral_impulse
 						v.lateral_impulse = new_total_impulse
 
-						apply_velocity_correction(body, applied_impulse * right, v.hit_point)
+						apply_velocity_correction(body, applied_impulse * body_right, v.hit_point)
 					}
 				} else {
 					v.lateral_impulse = 0
@@ -651,7 +693,7 @@ pgs_substep :: proc(sim_state: ^Sim_State, config: Solver_Config, dt: f32, inv_d
 				p := body_local_to_world(body, s.rel_pos)
 				hit_p := body_local_to_world(body, s.rel_pos + s.rel_dir * s.hit_t)
 				forward := wheel_get_forward_vec(body, s)
-				right := wheel_get_right_vec(body, s)
+				right := body_local_to_world_vec(body, Vec3{1, 0, 0})
 
 				apply_velocity_correction(
 					body,