% 16.210 - PS2 Matlab File.
% Prepared by Dora Tzianetopoulou / Spring 2000

%-----------------------------------------------------------------

% PROBLEM 1.

% We use the PVW to find the displacement at joint 3.

% Equilibrium equations in matrix form: A*F = Q
% The A matrix as we have already found is:
% F = [F12;F23;F34;F45;F56;F67;F17;F27;F24;F47;Ax;Ay;By] =
%	  [F1;F10;F9;F8;F7;F4;F2;F3;F6;F5;Ax;Ay;By]


clear all

syms l EA real;

A=zeros(13); 

x=1/sqrt(2);

A(1,1)=x; A(1,7)=1; A(1,11)=-1;
A(2,1)=x; A(2,12)=-1; A(3,1)=-x;
A(3,2)=x; A(3,9)=1; A(4,1)=-x;
A(4,2)=x; A(4,8)=-1; A(5,2)=x;
A(6,2)=x; A(6,3)=1; A(7,4)=x;
A(7,9)=-1; A(7,10)=-x; A(8,3)=1;
A(8,4)=-x; A(8,10)=-x; A(9,4)=x;
A(9,5)=1; A(10,4)=x; A(10,13)=1;
A(11,5)=1; A(11,6)=-1; A(12,6)=1;
A(12,7)=-1; A(12,10)=x; A(13,8)=1;A(13,10)=x;


% Q1 and Q2 are the external loads.
Q=sym('[0;0;0;Q2;Q1;0;0;0;0;0;0;0;0]');



% To get my results in symbolic form, I must have A 
% in symbolic form as well.
symA=sym(A);
invA=inv(symA);

% Force and Reaction Vector. 
F=invA*Q;

% I only need the internal member loads.
Fd=F(1:10,1);

% ----------------------------------------------------------------
% ----------------------------------------------------------------
%
% RESULTS
% 
% Real Loading System.
% Member Loads' Vector Fv1 due to the real applied loads Q1 and Q2.
% 
% 
% Fd =
%  
% [ -2/3*2^(1/2)*Q2+2/3*2^(1/2)*Q1]
% [                     2^(1/2)*Q1]
% [                            -Q1]
% [ -1/3*2^(1/2)*Q2-2/3*2^(1/2)*Q1]
% [                  1/3*Q2+2/3*Q1]
% [                  1/3*Q2+2/3*Q1]
% [                  2/3*Q2+1/3*Q1]
% [                 -1/3*Q2+1/3*Q1]
% [                 -2/3*Q2-1/3*Q1]
% [  1/3*2^(1/2)*Q2-1/3*2^(1/2)*Q1]


% Now we can proceed with the PVW. That is we can set
% in the above expressions for the member loads Q2 = 0 and Q1 = 1.

% Qvx=1 and Qvy=1 are the virtual loads.
% Qv1=sym('[0;0;0;0;1;0;0;0;0;0;0;0;0]'); I don't really use this as you 
% will see in a while.
Qv2=sym('[0;0;0;0;0;-1;0;0;0;0;0;0;0]');

% For the horizontal displacement all I must do is sustitute 
% in the member loads' matrix, Fd, Q1=1 and Q2=0.
  Q1=-1;Q2=0;
Fv1=subs(Fd);
Fv1=Fv1(1:10,:);

% ----------------------------------------------------------------
% ----------------------------------------------------------------
%
% RESULTS
% 
% Horizontal Virtual Loading case.
% Member Loads' Vector Fv1 due to the virtual load Qv1.
% 
% 
% Fv1 =
% 
%     0.9428
%     1.4142
%    -1.0000
%    -0.9428
%     0.6667
%     0.6667
%     0.3333
%     0.3333
%    -0.3333
%    -0.4714


% For the vertical displacment I must follow the same method 
% we used in PS1 to find the member loads. The geometry matrix 
% A is the same but this time we use Qv2 instead of Q for the 
% external loads' vector. So, we get the member loads' vector, 
% Fv2, with Q1=1 and Q2=0 and a unit vertical load applied at T.
  
Fv2=invA*Qv2;
Fv2=Fv2(1:10,:);

% ----------------------------------------------------------------
% ----------------------------------------------------------------
%
% RESULTS
% 
% Horizontal Virtual Loading case.
%  Member Loads' Vector Fv2 due to the virtual load Qv2.
% 
% Fv2 =
%  
% [ 1/3*2^(1/2)]
% [           0]
% [           1]
% [ 2/3*2^(1/2)]
% [        -2/3]
% [        -2/3]
% [        -1/3]
% [        -1/3]
% [         1/3]
% [ 1/3*2^(1/2)] 

% This is a diagonal vector with diagonal entries the lengths 
% of the respective member loads.
L=zeros(10,10);
L(1,1)=1/x; L(2,2)=1/x; L(3,3)=1; L(4,4)=1/x; L(5,5)=1;
L(6,6)=1; L(7,7)=1; L(8,8)=1; L(9,9)=1; L(10,10)=1/x;
L=l*L;

% Now I can use the formula from Handout 2.
% Horizontal displacent u:
h1=L*Fv1;
h2=Fd.*h1;
u=sum(h2(:,1))/EA
u=collect(u)

% Vertical displacent v:
h3=L*Fv2;
h4=Fd.*h3;
v=sum(h4(:,1))/EA;
v=collect(v)

% ----------------------------------------------------------------
% ----------------------------------------------------------------
%
% RESULTS
% 
% Horizontal displacement u.
%  
% u =
%  
% (-2/3*2^(1/2)*Q2+4*2^(1/2)*Q1+20/9*Q1+7/9*Q2)/EA*l^2
% 
%
% Vertical displacement v.
%  
% v =
%
% (-2/3*2^(1/2)*Q2-2/3*2^(1/2)*Q1-20/9*Q1-7/9*Q2)/EA*l^2
%  
%  These displacements are the same with the ones we found in PS1.

%-------------------------------------------------------------------
% 
% PROBLEM 2

% Again we use the PVW twice to find the horizontal and vertical 
% components of the displacement at T. In the handwriten part of 
% the solutions I derive the formula needed for this case.
% Also, we can use Fv1 and Fv2 from Problem 1.
% epsilon = alpha * DT
syms Alpha DT real;
constant=Alpha*DT;

% Horizontal displacement u for the temperature difference case.
sum1=sum(h1(:,1));
u=constant*sum1

% Vertical displacement v for the temperature difference case.
sum2=sum(h3(:,1));
v=constant*sum2

% ----------------------------------------------------------------
% ----------------------------------------------------------------
%
% RESULTS
% 
% Horizontal displacement u.
%  
% u =
%  
% 2*Alpha*DT*l
% 
%
% Vertical displacement v.
%  
% v =
%  
% 2*Alpha*DT*l

%-----------------------------------------------------------------