%
%  macm 202 -- 09 mar 03 -- djm
%
%  w10solve.m:  truss solve
%     - run w10mytruss#.m and w10setup.m first

%  truss equations -- make the matrix

matrix = zeros(Nstrut+3,Nstrut+3);
vector = zeros(Nstrut+3,1);

%  body forces -- b^j terms

matrix(1,1) = 1;
matrix(2,2) = 1;
matrix(4,3) = 1;

%  struts -- cos & sin alpha coefficients

for j = 1:Nnodes
    for k = 1:Nnodes
        if(struts(j,k)~=0)
            matrix(2*j-1,3+struts(j,k)) = cos(angles(j,k));
            matrix(2*j  ,3+struts(j,k)) = sin(angles(j,k));
            
            vector(2*j-1) = loads(j,1);
            vector(2*j  ) = loads(j,2);
        end
    end
end

%  do the math -- matrix*solution = -vector !!

solution = -matrix\vector;

%  plot, plot, plot some stuff

figure(2);clf
plot(nodes(:,1),nodes(:,2),'*')
hold on
plot(nodes(1:2,1),nodes(1:2,2),'o')
axis([min(nodes(:,1))-1.5 max(nodes(:,1))+1 ...
        min(nodes(:,2))-1.5 max(nodes(:,2))+1])
axis equal
shift = 1/8;
for j = 1:Nnodes
    text(nodes(j,1)+shift,nodes(j,2)+shift,num2str(j),'color','k')
end
    
for j = 1:Nnodes
    for k = j:Nnodes
        if (join(j,k)==1)
            plot([nodes(j,1) nodes(k,1)],[nodes(j,2) nodes(k,2)],'k')
            label = num2str(solution(3+struts(j,k)),'%5.2f');
            xmid = (nodes(k,1)+nodes(j,1))/2;
            ymid = (nodes(k,2)+nodes(j,2))/2;
            if (solution(3+struts(j,k))>0)
                text(xmid+shift,ymid+shift,label,'color','b')
            else
                text(xmid+shift,ymid+shift,label,'color','r')
            end
        end    
    end
end

%  support nodes

sn = [1 1 2];  snang = [-pi 3*pi/2 3*pi/2];  sh = [-2 1 1]*shift;
for j=1:3
    plot(nodes(sn(j),1)+[0 6*shift*cos(snang(j))], ...
         nodes(sn(j),2)+[0 6*shift*sin(snang(j))],'k');
    label = num2str(solution(j),'%5.2f');
    
    if (solution(j)>=0)
        text(nodes(sn(j),1)+6*shift*cos(snang(j))+sh(j), ...
             nodes(sn(j),2)+6*shift*sin(snang(j))-shift,label,'color','b');
    else
        text(nodes(sn(j),1)+6*shift*cos(snang(j))+sh(j), ...
             nodes(sn(j),2)+6*shift*sin(snang(j))-shift,label,'color','r');
    end
end

%  loads

for j = 1:Nnodes
    if (loadR(j)~=0)
        label = ['(' num2str(loads(j,1),'%5.2f') ',' ...
                     num2str(loads(j,2),'%5.2f') ')'];
        plot(nodes(j,1)+[0 3*shift*cos(loadA(j))], ...
             nodes(j,2)+[0 3*shift*sin(loadA(j))],'g');
        text(nodes(j,1)+  shift, ...
             nodes(j,2)-4*shift,label,'color','g');
    end
end

title('\bf solved truss (blue=tension, red=compression)')