//
// Esercitazione 6
//
// ***************************
//
//     Procedure SI ed SA
//
// ***************************
//
function x = SI(T,c)
//
// Per T matrice N x N triangolare superiore invertibile e c colonna
// di N componenti, calcola la colonna x (di N numeri reali) che approssima la
// soluzione del sistema di equazioni lineari:
//
//       T z = c
//
N = size(T,'r');
//
x = zeros(N,1);
//
x(N) = c(N) / T(N,N);
for k = N-1:-1:1,
    x(k) = (c(k) - T(k,k+1:N)*x(k+1:N)) / T(k,k);
end;
endfunction
//
function x = SA(T,c)
//
// Per T matrice N x N triangolare inferiore invertibile e c colonna
// di N componenti, calcola la colonna x (di N numeri reali) che approssima la
// soluzione del sistema di equazioni lineari:
//
//       T z = c
//
N = size(T, 'r');
//
x = zeros(N,1);
//
x(1) = c(1) / T(1,1);
for k = 2:N,
    x(k) = (c(k) - T(k,1:k-1)*x(1:k-1)) / T(k,k);
end
endfunction
//
// *************************************************
//
//     Procedure EquazioniNodali e SoluzioneRete
//
// *************************************************
//
function [G,J] = EquazioniNodali(Rete)
    //
    // Determina la matrice G ed il termine noto J del sistema delle
    // equazioni nodali per il circuito di sole resistenze e generatori di
    // corrente, con N + 1 nodi di numero 0,...,N, descritto nella lista Rete.
    //
    // Rete: lista di due matrici che descrive il circuito.
    // Rete(1): matrice NRes x 3 ad elementi reali con righe, una per ciascuna
    //          resistenza del circuito, di formato
    //               valore resistenza , nodo 1 , nodo 2
    // Rete(2): matrice NGen x 3 ad elementi reali con righe, una per ciascun
    //          generatore di corrente del circuito, di formato
    //               valore generatore , nodo 1 , nodo 2
    //
    // G: matrice N x N ad elementi reali.
    // J: colonna di N numeri reali.
    //
    // Determina il numero di nodi della rete, di resistenze e di generatori
    [Res, Gen] = Rete(:);
    Nodi = union(Res(:,2:3), Gen(:,2:3));
    N = max(Nodi);
    NRes = size(Res,'r');
    NGen = size(Gen,'r')
    // Inizializza G e J
    G = zeros(N,N);
    J = zeros(N,1);
    // Costruisce G analizzando Res una riga alla volta
    for k = 1:NRes,
        riga = Res(k,:);
        valore = riga(1); nodo1 = riga(2); nodo2 = riga(3);
        if nodo1 ~= 0 & nodo2 ~= 0 then
           G([nodo1,nodo2],[nodo1,nodo2]) = ...
           G([nodo1,nodo2],[nodo1,nodo2]) + ...
               [1,-1 ;...
               -1, 1] / valore;
           else if nodo1 ~= 0 then
                   G(nodo1,nodo1) = G(nodo1,nodo1) + 1 / valore; end;
                if nodo2 ~= 0 then
                   G(nodo2,nodo2) = G(nodo2,nodo2) + 1 / valore; end;
        end;
    end;
    // Costruisce J analizzando Gen una riga alla volta
    for k = 1:NGen,
        riga = Gen(k,:);
        valore = riga(1); nodo1 = riga(2); nodo2 = riga(3);
        if nodo1 ~= 0 & nodo2 ~= 0 then
           J([nodo1,nodo2]) = ...
               J([nodo1,nodo2]) +  [-1;1] * valore;
           else if nodo1 ~= 0 then
                   J(nodo1) = J(nodo1) - valore; end;
                if nodo2 ~= 0 then
                   J(nodo2) = J(nodo2) + valore; end;
        end;
    end
endfunction
//
function [E,V,I] = SoluzioneRete(Rete)
    //
    // Determina tutte le tensioni e correnti del circuito, con N + 1
    // nodi di numero 0,...,N ed R rami, descritto nella lista Rete.
    //
    // Rete: lista di due matrici che descrive il circuito.
    // Rete(1): matrice NR x 3 ad elementi reali con righe, una per ciascuna
    //          resistenza del circuito, di formato
    //               valore resistenza , nodo 1 , nodo 2.
    // Rete(2): matrice NG x 3 ad elementi reali con righe, una per ciascun
    //          generatore di corrente del circuito, di formato
    //               valore generatore , nodo 1 , nodo 2
    //
    // I rami del circuito sono numerati considerando, in ordine, per primi gli
    // NR rami contenenti resistenze poi gli NG rami contenenti generatori di
    // corrente.
    //
    // E: colonna di N numeri reali. E' il vettore delle tensioni di nodo.
    // V: colonna di R numeri reali. E' il vettore delle tensioni di ramo.
    // I: colonna di R numeri reali. E' il vettore delle correnti di ramo.
    //
    // Determina i dati delle equazioni nodali
    [G, J] = EquazioniNodali(Rete);
    // (1) fattorizzazione di G
    [S,D,P] = lu(G);
    // (2) soluzione del sistema S c = PJ
    c = SA(S, P*J);
    // (3) soluzione del sistema D x = c
    E = SI(D,c);
    // Determina le grandezze di ramo del circuito.
    [Res, Gen] = Rete(:);
    Nodi = union(Res(:,2:3), Gen(:,2:3));
    N = max(Nodi);
    NRes = size(Res,'r');
    NGen = size(Gen,'r')
    R = NRes + NGen;
    // Inizializza V
    V = zeros(R,1);
    // Costruisce V analizzando Res e Gen una riga alla volta
    for k = 1:NRes,
        riga = Res(k,:);
        nodo1 = riga(2); nodo2 = riga(3);
        if nodo1 ~= 0 & nodo2 ~= 0 then
           V(k) = E(nodo1) - E(nodo2);
           else if nodo1 ~= 0 then V(k) = E(nodo1); end;
                if nodo2 ~= 0 then V(k) = - E(nodo2); end;
        end;
    end;
    for k = 1:NGen,
        riga = Gen(k,:);
        nodo1 = riga(2); nodo2 = riga(3);
        if nodo1 ~= 0 & nodo2 ~= 0 then
           V(NRes + k) = E(nodo1) - E(nodo2);
           else if nodo1 ~= 0 then V(NRes + k) = E(nodo1); end;
                if nodo2 ~= 0 then V(NRes + k) = - E(nodo2); end;
        end;
    end;
    // Inizializza I
    I = zeros(R,1);
    // Costruisce I analizzando Res e Gen una riga alla volta
    for k = 1:NRes,
        riga = Res(k,:);
        valore = riga(1);
        I(k) = V(k) / valore;
    end;
    for k = 1:NGen,
        riga = Gen(k,:);
        valore = riga(1);
        I(NRes + k) = valore;
    end;
endfunction
