% Two-element truss, truss1,  3 nodes
clear                       % variables
format short e
% format long               % double precision

dtor = pi/180;              % convert degrees to radians
force = 10e3;               % applied force
fang = (180+30)*dtor;       % angle of force
px = (force*cos(fang))
py = (force*sin(fang))
c1 = cos(0*dtor);
s1 = sin(0*dtor);
c2 = cos((90+45)*dtor);
s2 = sin((90+45)*dtor);

ex = [ 29 11 ]*1e6          % elastic modulus
area = [ 2 1.5 ]
len = [ 4  4*sqrt(2) ]*12
k12 = area(1)*ex(1)/len(1)  % element 1
k12c = k12*c1^2
k12s = k12*s1^2
k12cs = k12*c1*s1
k23 = area(2)*ex(2)/len(2)  % element 2
k23c = k23*c2^2
k23s = k23*s2^2
k23cs = k23*c2*s2
k33 = k12c + k23c
k44 = k12s + k23s
k34 = k12cs + k23cs
beta = k12 * 1e7
stiff = [ k12c   k12cs  -k12c   -k12cs    0      0  % stiff matrix
          k12cs  k12s   -k12cs  -k12s     0      0
         -k12c  -k12cs   k33     k34    -k23c  -k23cs
         -k12cs -k12s    k34     k44    -k23cs -k23s
          0        0    -k23c   -k23cs   k23c   k23cs
          0        0    -k23cs  -k23s    k23cs  k23s  ]
df = [ 0*beta; 0*beta; px; py; 0*beta; 0*beta ] % combined d, F
zienk = stiff;
zienk(1,1)=zienk(1,1)+beta;
zienk(2,2)=zienk(2,2)+beta;
zienk(5,5)=zienk(5,5)+beta;
zienk(6,6)=zienk(6,6)+beta
format short e
displ = inv(zienk) * df;     % displacements
force = stiff * displ;       % force vector

% print the results
fprintf('\n    Two-element truss 6x6, truss1')
fprintf('\n    Your name, MENG 421, ')
disp(date)
fprintf('\n    Node   X Displ      FX        Y Displ      FY  ')
fprintf('\n           inches       Kip       inches       Kip \n')
for i = 1:2:length(df)
  disp(sprintf('    %2.0f   %10.2e %7.2f    %11.2e %7.2f',...
   (i+1)/2,displ(i),force(i)/1000,displ(i+1),force(i+1)/1000))
end

%    Two-element truss 6x6, truss1
%    Your name, MENG 421, 05-Feb-2004

%    Node   X Displ      FX        Y Displ      FY
%           inches       Kip       inches       Kip
%     1   -1.13e-009   13.66      0.00e+000    0.00
%     2   -1.13e-002   -8.66     -5.24e-002   -5.00
%     3    4.14e-010   -5.00     -4.14e-010    5.00