fourbar.m

function [angles] = fourbar(lengths, t1, t2, sig)

t2s = t2 - t1;

a = lengths(1);
b = lengths(2);
c = lengths(3);
d = lengths(4);

K1 = d/a;
K2 = d/c;
K3 = (a*a - b*b + c*c + d*d)/2/a/c;

A = cosd(t2s) - K1 - K2*cosd(t2s) + K3;
B = -2*sind(t2s);
C = K1 - (K2 + 1)*cosd(t2s) + K3;

t4s = 2*atand( (-B - sig*sqrt(B.*B - 4*A.*C)) /2./A); %multiplying by negative sig eliminates an if statement
t3s = asind( (c.*sind(t4s) - a.*sind(t2s) ) / b);

t4 = t4s + t1;
t3 = t3s + t1;

angles = [t3; t4];

end