%% Analisi di una struttura reticolare isostatica
%
%
%% Ordinamento delle variabili (solo promemoria)
% Le variabili in gioco saranno ordinate come segue:
%     F  = [FBx; FBy; FCx; FCy; FDx; FDy; FEx];
%     uf = [uBx; uBy; uCx; uCy; uDx; uDy; uEx];
%     N  = [N1;  N2;  N3;  N4;  N5;  N6;  N7];
%     d  = [d1;  d2;  d3;  d4;  d5;  d6;  d7]; <-- allungamenti nelle aste
%     R = [HA; VA; VE]; <-- reazioni vincolari
%     D = [DAx; DAy; DEy]; <-- cedimenti vincolari
%
%% 1. Definizione delle azioni/matrici di congruenza/matrice di rigidezza
%% Azioni

syms F % si indica genericamente con F l'intenzità della forza
       % il valore di riferiemnto del carico applicato
Forze  = [0; -F; 0; 0; F/2; -F; 0];

DAnelast = [0; 0; 0; 0; 0; 0; 0]; % allungamento anelastico delle aste (dist. termica)
DeltaEst  = [0; 0; 0];   % cedimenti vincolari

%% Matrici di congruenza

    syms h L El Ar% El = modulo elastico, Ar = area
    assume([0 0]<=[h L])
%     h = 2*L;
% coordinate dei punti (componenti per colonna)
    A = [0; 0]; B = [L; h]; C = [2*L; 0];
    D = [3*L; h]; E = [4*L; 0];
% Matrice Cf
    % funzione per il calcolo del versore B-A
    vers = @(B,A) (B-A) / sqrt((B-A)'*(B-A));
    vBA = vers(B,A); vCA = vers(C,A); vCB = vers(C,B);
    vDB = vers(D,B); vDC = vers(D,C); vEC = vers(E,C);
    vED = vers(E,D);
    vAy = [0; 1]; vEy = [0; 1]; vEx = [1; 0];
    zero = [0; 0];
    Cf = [vBA', zero', zero', 0;
          zero', vCA', zero', 0;
          -vCB', vCB', zero', 0;
          -vDB', zero', vDB', 0;
          zero', -vDC', vDC', 0;
          zero', -vEC', zero', vEC(1);
          zero', zero', -vED', vED(1)]
    Cv = [-vBA', 0;
          -vCA', 0;
          zero', 0;
          zero', 0;
          zero', 0;
          zero', vEC(2);
          zero', vED(2)]

%% Matrice di rigidezza assiale

    % funzione per il calcolo del modulo B-A
    modulo = @(B,A) sqrt((B-A)'*(B-A));
    mBA = modulo(B,A); mCA = modulo(C,A); mCB = modulo(C,B);
    mDB = modulo(D,B); mDC = modulo(D,C); mEC = modulo(E,C);
    mED = modulo(E,D);
    kel = El*Ar*diag([1/mBA, 1/mCA, 1/mCB, 1/mDB, 1/mDC, 1/mEC, 1/mED])
clear mBA mCA mCB mDB mDC mEC mED
clear vBA vCA vCB vDB vDC vEC vED

%
%% Risoluzione della struttura
%% 2.1 Calcolo degli sforzi normali

SfNorm = inv(Cf')*Forze;

%% 2.2 Calcolo delle reazioni vincolari esterne

Reaz = Cv'*SfNorm;

%% 3. Calcolo degli allungamenti assiali

dassiali = inv(kel)*SfNorm + DAnelast;

%% 4. Calcolo degli spostamenti nodali

uf = inv(Cf)*(dassiali - Cv*DeltaEst);

%
%% Cacolo

DaSost = [El,Ar,L,h,F];
Valori=[210000,10,1000,2000,1000];

table({'N_1=';'N_2=';'N_3=';'N_4=';'N_5=';'N_6=';'N_7='},...
          double(subs(SfNorm,DaSost,Valori)),...
          {'N';'N';'N';'N';'N';'N';'N'})
table({'H_A=';'V_A=';'V_E='},...
       double(subs(Reaz,DaSost,Valori)),...
       {'N';'N';'N'})
table({'delta_1=';'delta_2=';'delta_3=';'delta_4=';'delta_5=';'delta_6=';'delta_7='},...
    double(subs(dassiali,DaSost,Valori)),...
    {'mm';'mm';'mm';'mm';'mm';'mm';'mm'})
table({ 'uBx'; 'uBy'; 'uCx'; 'uCy'; 'uDx'; 'uDy'; 'uEx'},...
    double(subs(uf,DaSost,Valori)),...
    {'mm';'mm';'mm';'mm';'mm';'mm';'mm'})