local Y=.0001 local C=math.exp local O=math.cos local J=math.sin local n=math.sqrt local D={}local e={__index=D}function D.new(Y,C,O,J,n)local D={}D.p=Y D.v=C D.target=O D.angularFrequency=J or 10 D.dampingRatio=n or 1 return setmetatable(D,e)end function D.update(e,D)local p=e.angularFrequency local i=e.dampingRatio if p<Y then return end if i<0 then i=0 end local G=e.target local R=e.p-G local x=e.v if i>1+Y then local Y=-p*i local O=p*n(i*i-1)local J=Y-O local L=Y+O local I=C(J*D)local M=C(L*D)local F=(x-R*L)/(-2*O)local H=R-F e.p=(G+F*I)+H*M e.v=(F*J)*I+(H*L)*M elseif i>1-Y then local Y=C(-p*D)local O=x+p*R local J=R local n=(O*D+J)*Y e.p=G+n e.v=O*Y-n*p else local Y=p*i local L=p*n(1-i*i)local I=C(-Y*D)local M=O(L*D)local F=J(L*D)local H=R local h=(x+Y*R)/L e.p=G+I*(H*M+h*F)e.v=-I*((H*Y-h*L)*M+(H*L+h*Y)*F)end end return D


local EPSILON = 0.0001;

local exp = math.exp;
local cos = math.cos;
local sin = math.sin;
local sqrt = math.sqrt;

local spring = {};
local spring_mt = {__index = spring};

function spring.new(p0, v0, target, angularFrequency, dampingRatio)
	local self = {};
	self.p = p0;
	self.v = v0;
	self.target = target;
	self.angularFrequency = angularFrequency or 10;
	self.dampingRatio = dampingRatio or 1;
	return setmetatable(self, spring_mt);
end

function spring:update(dt)
	local aF = self.angularFrequency;
	local dR =  self.dampingRatio;
	
	if (aF < EPSILON) then return; end;
	if (dR < 0) then dR = 0; end;
	
	local epos = self.target;
	local dpos = self.p - epos;
	local dvel = self.v;
	
	if (dR > 1 + EPSILON) then
		local za = -aF * dR;
		local zb = aF * sqrt(dR*dR - 1);
		local z1 = za - zb;
		local z2 = za + zb;
		local expTerm1 = exp(z1 * dt);
		local expTerm2 = exp(z2 * dt);
		
		local c1 = (dvel - dpos*z2)/(-2*zb);
		local c2 = dpos - c1;
		self.p = epos + c1*expTerm1 + c2*expTerm2;
		self.v = c1*z1*expTerm1 + c2*z2*expTerm2;
	elseif (dR > 1 - EPSILON) then
		local expTerm = exp(-aF * dt);
		
		local c1 = dvel + aF*dpos;
		local c2 = dpos;
		local c3 = (c1*dt + c2)*expTerm;
		
		self.p = epos + c3;
		self.v = (c1*expTerm) - (c3*aF);
	else
		local omegaZeta = aF*dR;
		local alpha = aF*sqrt(1 - dR*dR);
		local expTerm = exp(-omegaZeta*dt);
		local cosTerm = cos(alpha*dt);
		local sinTerm = sin(alpha*dt);
		
		local c1 = dpos;
		local c2 = (dvel + omegaZeta*dpos) / alpha;
		self.p = epos + expTerm*(c1*cosTerm + c2*sinTerm);
		self.v = -expTerm*((c1*omegaZeta - c2*alpha)*cosTerm + (c1*alpha + c2*omegaZeta)*sinTerm);
	end
end

return spring;
