-- Boids
-- Kartik Agaram, 2025-04-26
-- Requires LÖVE; tested on v11.5

Boids = {}
Personal_space = 12
Perceptual_distance = 70

function love.load()
	Window_side = 800
	love.window.setMode(Window_side, Window_side)
	love.graphics.setBackgroundColor(0.8,0.8,0.8)
	for _ = 1,80 do
		table.insert(Boids, make_boid(math.random(Window_side), math.random(Window_side)))
	end
end

function make_boid(x, y)
	return {pos={x=x, y=y}, velocity={x=0, y=0}}
end

function love.draw()
	love.graphics.setColor(0,0,0)
	draw_boids(Boids)
	draw_hud_stats()
end

function draw_boids(boids)
	local boid_width, boid_length = Personal_space/2, Personal_space
	for _, boid in ipairs(boids) do
		local pos, velocity = boid.pos, boid.velocity
		local dir = vnorm(velocity) or {x=0, y=-1}
		local p1 = vadd(pos, vscale(dir, boid_length/2))
		local p2 = vadd(
				pos,
				vscale(dir, -boid_length/2),
				vscale(vperpendicular_anticlockwise(dir),
						boid_width/2))
		local p3 = vadd(
				pos,
				vscale(dir, -boid_length/2),
				vscale(vperpendicular_clockwise(dir),
						boid_width/2))
		love.graphics.polygon('fill', p1.x, p1.y, p2.x, p2.y, p3.x, p3.y)
	end
end

function draw_hud_stats()
	local count = 0
	for _, boid in ipairs(Boids) do
		local p = boid.pos
		if p.x >= 0 and p.x <= Window_side
			and p.y >= 0 and p.y <= Window_side
		then
			count = count+1
		end
	end
	love.graphics.setColor(1,1,1)
	love.graphics.print(('visible: %d'):format(count), 100, 100)
end

function love.update(dt)
	-- velocities and accelerations are per frame
	local max_speed = 0.5
	local max_accel = 20
	local nan = tostring(0/0)
	for _,boid in ipairs(Boids) do
		boid.pos = vadd(boid.pos, boid.velocity)
		assert(tostring(boid.pos.x) ~= nan)
		assert(tostring(boid.pos.y) ~= nan)
	end
	for i,boid in ipairs(Boids) do
		local accel = {x=0, y=0}
		-- each rule returns a normalized vector
		-- all weighting occurs here
		accel = vadd(accel,
			vscale(avoid_others(boid, Boids), 20*dt))
		accel = vadd(accel,
			vscale(seek_others(boid, Boids), 10*dt))
		accel = vadd(accel,
			vscale(align_with_others(boid, Boids), 10*dt))
		accel = vadd(accel,
			vscale(remain_within_viewport(boid), 40*dt))
		local curr_heading = vnorm(boid.velocity)  -- could be nil
		accel = vclamp2(accel, max_accel*dt, curr_heading, math.pi/6*dt)
		boid.velocity = vadd(boid.velocity, accel)
		boid.velocity = vclamp(boid.velocity, max_speed)
	end
end

function avoid_others(boid, boids)
	local result = {x=0, y=0}
	local dr = 0
	for i, boid2 in ipairs(boids) do
		if boid2 ~= boid then
			local delta = vsub(boid.pos, boid2.pos)
			local dist = vlength(delta)
			if dist < Personal_space then
				result = vadd(result, vnorm(delta))
				dr = dr + 1
			end
		end
	end
	if dr == 0 then
		return result
	end
	return vnorm(result)
end

function seek_others(boid, boids)
	-- fly towards the center of neighbors
	local center = {x=0, y=0}
	local dr = 0
	for _, boid2 in ipairs(boids) do
		if boid2 ~= boid and vdist(boid.pos, boid2.pos) < Perceptual_distance then
			center = vadd(center, boid2.pos)
			dr = dr + 1
		end
	end
	if dr == 0 then
		return {x=0, y=0}
	end
	center = vscale(center, 1/dr)
	return vnorm(vsub(center, boid.pos)) or {x=0, y=0}
end

function align_with_others(boid, boids)
	-- aim in the average heading of neighbors
	local result = {x=0, y=0}
	for _, boid2 in ipairs(boids) do
		if boid2 ~= boid and vdist(boid.pos, boid2.pos) < Perceptual_distance then
			result = vadd(result, boid2.velocity)
		end
	end
	return vnorm(result) or {x=0, y=0}
end

function remain_within_viewport(boid)
	local result = {x=0, y=0}
	local margin = 50
	local x, y = boid.pos.x, boid.pos.y
	if x < margin then
		result.x = result.x + (margin-x)/margin
	end
	if x > Window_side - margin then
		result.x = result.x + (Window_side-margin-x)/margin
	end
	if y < margin then
		result.y = result.y + (margin-y)/margin
	end
	if y > Window_side - margin then
		result.y = result.y + (Window_side-margin-y)/margin
	end
	return vnorm(result) or result
end

-- make sure magnitude of v is max or less
function vclamp(v, max)
	local l = vlength(v)
	if l <= max then return v end
	return vscale(v, max/l)
end

-- make sure:
--  * magnitude of v is mag or less, and
--  * angle of v is no more than `ang` away from vector `heading`
function vclamp2(v, mag, heading, ang)
	assert(ang >= 0)
	if heading == nil then  -- we have no basis for preferring a direction
		return v
	end
	local r = vlength(v)
	if r == 0 then return v end  -- we have no direction
	local theta = math.atan2(v.y, v.x)
	local base = math.atan2(heading.y, heading.x)
	local dtheta = canonicalize_angle(theta - base)
	dtheta = math.max(-ang, math.min(dtheta, ang))
	theta = base + dtheta
	return {x=r*math.cos(theta), y=r*math.sin(theta)}
end

function canonicalize_angle(theta)
	-- fix angle between -pi and +pi
	if theta > math.pi then
		return - (2*math.pi - theta)  -- must be negative
	end
	if theta <= -math.pi then
		return 2*math.pi + theta  -- must be positive
	end
	return theta
end

-- length of a vector
function vlength(v)
	return (v.x^2 + v.y^2)^0.5
end

-- scalar distance between two vectors
function vdist(v1, v2)
	return vlength(vsub(v1, v2))
end

-- add multiple vectors together
function vadd(...)
	local args = {...}
	local x, y = 0, 0
	for _, v in ipairs(args) do
		x = x+v.x
		y = y+v.y
	end
	return {x=x, y=y}
end

-- subtract v2 from v1
function vsub(v1, v2)
	return {x=v1.x-v2.x, y=v1.y-v2.y}
end

-- scale vector to magnitude 1
-- can return nil
function vnorm(v)
	local l = vlength(v)
	if l == 0 then return end
	return vscale(v, 1/l)
end

-- scale vector by factor of k
function vscale(v, k)
	return {x=k*v.x, y=k*v.y}
end

-- turn v 90° counter-clockwise
function vperpendicular_anticlockwise(v)
	return {x=-v.y, y=v.x}
end

-- turn v 90° clockwise
function vperpendicular_clockwise(v)
	return {x=v.y, y=-v.x}
end

-- print a vector to terminal
function vprint(v)
	print('vector:', v.x, v.y)
end