%This file to compute the static response of a truss
% Define a vector of indices for each spring
% Indices represent node numbers
clear;
close all;
%Define springs by node ends
link = [1 5;1 6;2 6;2 7;3 7;3 8;4 8;5 6;5 9;7 9;9 10;8 10;8 11;10 11;10 12;11 12];
% Define nodes by (x,y) coordinates
a = 1/2;
b = 3^0.5/2;

nodes=[0 0;1 0;2 0;3 0;-a -b;a -b;a+1 -b;a+2 -b;a -3*b;a+1 -3*b;a+3 -3*b;a+2 -5*b];

[u,v] = size(link);
[w,x] = size(nodes);


 
SpringV = zeros(u,v);
%Plot nodes, and define spring vectors
for i = 1:u,
    SpringV(i,:) = nodes(link(i,2),:)-nodes(link(i,1),:)
    plot([nodes(link(i,1),1);nodes(link(i,1),1)+SpringV(i,1)],[nodes(link(i,1),2);nodes(link(i,1),2)+SpringV(i,2)]);hold on;
end;



% Implement force balance equation for each node;

FB = zeros(2*w,u);
for j = 1:w,
    for k=1:u,
        if link(k,1)==j,
            FB(j,k)=FB(j,k)+SpringV(k,1)/(SpringV(k,:)*SpringV(k,:)')^0.5;
            FB(j+w,k)=FB(j+w,k)+SpringV(k,2)/(SpringV(k,:)*SpringV(k,:)')^0.5;
        end;
        if link(k,2)==j,
            FB(j,k)=FB(j,k)-SpringV(k,1)/(SpringV(k,:)*SpringV(k,:)')^0.5;
            FB(j+w,k)=FB(j+w,k)-SpringV(k,2)/(SpringV(k,:)*SpringV(k,:)')^0.5;
        end;
    end;
end;

%Now Force vectors for attachment points.
FA = zeros(2*w,8);

for j=1:4,
    FA(j,j) = 1;
    FA(j+w,j+4)=1;
end;

FB = [FB FA];

%Forces

F = zeros(2*w,1);
F(5) = -1;
F(9) = -2^0.5/2;
F(12+9) =-2^0.5/2;
F(11)=1;
F(24)=-1;

U = -FB\F;