//
// Esercitazione 7
//
// ***************************
//
//     Procedure 
//        *  SI
//        *  SA
//
// ***************************
//
function x = SI(T,c)
//
// Per T matrice N x N triangolare superiore invertibile e c colonna di N
// componenti, calcola un'approssimazione x della colonna soluzione del sistema
// di equazioni lineari T z = c.
//
N = size(T,'r');
//
// *** Verifica se T invertibile. 
//
if find(diag(T)==0) ~= [] then
   error("elementi nulli sulla diagonale");
   else
       //
       // *** Calcolo soluzione
       //
       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;
end;
endfunction
//
function x = SA(T,c)
//
// Per T matrice N x N triangolare inferiore invertibile e c colonna di N
// componenti, calcola un'approssimazione x della colonna soluzione del sistema
// di equazioni lineari T z = c
//
N = size(T, 'r');
//
// *** Verifica se T invertibile. 
//
if find(diag(T)==0) ~= [] then
   error("elementi nulli sulla diagonale");
   else
       //
       // *** Calcolo soluzione
       //
       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;
end
endfunction
//
// *********************************************************
//
//     Procedure:
//        *  EquazioniNodaliFormaSimbolica
//        *  SoluzioneReteMetodoFasoriale
//
// *********************************************************
//
function [Y,J] = EquazioniNodaliFormaSimbolica(Rete)
    //
    // Determina la matrice delle ammettenze Y ed il termine noto J del
    // sistema delle equazioni nodali per il circuito di sole impedenze
    // 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 NZ x 3 ad elementi complessi con righe, una per ciascuna
    //          impedenza del circuito, di formato
    //               valore impedenza , nodo 1 , nodo 2
    // Rete(2): matrice NGen x 3 ad elementi complessi con righe, una per 
    //          ciascun generatore di corrente del circuito, di formato
    //               valore generatore , nodo 1 , nodo 2
    //
    // Z: matrice N x N ad elementi complessi.
    // J: colonna di N numeri complessi.
    //
    // Determina il numero di nodi della rete, di impedenze e di generatori
    [Z, Gen] = Rete(:);
    Nodi = union(Z(:,2:3), Gen(:,2:3));
    N = max(Nodi);
    NZ = size(Z,'r');
    NGen = size(Gen,'r')
    // Inizializza Y e J
    Y = zeros(N,N);
    J = zeros(N,1);
    // Costruisce Y analizzando Z una riga alla volta
    for k = 1:NZ,
        riga = Z(k,:);
        valore = riga(1); nodo1 = riga(2); nodo2 = riga(3);
        if nodo1 ~= 0 & nodo2 ~= 0 then
           Y([nodo1,nodo2],[nodo1,nodo2]) = ...
           Y([nodo1,nodo2],[nodo1,nodo2]) + ...
               [1,-1 ;...
               -1, 1] / valore;
           else if nodo1 ~= 0 then
                   Y(nodo1,nodo1) = Y(nodo1,nodo1) + 1 / valore; end;
                if nodo2 ~= 0 then
                   Y(nodo2,nodo2) = Y(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] = SoluzioneReteMetodoFasoriale(Rete)
    //
    // Determina tutti i fasori delle 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 NZ x 3 ad elementi complessi con righe, una per ciascuna
    //          impedenza del circuito, di formato
    //               valore impedenza , 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
    //
    // I rami del circuito sono numerati considerando, in ordine, per primi gli
    // NZ rami contenenti impedenze poi gli NGEN rami contenenti generatori di
    // corrente.
    //
    // E: colonna di N numeri complessi. E' il vettore dei fasori delle tensioni di nodo.
    // V: colonna di R numeri complessi. E' il vettore dei fasori delle tensioni di ramo.
    // I: colonna di R numeri complessi. E' il vettore dei fasori delle correnti di ramo.
    //
    // Determina i dati delle equazioni nodali
    [Y, J] = EquazioniNodaliFormaSimbolica(Rete);
    // (1) fattorizzazione di Y
    [S,D,P] = lu(Y);
    // (2) soluzione del sistema S x = PJ
    c = SA(S, P*J);
    // (3) soluzione del sistema D x = c
    E = SI(D,c);
    // Determina le grandezze di ramo del circuito.
    [Z, Gen] = Rete(:);
    Nodi = union(Z(:,2:3), Gen(:,2:3));
    N = max(Nodi);
    NZ = size(Z,'r');
    NGen = size(Gen,'r')
    NRami = NZ + NGen;
    // Inizializza V
    V = zeros(NRami,1);
    // Costruisce V analizzando Z e Gen una riga alla volta
    for k = 1:NZ,
        riga = Z(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(NZ + k) = E(nodo1) - E(nodo2);
           else if nodo1 ~= 0 then V(NZ + k) = E(nodo1); end;
                if nodo2 ~= 0 then V(NZ + k) = - E(nodo2); end;
        end;
    end;
    // Inizializza I
    I = zeros(NRami,1);
    // Costruisce I analizzando Z e Gen una riga alla volta
    for k = 1:NZ,
        riga = Z(k,:);
        valore = riga(1);
        I(k) = V(k) / valore;
    end;
    for k = 1:NGen,
        riga = Gen(k,:);
        valore = riga(1);
        I(NZ + k) = valore;
    end;
endfunction
//
// *********************************************************
//
//     Procedure:
//       *  Fasore
//       *  FormaPolare
//       *  Grafico
//
// *********************************************************
//
function [] = Fasore(z)
    m = max(abs(z));
    plot([zeros(1,length(z));real(z)'],[zeros(1,length(z));imag(z)'],'thickness',2);
    plot2d(0,rect=[-m,-m,m,m],frameflag = 3, strf='034'); xgrid();
    legenda = [];
    for i = 1:length(z),
        legenda(i) = '('+string(i)+')';
    end;
    legend(legenda, 'out_upper_right');
endfunction
//
function [Modulo, Fase] = FormaPolare(z)
    // z: vettore a componenti numeri complessi in forma cartesiana;
    // modulo(k): |z(k)|
    // fase(k): argomento di z(k) in (-pi, pi].
    Modulo = abs(z);
    Fase = atan(imag(z), real(z)) * 180/%pi;  // in gradi
    for j = 1:length(z),
        mprintf('(%d)   %.2g < %.2f\n', j, Modulo(j), Fase(j));
    end
endfunction
//
function [] = Grafico(z)
    // Disegna, su uno stesso piano cartesiano, un periodo delle 
    // funzioni
    //           Re( z exp(i omega t) )
    Modulo = abs(z);
    Argomento = atan(imag(z), real(z));
    x = linspace(0,2*%pi,300)';
    y = zeros(300,length(z));
    for i = 1:length(z),
        y(:,i) = Modulo(i) * cos(x + Argomento(i));
    end
    plot(x, y,'thickness',2); xgrid();
    Assi = gca();
    Assi.tight_limits(1) = 'on'
    for i = 1:length(z),
        legenda(i) = '('+string(i)+')';
    end;
    legend(legenda, 'lower_caption');
    xlabel('w t');
endfunction
