temp=size(joint);
njoint=temp(1,1);
temp=size(member);
nmember=temp(1,1);
temp=size(constraints);
nconst=temp(1,1);
temp=size(force);
nforce=temp(1,1);
if (njoint*2)==(nconst+nmember)
cmat=zeros(njoint*2,njoint*2);
fvec=zeros(njoint*2,1);
for I=1:nmember
 j1=member(I,1);
 j2=member(I,2);
 x1=joint(j1,1);
 y1=joint(j1,2);
 x2=joint(j2,1);
 y2=joint(j2,2);
 length=sqrt((x2-x1)^2+(y2-y1)^2);
 dcx=(x2-x1)/length;
 dcy=(y2-y1)/length;
 cmat(2*(j1-1)+1,I)=dcx;
 cmat(2*(j1-1)+2,I)=dcy;
 cmat(2*(j2-1)+1,I)=-dcx;
 cmat(2*(j2-1)+2,I)=-dcy;
end
for I=1:nconst
 jc=constraints(I,1);
 xyc=constraints(I,2);
 cmat(2*(jc-1)+xyc,nmember+I)=1;
end
for I=1:nforce
 jf=force(I,1);
 xf=force(I,2);
 yf=force(I,3);
 fvec(2*(jf-1)+1,1)=-xf;
 fvec(2*(jf-1)+2,1)=-yf;
end
if(det(cmat)==0)
   fprintf('Structure is a mechanism');
   fprintf('')
else
   solution=inv(cmat)*fvec;
end
clear I j1 j2 jc jf x1 y1 x2 y1
clear nforce njoint nconst xyc
clear dcx dcy length nmember temp
clear xf y2 yf
else
fprintf('Structure is either unstable or overconstrained.')
end