%% Analisi di una struttura reticolare iperstatica
%
%
%% Ordinamento delle variabili (solo promemoria)
% Le variabili in gioco saranno ordinate come segue:
%     F  = [FBy; FCx; FDx];
%     uf = [uBy; uCx; uDx];
%     N  = [N1;  N2;  N3;  N4;  N5;  N6];
%     d  = [d1;  d2;  d3;  d4;  d5;  d6]; <-- allungamenti nelle aste
%     R = [HA;  VA;  HB;  VC;  VD]; <-- reazioni vincolari
%     D = [DAx;  DAy;  DBx;  DCy;  DDy]; <-- cedimenti vincolari
%

% funzione per il calcolo del versore B-A
    vers = @(B,A) (B-A) / sqrt((B-A)'*(B-A));
% funzione per il calcolo del modulo B-A
    modulo = @(B,A) sqrt((B-A)'*(B-A));

%% 1. Definizione delle azioni/matrici di congruenza/matrice di rigidezza
%% Geometria

    syms h L1 L2 L El Ar % El = modulo elastico, Ar = area
    assume([0 0 0 0 0 0]<=[h L1 L2 L El Ar])
    % per verifica power point
        L1=L;
        L2=L;
        h=L;
% coordinate dei punti (componenti per colonna)
    A = [0; 0]; B = [0; h]; C = [L1+L2; h]; D = [L1; 0];

%% Azioni

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

syms alpha Dt % coeff. di dilatazione termica e Delta temperatura
DAnelast = alpha*Dt*[0; 0; modulo(D,B); 0; 0; -modulo(C,D)]; % allungamento anelastico delle aste (dist. termica)
DeltaEst  = [0; 0; 0; 0; 0];   % cedimenti vincolari esterni

%
%% Matrici di congruenza

% Matrice Cf
    vBA = vers(B,A); vDA = vers(D,A); vDB = vers(D,B);
    vCB = vers(C,B); vCA = vers(C,A); vCD = vers(C,D);
    vAx = [1; 0]; vAy = [0; 1]; vBx = [1; 0];
    vCy = [0; 1]; vDy = [0; 1];
    
    zero = [0; 0];
    
    Ctot = [-vBA',  vBA', zero', zero';...
            -vDA', zero', zero', vDA';...
            zero', -vDB', zero', vDB';...
            zero', -vCB',  vCB', zero';...
            -vCA', zero',  vCA', zero';...
            zero', zero',  vCD', -vCD'];
    IndLib  = [4, 5, 7];
    IndVinc = [1, 2, 3, 6, 8];
    Cf = Ctot(:,[IndLib])
    Cv = Ctot(:,[IndVinc])

%% Matrice di rigidezza assiale

    mBA = modulo(B,A); mDA = modulo(D,A); mDB = modulo(D,B);
    mCB = modulo(C,B); mCA = modulo(C,A); mCD = modulo(C,D);
    kel = El*Ar*diag([1/mBA, 1/mDA, 1/mDB, 1/mCB, 1/mCA, 1/mCD])
clear mBA mDA mDB mCB mCA mCD
clear vBA vDA vDB vCB vCA vCD vAx vAy vBx vCy vDy
clear Ctot IndVinc IndLib

%
%% Calcolo attraverso le equazioni del metodo degli spostamenti
%% 2. Calcoli preliminari (metodo degli spostamenti)

Fhat = Forze - Cf'*kel*Cv*DeltaEst + Cf'*kel*DAnelast
Krig = Cf'*kel*Cf

%% 3. Calcolo degli spostamenti nodali (metodo degli spostamenti)

uf = inv(Krig)*Fhat;

%% 4. Calcolo degli allungamenti assiali (metodo degli spostamenti)

dassiali = Cf*uf + Cv*DeltaEst;

%% 5. Calcolo degli sforzi normali (metodo degli spostamenti)

SfNorm = kel*(dassiali-DAnelast);

%% 6. Calcolo delle reazioni vincolari (metodo degli spostamenti)

Reaz = Cv'*SfNorm;

%
%% Cacolo

% sostituzione per lunghezze generiche
DaSost = [El,     Ar,  L1,   L2,   h,    L,    alpha,        Dt, P];
Valori = [210000, 100, 1000, 1000, 1000, 1000, 11.6*10^(-6), 10, 10000];

table({'N_1=';'N_2=';'N_3=';'N_4=';'N_5=';'N_6='},...
          double(subs(SfNorm,DaSost,Valori)),...
          {'N';'N';'N';'N';'N';'N'})
table({'H_A=';'V_A=';'H_B=';'V_C=';'V_D='},...
       double(subs(Reaz,DaSost,Valori)),...
       {'N';'N';'N';'N';'N'})
table({'delta_1=';'delta_2=';'delta_3=';'delta_4=';'delta_5=';'delta_6='},...
    double(subs(dassiali,DaSost,Valori)),...
    {'mm';'mm';'mm';'mm';'mm';'mm'})
table({ 'uBy'; 'uCx'; 'uDx'},...
    double(subs(uf,DaSost,Valori)),...
    {'mm';'mm';'mm'})