local z=.001 local a=.001 local h=.0001 local n={}n.__index=n function n.new(z,a)assert(z,"\077\105\115\115\105\110\103 \097\114\103\117\109\101\110\116 \035\049\058 \116\097\114\103\101\116\086\097\108\117\101")a=a or{}return setmetatable({_targetValue=z;_frequency=a.frequency or 4;_dampingRatio=a.dampingRatio or 1},n)end function n.step(f,n,P)local M=f._dampingRatio local B=(f._frequency*2)*math.pi local s=f._targetValue local U=n.value local S=n.velocity or 0 local F=U-s local k=math.exp((-M*B)*P)local r,i if M==1 then r=(F*(1+B*P)+S*P)*k+s i=(S*(1-B*P)-F*((B*B)*P))*k elseif M<1 then local z=math.sqrt(1-M*M)local a=math.cos((B*z)*P)local n=math.sin((B*z)*P)local f if z>h then f=n/z else local a=P*B f=a+(((((a*a)*(z*z))*(z*z))/20-z*z)*((a*a)*a))/6 end local U if B*z>h then U=n/(B*z)else local a=B*z U=P+(((((P*P)*(a*a))*(a*a))/20-a*a)*((P*P)*P))/6 end r=(F*(a+M*f)+S*U)*k+s i=(S*(a-f*M)-F*(f*B))*k else local z=math.sqrt(M*M-1)local a=-B*(M-z)local h=-B*(M+z)local n=(S-F*a)/((2*B)*z)local f=F-n local U=f*math.exp(a*P)local k=n*math.exp(h*P)r=(U+k)+s i=U*a+k*h end local c=math.abs(i)<z and math.abs(r-s)<a return{complete=c;value=c and s or r;velocity=i}end return n
local VELOCITY_THRESHOLD = 0.001
local POSITION_THRESHOLD = 0.001

local EPS = 0.0001

local Spring = {}
Spring.__index = Spring

function Spring.new(targetValue, options)
	assert(targetValue, "Missing argument #1: targetValue")
	options = options or {}

	return setmetatable({
		_targetValue = targetValue,
		_frequency = options.frequency or 4,
		_dampingRatio = options.dampingRatio or 1,
	}, Spring)
end

function Spring:step(state, dt)
	-- Copyright 2018 Parker Stebbins (parker@fractality.io)
	-- github.com/Fraktality/Spring
	-- Distributed under the MIT license

	local d = self._dampingRatio
	local f = self._frequency * 2 * math.pi
	local g = self._targetValue
	local p0 = state.value
	local v0 = state.velocity or 0

	local offset = p0 - g
	local decay = math.exp(-d * f * dt)

	local p1, v1

	if d == 1 then -- Critically damped
		p1 = (offset * (1 + f * dt) + v0 * dt) * decay + g
		v1 = (v0 * (1 - f * dt) - offset * (f * f * dt)) * decay
	elseif d < 1 then -- Underdamped
		local c = math.sqrt(1 - d * d)

		local i = math.cos(f * c * dt)
		local j = math.sin(f * c * dt)

		-- Damping ratios approaching 1 can cause division by small numbers.
		-- To fix that, group terms around z=j/c and find an approximation for z.
		-- Start with the definition of z:
		--    z = sin(dt*f*c)/c
		-- Substitute a=dt*f:
		--    z = sin(a*c)/c
		-- Take the Maclaurin expansion of z with respect to c:
		--    z = a - (a^3*c^2)/6 + (a^5*c^4)/120 + O(c^6)
		--    z ≈ a - (a^3*c^2)/6 + (a^5*c^4)/120
		-- Rewrite in Horner form:
		--    z ≈ a + ((a*a)*(c*c)*(c*c)/20 - c*c)*(a*a*a)/6

		local z
		if c > EPS then
			z = j / c
		else
			local a = dt * f
			z = a + ((a * a) * (c * c) * (c * c) / 20 - c * c) * (a * a * a) / 6
		end

		-- Frequencies approaching 0 present a similar problem.
		-- We want an approximation for y as f approaches 0, where:
		--    y = sin(dt*f*c)/(f*c)
		-- Substitute b=dt*c:
		--    y = sin(b*c)/b
		-- Now reapply the process from z.

		local y
		if f * c > EPS then
			y = j / (f * c)
		else
			local b = f * c
			y = dt + ((dt * dt) * (b * b) * (b * b) / 20 - b * b) * (dt * dt * dt) / 6
		end

		p1 = (offset * (i + d * z) + v0 * y) * decay + g
		v1 = (v0 * (i - z * d) - offset * (z * f)) * decay
	else -- Overdamped
		local c = math.sqrt(d * d - 1)

		local r1 = -f * (d - c)
		local r2 = -f * (d + c)

		local co2 = (v0 - offset * r1) / (2 * f * c)
		local co1 = offset - co2

		local e1 = co1 * math.exp(r1 * dt)
		local e2 = co2 * math.exp(r2 * dt)

		p1 = e1 + e2 + g
		v1 = e1 * r1 + e2 * r2
	end

	local complete = math.abs(v1) < VELOCITY_THRESHOLD and math.abs(p1 - g) < POSITION_THRESHOLD

	return {
		complete = complete,
		value = complete and g or p1,
		velocity = v1,
	}
end

return Spring
