local m={}local o={__index=m}function m.new(m,j,h,F,y)local w={}w.p=m w.v=j w.a=0 w.target=h w.angularFrequency=F or.1 w.dampingRatio=y or.1 return setmetatable(w,o)end function m.update(o,m)local j=o.angularFrequency local h=o.dampingRatio if j<EPSILON then return end if h<0 then h=0 end local F=o.target local y=o.p-F local w=o.v if h>1+EPSILON then local b=-j*h local u=j*sqrt(h*h-1)local C=b-u local p=b+u local n=exp(C*m)local c=exp(p*m)local K=(w-y*p)/(-2*u)local Q=y-K o.p=(F+K*n)+Q*c o.v=(K*C)*n+(Q*p)*c elseif h>1-EPSILON then local h=exp(-j*m)local b=w+j*y local u=y local C=(b*m+u)*h o.p=F+C o.v=b*h-C*j else local b=j*h local u=j*sqrt(1-h*h)local C=exp(-b*m)local p=cos(u*m)local n=sin(u*m)local c=y local K=(w+b*y)/u o.p=F+C*(c*p+K*n)o.v=-C*((c*b-K*u)*p+(c*u+K*b)*n)end end return m


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

function spring.new(p0, v0, target, angularFrequency, dampingRatio)
	local self = {};
	self.p = p0;
	self.v = v0;
	self.a = 0;
	self.target = target;
	self.angularFrequency = angularFrequency or 0.1;
	self.dampingRatio = dampingRatio or 0.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
		-- critical
		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;
