%% Lab 3 Övning
% Table of Contents
% Uppgift 1 - Trusses (fackverk) as linear system of equations
% Uppgift 3
% Uppgift 4 -
%% Uppgift 1 - Trusses (fackverk) as linear system of equations
% Each edge and node is enumerated. For each node, set up equations for external forces in vertical and horizontal directions. All the forces in the structure must be balanced, and we thus solve for the net force in each edge (tension or compression), such that the structure is static.
% The node at the fixed anchor point on the left is not included (unless we make it a free-body problem, inwhich normal-forces are counted). The roller on the right has it's normal-force in the vertican direction, so that isn't counted, only the horizontal forces exerted by the structure.
%
% To solve this efficiently in Matlab, we first find the equations, then fill a sparse matrix and use backslash \.

close all

u = 1/sqrt(2); 
alpha = atan(1/2);
rad=[1 1 2 3 3 4 4 5 5 5 6 6 6 7 7 8 9 9];
kol=[2 6 3 1 5 1 3 4 5 8 4 7 8 6 9 7 8 9];
ele=[-1 1 1 -u 1 -u -1 -cos(alpha) -1 u -sin(alpha) -1 -u -1 1 1 -u -1];
A = sparse(rad,kol,ele);


blep = [0, 0, 0, 0, 0, 0, 0, 0, 0;
0, 0, 0, 0, 0, 0, 0, 0, 0;
0, 0, 0, 0, 0, 0, 0, 0, 0;
0, 0, 0, 0, 0, 0, 0, 0, 0;
0, 0, 0, 0, 0, 0, 0, 0, 0;
0, 0, 0, 0, 0, 0, 0, 0, 0;
0, 0, 0, 0, 0, 0, 0, 0, 0;
0, 0, 0, 0, 0, 0, 0, 0, 0;
0, 0, 0, 0, 0, 0, 0, 0, 0]

b = [0 10 0 0 0 0 0 15 0]';
x = A\b
spy(A)



%
%% Uppgift 3

load('avkylningseffekt.mat')



%
%% Uppgift 4 -
%
% We see that q(x) and I(x) are functions only of x.

close all
clear all

n = 30; % INNER grid points
L = 1;

E = 2e11; % Given

I = @(x) (2-x).*10^(-7);
q = @(x) (10 - 9.*x).*500;


h = L/(n+1); % n+1 intervals in the system when n is the INNER grid.

x = ((1:n).*h)'; % x positions of the INNER grid.

qv = q(x);
e = ones(n, 1);
A = spdiags([-e, 2.*e, -e], [-1, 0, 1], n,n);
Afull = full(A);
Mv = A\(h^2.*qv)

vv = A\((h^2)./(E.*I(x)).*Mv)

V = [0;vv;0]

M = [0; Mv; 0]

X = [0; x; L]

plot(X, V)
plot(X, M)

% Official solution
C = chol(A);
qv = h^2*q(x);
M = C\(C'\qv);
r = h^2*M./(E*I(x));
v = C\(C'\r);

V = [0; v; 0];

plot(X, V)