%
% lec2A.m
% generate the A matrix for the truss of figure 3
% proceed to construct, solve, and plot A^TKA x = f
%

data = [			% one row of data for each fiber, the
1 	4	-pi/2		% first two columns are node numbers
1	5	-pi/4		% while the third is the fiber angle
1	2	0
2	4	-3*pi/4
2	5	-pi/2
2	6	-pi/4
2	3	0
3	5	-3*pi/4
3	6	-pi/2
4	7	-pi/2
4	8	-pi/4
4	5	0
5	7	-3*pi/4
5	8	-pi/2
5	9	-pi/4
5	6	0
6	8	-3*pi/4
6	9	-pi/2
7	8	0
8	9	0];

A = zeros(20,18);		% initialize the A matrix

for j=1:20,

    m = data(j,1);
    n = data(j,2);
    theta = data(j,3);

    A(j,2*n-1) = cos(theta);	% direct implementation of general formula
    A(j,2*m-1) = -cos(theta);
    A(j,2*n) = sin(theta);
    A(j,2*m) = -sin(theta);

end

syms a cA	% this A is too big to view conveniently, so round the
		% zeros and ones and replace the bulky .7071 with `a'
for i=1:20,
for j=1:18,
    
    if abs(abs(A(i,j))-1/sqrt(2)) > .1,
       cA(i,j) = round(A(i,j));
    else
       if A(i,j) > .5,
          cA(i,j) = a;
       else
          cA(i,j) = -a;
       end
    end

end
end

disp(cA)		% display the cleaned A matrix (see the diary)
 
K = eye(20);
 
S = A'*K*A;

f = [-1 1 0 1 1 1 -1 0 0 0 1 0 -1 -1 0 -1 1 -1]';

pS = pinv(S);

x = pS*f;

% draw the nodes of the original and deformed specimens

t = 0:.01:1;		% a reference circle
mx = cos(2*pi*t);
my = sin(2*pi*t);

for h=1:3,		% translate and plot the reference circle

    plot(mx,my+(h-1)*6)			% original
    hold on
    plot(mx+x(1+(h-1)*6),my+(h-1)*6 - x(2+(h-1)*6),'--') 	% deformed
    plot(mx+6,my+(h-1)*6)
    plot(mx+6+x(3+(h-1)*6),my+(h-1)*6 - x(4+(h-1)*6),'--')
    plot(mx+12,my+(h-1)*6)
    plot(mx+12+x(5+(h-1)*6),my+(h-1)*6 - x(6+(h-1)*6),'--')

end

axis equal
axis off
hold off

