local Q=.001 local R=.001 local w=.0001 local S={}S.__index=S function S.new(Q,R)assert(Q,"\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")R=R or{}return setmetatable({_targetValue=Q;_frequency=R.frequency or 4;_dampingRatio=R.dampingRatio or 1},S)end function S.step(l,S,u)local h=l._dampingRatio local o=(l._frequency*2)*math.pi local V=l._targetValue local Z=S.value local s=S.velocity or 0 local B=Z-V local y=math.exp((-h*o)*u)local W,L if h==1 then W=(B*(1+o*u)+s*u)*y+V L=(s*(1-o*u)-B*((o*o)*u))*y elseif h<1 then local Q=math.sqrt(1-h*h)local R=math.cos((o*Q)*u)local S=math.sin((o*Q)*u)local l if Q>w then l=S/Q else local R=u*o l=R+(((((R*R)*(Q*Q))*(Q*Q))/20-Q*Q)*((R*R)*R))/6 end local Z if o*Q>w then Z=S/(o*Q)else local R=o*Q Z=u+(((((u*u)*(R*R))*(R*R))/20-R*R)*((u*u)*u))/6 end W=(B*(R+h*l)+s*Z)*y+V L=(s*(R-l*h)-B*(l*o))*y else local Q=math.sqrt(h*h-1)local R=-o*(h-Q)local w=-o*(h+Q)local S=(s-B*R)/((2*o)*Q)local l=B-S local Z=l*math.exp(R*u)local y=S*math.exp(w*u)W=(Z+y)+V L=Z*R+y*w end local j=math.abs(L)<Q and math.abs(W-V)<R return{complete=j;value=j and V or W,velocity=L}end return S
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
