% This script reads in a truss specified by a file
% containing the following information
% 
% nn  - the number of nodes
% x y - nn lines containing the (x,y) coordinates of the nodes
% nb  - the number of bars
% a b - the n lines containing the pair of nodes defining each bar.

% Daniel D. Warner
% September 2, 2007

name = input('Enter the name of the file defining the truss: ');
fid = fopen(name);
% Read in the number of nodes
nn = fscanf(fid,'%d',1)
% Read in the list of nodes
nodes = fscanf(fid,'%f %f',[2 nn])

% Read in the number of bars
nb = fscanf(fid,'%d',1)
% Read in the list of bars
bars = fscanf(fid,'%d %d',[2 nb])

% Plot the nodes
xn = nodes(1,:);
yn = nodes(2,:);
plot(xn,yn,'ok')
hold on
% Set the view
xmin = min(xn)-5;
xmax = max(xn)+5;
ymin = min(yn)-5;
ymax = max(yn)+5;
axis([xmin xmax ymin ymax]);

% Plot the bars
for k=1:nb
    % Get the index of the source node
    s = bars(1,k);
    % Get the index of the target node
    t = bars(2,k);
    plot([xn(s) xn(t)],[yn(s) yn(t)],'k')
end
pause

% Set up the structural matrix (the bar-node matrix)
A0 = zeros(nb,2*nn);
for k=1:nb
    % Get the index of the source node
    s = bars(1,k);
    % Get the index of the target node
    t = bars(2,k);
    % Calculate the length of the bar
    Lb = sqrt((xn(t)-xn(s))^2 + (yn(t)-yn(s))^2);
    % Calculate cos(theta) for the bar
    cb = (xn(t)-xn(s))/Lb;
    % Calculate sin(theta) for the bar
    sb = (yn(t)-yn(s))/Lb;
    A0(k,2*s-1) = cb;
    A0(k,2*s)   = sb;
    A0(k,2*t-1) = -cb;
    A0(k,2*t)   = -sb;
end
A0;

% Read in the number of constraints
% These are displacements that must be zero.
nc = fscanf(fid,'%d',1)
% Read in the list of constraints
% The first number is the node the second number
% indicates whether the x (0) or y (1) displacement
% is constrained to be zero.
cd = fscanf(fid,'%d %d',[2 nc])
% Record the active displacements
active = ones(1,2*nn);
for k=1:nc
    indx = 2*cd(1,k)+cd(2,k)-1;
    active(indx) = 0;
end
active;

% Read in the material constant
c = fscanf(fid,'%f',1);
C = c*eye(nb);

% Read in the number of applied loads
nf = fscanf(fid,'%d',1)
% Read in the list of applied loads
% The first number is the node the second and third
% numbers are the loads in the x and y directions,
% respectively.
ff = fscanf(fid,'%d %f %f',[3 nf])
% Set up the vector of applied loads and check for
% consistency
f0 = zeros(2*nn,1);
for k=1:nf
    indx = 2*ff(1,k)-1;
    f0(indx)   = ff(2,k);
    f0(indx+1) = ff(3,k);
    fprintf('Node: %d %12.4e %12.4e',ff(:,k));
end
pause
clc
% Delete the columns in the structural matrix
% and the components in the vector of applied
% loads that are irrelevant because of the
% constraints at the fixed nodes.
A = zeros(nb,2*nn-nc);
f = zeros(2*nn-nc,1);
j = 1;
for k=1:2*nn
    if active(k)
        A(:,j) = A0(:,k);
        f(j) = f0(k);
        j = j+1;
    end
end
A
% f'
tal = norm(f,1);
% fprintf('The total applied load is: %12.3e\n',tal);
fprintf('The stiffness matrix is:\n');
K = A'*C*A
det(K)
fprintf('The displacements are:\n')
d = K\f
e = -A*d
y = C*e
% Redraw the truss with the nodes displaced
% and the bars under tension in red
% and bars under compression in blue.
dnodes = nodes;
j = 1;
for k=1:nn
    kx = 2*k-1;
    ky = 2*k;
    if active(kx)
        dnodes(1,k) = dnodes(1,k) + d(j);
        j = j+1;
    end
    if active(ky)
        dnodes(2,k) = dnodes(2,k) + d(j);
        j = j+1;
    end
end
dnodes
pause
% Plot the displaced nodes
xn = dnodes(1,:);
yn = dnodes(2,:);
plot(xn,yn,'og')

% Plot the bars
for k=1:nb
    % Get the index of the source node
    s = bars(1,k);
    % Get the index of the target node
    t = bars(2,k);
    if (0.001*tal < y(k))
        % Bar is under tension
        plot([xn(s) xn(t)],[yn(s) yn(t)],'r')
    elseif (y(k) < -0.001*tal)
        % Bar is under compression
        plot([xn(s) xn(t)],[yn(s) yn(t)],'b')
    else
        plot([xn(s) xn(t)],[yn(s) yn(t)],'g')
    end 
end

% Calculate the reactive forces
f0r = (A0')*y;
% Contains the external loads and the reactive
% forces at the fixed nodes
fr = zeros(nc,1);
j = 1;
for k=1:2*nn
    if ~active(k)
        fr(j) = -f0r(k);
        j = j+1;
    end
end
fprintf('The reactive forces\n')
fr

% clean up
hold off
fclose(fid);