function [F] = calcIBForce(XIB,LINKS,P1,P2,STIFF,REST,NIB,HINGES,B1,B2,...
               B3,BENDSTIFF,C,YT,TARGETSTIFF);
%
% [F] = calcIBForce(XIB,LINKS,P1,P2,STIFF,REST,NIB,HINGES,B1,B2,
%       B3,BENDSTIFF,C,YT,TARGETSTIFF);
%
%  Calculates the elastic force at each IB point
%  
%  Returns:
%
%     F      = total elastic force on each IB point
%    
%  Input:
%     XIB     = IB point positions
%     LINKS   = list of number of Lagrangian points
%     P1,P2   = indices of pairs of IB points joined by link l
%     STIFF   = stiffness of link l
%     REST    = rest length of link l
%     NIB     = number of IB points
%     dq      = IB point spacing
%     YT      = target point locations
%
%
      F = zeros(NIB,2);
% 'Link' spring forces
%
      tmp = XIB(P2(LINKS),:)-XIB(P1(LINKS),:);
% compute link lengths
      lnths = sqrt(dot(tmp,tmp,2));
% compute unit tangent vector
      tau   = [tmp(:,1)./lnths tmp(:,2)./lnths];
% stiffnesses
      stf   = STIFF(LINKS);
% resting lengths
      rst   = REST(LINKS);
% link force values
      f1tmp = stf.*(lnths-rst).*tau(:,1);
      f2tmp = stf.*(lnths-rst).*tau(:,2);
% store link forces in appropriate locations of F
      F(P1(LINKS),1) = F(P1(LINKS),1) + f1tmp;
      F(P1(LINKS),2) = F(P1(LINKS),2) + f2tmp;
      F(P2(LINKS),1) = F(P2(LINKS),1) - f1tmp;
      F(P2(LINKS),2) = F(P2(LINKS),2) - f2tmp;
%
% 'Hinge' forces
%
      tmp23 = XIB(B3(HINGES),:)-XIB(B2(HINGES),:);
      tmp12 = XIB(B2(HINGES),:)-XIB(B1(HINGES),:);
      bstf = BENDSTIFF(HINGES);
      curv = C(HINGES);
      factor= bstf.*( tmp23(:,1).*tmp12(:,2) - tmp12(:,1).*tmp23(:,2) - curv);
%
% store hinge forces in appropriate locations of F
      F(B1(HINGES),1)= F(B1(HINGES),1) + factor.*(- tmp23(:,2));
      F(B1(HINGES),2)= F(B1(HINGES),2) + factor.*(  tmp23(:,1));
      F(B2(HINGES),1)= F(B2(HINGES),1) + factor.*(  tmp23(:,2)+tmp12(:,2));
      F(B2(HINGES),2)= F(B2(HINGES),2) + factor.*(-(tmp23(:,1)+tmp12(:,1)));
      F(B3(HINGES),1)= F(B3(HINGES),1) + factor.*(- tmp12(:,2));
      F(B3(HINGES),2)= F(B3(HINGES),2) + factor.*(  tmp12(:,1));
%
% calculate target forces
%
      F(:,1) = F(:,1) + TARGETSTIFF.*(YT(:,1) - XIB(:,1));
      F(:,2) = F(:,2) + TARGETSTIFF.*(YT(:,2) - XIB(:,2));