nb = 16;
nv = 8;

barras=load('barras.txt')
nos=load('nos.txt')
xx=zeros(2*nv,1);
restric=zeros(2*nv,1);

% construcao de xx, vetor posicao inicial
% construcao de restric, informando restricoes

for i=1:nv
 xx(2*i-1)=nos(i,1);
 xx(2*i)=nos(i,2);
 massa(i)=nos(i,3);
 kode=nos(i,4);
 if (kode >= 10) 
   restric(2*i)=1;
   if (kode == 1 | kode == 11) 
    restric(2*i-1)=1;
   end
 end
end

%% plota estado atual
for i=1:nb
	xyloc(1,1) = xx(2*barras(i,1)-1);
	xyloc(1,2) = xx(2*barras(i,1));
	xyloc(2,1) = xx(2*barras(i,2)-1);
	xyloc(2,2) = xx(2*barras(i,2));
        plot(xyloc(:,1),xyloc(:,2));
        hold on;
end
axis([-1 2 -1 5]);
pause;

%%% calcula forcas internas no estado xx
intforces=zeros(nv*2,1);
barforces=zeros(2*2,1);
noloc=zeros(2,1);
%% compute internal forces
tensao=zeros(nb,1);
for i=1:nb
    noloc(1)=barras(i,1);noloc(2)=barras(i,2);
    xyloc(1,1) = xx(2*barras(i,1)-1);
    xyloc(1,2) = xx(2*barras(i,1));
    xyloc(2,1) = xx(2*barras(i,2)-1);
    xyloc(2,2) = xx(2*barras(i,2));
    l0=barras(i,3);
    dd=xyloc(1,:)-xyloc(2,:);
    ll=norm(dd);
    k=barras(i,4);
    auxcos=dd/ll;
    tens=k*(ll-l0);
    barforces(1)=-tens*auxcos(1);
    barforces(2)=-tens*auxcos(2);
    barforces(3)=-barforces(1);
    barforces(4)=-barforces(2);
    tensao(i)=tens;
    for ino=1:2
     for ico=1:2
       iunk=(noloc(ino)-1)*2+ico;
       iloc=(ino-1)*2+ico;
       intforces(iunk)=intforces(iunk)+barforces(iloc);
     end
    end
end

%%%
%% External forces
extforces=zeros(2*nv,1);
for i=1:nv
    extforces(2*i)=massa(i)*(-10);
end
extforces(15)=70;
%%%

%% 
  totforces=intforces+extforces;
%%
%% Restricoes
for i=1:2*nv
    if (restric(i) != 0) totforces(i)=0; end
end

% Aqui totforces tem a resultante total de forcas