%--------------------------------------------------------------------------
%      SIMULAZIONE DELLA DIFFUSIONE DEL CALORE ALL'INTERNO DI UN
%              PARALLELEPIPEDO DI MATERIALE NON OMOGENEO
%--------------------------------------------------------------------------
%
%                      Parametri fisici del sistema
%        __________       
%       |\ ________|\      -> L'oggetto ha una capacita' termica variabile c(x,y,z)
%       | |        | |     -> L'oggetto ha una densità variabile r(x,y,z)
%  (z)h | |        | |     -> L'oggetto ha una conducibilita' termica variabile k(x,y,z) 
%       | |        | |     -> Gli scambi di calore per convezione con l'esterno avvengono
%    (x)w\|_________\|        secondo la legge q(x,y,x)=h(x,y,z)*(T-T_ext)
%              l              dove h e' un parametro che dipende dal
%             (y)             punto sulla superficie e T_ext e' la
%                             temperatura esterna che puo' variare da punto
%                             a punto e nel tempo.
%                          -> Dal sole arriva anche un flusso di calore
%                             q(x,y,z) [W/m^2] assorbito dall'edificio
%                       
%
% L'oggetto e' cosi' suddiviso:
% Si crea una griglia di punti regolare, lungo la direzione h ci sono nh
% punti, lungo w ce ne sono nw e lungo l nl. In tutto il numero di nodi N
% e' dato da N=nh*nl*nw

clear all;
close all;
clc;
warning off;
Parametri_Globali;

SalvaSimulazione=1;
FileSalvataggio='DatiSim.mat';
CaricaSimulazione=0;
FileCaricamento='DatiSim.mat';
GraficoSezioneTemp=1;
TempoSezione=20;
PassoSezione=1;
hSezione1=-0.25;
hSezione2=0;
hSezione3=0.25;

if (CaricaSimulazione==0)

    nh=3;
    nw=3;
    nl=3;
    N=nh*nw*nl;

    GraficoMesh=1;
    Passi=2000;  % Passi della simulazione
    DeltaT=15;   % Intervallo campionamento simulazione
    TempoSimComplessivo=0;
    T0=273;
    Temp_ini=T0*ones(N,1);   
       
    %generiamo i vertici con le loro coordinate... il baricentro del parallelepipedo
    %e' nell'origine
    Inizio = clock;
    fprintf('Generazione mesh...');
    
    DensitaMuro=5000;
    CapacitaMuro=440;
    ConducibilitaMuro=7730;
    
    %Parete di sinistra
%     [x,y,z] = ndgrid(-1:0.05:-0.8,-(2*nl-1)/2:(2*nl-1)/2,-(2*nh-1)/2:(2*nh-1)/2);
%     y=y.*(L/(2*nl-1)); z=z.*(H/(2*nh-1));
%     x=x(:);
%     y=y(:);
%     z=z(:);
%     
%     N=length(x);
    
%     %  DUE TETRAEDRI
%     x=[1; 0; 1; 1; 2];
%     y=[0; 0; 1; 0; 0];
%     z=[1; 0; 0; 0; 0];
%     
%     N=5;

%     % UN CUBO
%     x=[0; 0; 0; 0; 1; 1; 1; 1;];
%     y=[0; 0; 1; 1; 0; 0; 1; 1;];
%     z=[0; 1; 0; 1; 0; 1; 0; 1;];
%     
%     N=8;

    % DUE CUBI
%     x=[0; 0; 0; 0; 1; 1; 1; 1; 2; 2; 2; 2];
%     y=[0; 0; 1; 1; 0; 0; 1; 1; 0; 0; 1; 1];
%     z=[0; 1; 0; 1; 0; 1; 0; 1; 0; 1; 0; 1];
%     
%    N=12;

    % QUATTRO CUBI
%     x=[0; 0; 0; 0; 1; 1; 1; 1; 2; 2; 2; 2; 0; 0; 1; 1; 2; 2];
%     y=[0; 0; 1; 1; 0; 0; 1; 1; 0; 0; 1; 1; 1; 0; 1; 0; 1; 0];
%     z=[0; 1; 0; 1; 0; 1; 0; 1; 0; 1; 0; 1; 2; 2; 2; 2; 2; 2];
%     x=x.*(W/(nw-1)); y=y.*(L/(nl-1)); z=z.*(H/(nh-1));
%     x=x(:);
%     y=y(:);
%     z=z(:);

    % OTTO CUBI
    x=[0; 0; 0; 0; 1; 1; 1; 1; 2; 2; 2; 2; 0; 0; 1; 1; 2; 2; 0; 0; 0; 1; 1; 1; 2; 2; 2];
    y=[0; 0; 1; 1; 0; 0; 1; 1; 0; 0; 1; 1; 1; 0; 1; 0; 1; 0; -1;-1;-1;-1;-1;-1;-1;-1;-1];
    z=[0; 1; 0; 1; 0; 1; 0; 1; 0; 1; 0; 1; 2; 2; 2; 2; 2; 2; 0; 1; 2; 0; 1; 2; 0; 1; 2];
    x=x.*(W/(nw-1)); y=y.*(L/(nl-1)); z=z.*(H/(nh-1));
    x=x(:);
    y=y(:);
    z=z(:);
    
%     [x,y,z] = meshgrid(-(nw-1)/2:(nw-1)/2,-(nl-1)/2:(nl-1)/2,-(nh-1)/2:(nh-1)/2);
%     x=x.*(W/(nw-1)); y=y.*(L/(nl-1)); z=z.*(H/(nh-1));
%     x=x(:);
%     y=y(:);
%     z=z(:);
    
    IndiciInterni=intersect(intersect(find(x>-W/2+TOLLERANZA & x<=W/2-TOLLERANZA),find(y>-L/2+TOLLERANZA & y<L/2-TOLLERANZA)),find(z>-H/2+TOLLERANZA & z<H/2-TOLLERANZA));
    
    N=length(x);
    
    x_perturb=x;
%    x_perturb(IndiciInterni)=x_perturb(IndiciInterni)+0.1*randn(length(IndiciInterni),1);
    y_perturb=y;
%    y_perturb(IndiciInterni)=y_perturb(IndiciInterni)+0.1*randn(length(IndiciInterni),1);
    z_perturb=z;
%    z_perturb(IndiciInterni)=z_perturb(IndiciInterni)+0.1*randn(length(IndiciInterni),1);
    
    T = uint32(delaunay3(x_perturb, y_perturb, z_perturb, {'Qt', 'Qbb', 'Qc', 'Qz'}));
    
    [Ne App]=size(T);  %Ne contiene il numero complessivo di elementi
    PT = zeros(Ne,4);
    ControllaVolumi;
    PT(:,2)=DensitaMuro*ones(Ne,1);
    PT(:,3)=CapacitaMuro*ones(Ne,1);
    PT(:,4)=ConducibilitaMuro*ones(Ne,1);
    GeneraFacce;
    
    for i=1:N
        if (x(i)^2+y(i)^2+z(i)^2<0.2^2)
            Temp_ini(i)=400;
        end
    end  
    
    % Cosi' come le equazioni differenziali ordinarie hanno bisogno di una
    % discretizzazione temporale per poter essere simulate, le equazioni alle
    % derivate parziali hanno bisogno di una discretizzazione sia spaziale che
    % temporale che sono considerate in modo diverso. La discretizzazione
    % spaziale si ottiene suddividendo l'oggetto in parti piu' piccole (gli
    % elementi) che, nel caso piu' semplice in 3 dimensioni, sono tetraedri.
    % Quindi bisogna unire i vertici appena ottenuti 4 a 4 in modo da
    % ottenere una copertura dell'intero solido costituita da elementi. 
    % Questa copertura utilizzando tetraedri puo' sempre essere fatta ma non e'
    % univoca... l'algoritmo utilizzato e' quello di delaunay

    %T = delaunay3(x,y,z);
    %[Ne App]=size(T);  %Ne contiene il numero complessivo di elementi
    TempoSimComplessivo=TempoSimComplessivo+etime(clock,Inizio);
    fprintf(strcat(strcat('fatto ( ',num2str(etime(clock,Inizio))),'sec)\n'));
    if (GraficoMesh)
        handle=figure;
        plot3(x(1),y(1),z(1),'r');
        hold;
        for i=2:N
            plot3(x(i),y(i),z(i),'r');
        end
        for i=1:Ne
            plot3([x(T(i,1)); x(T(i,2))],[y(T(i,1)); y(T(i,2))],[z(T(i,1)); z(T(i,2))],'r');
            plot3([x(T(i,1)); x(T(i,3))],[y(T(i,1)); y(T(i,3))],[z(T(i,1)); z(T(i,3))],'r');
            plot3([x(T(i,1)); x(T(i,4))],[y(T(i,1)); y(T(i,4))],[z(T(i,1)); z(T(i,4))],'r');
            plot3([x(T(i,2)); x(T(i,3))],[y(T(i,2)); y(T(i,3))],[z(T(i,2)); z(T(i,3))],'r');
            plot3([x(T(i,2)); x(T(i,4))],[y(T(i,2)); y(T(i,4))],[z(T(i,2)); z(T(i,4))],'r');
            plot3([x(T(i,3)); x(T(i,4))],[y(T(i,3)); y(T(i,4))],[z(T(i,3)); z(T(i,4))],'r');
        end
    end

    return;
    
    % L'equazione che prenderemo in esame sarà
    %       _   _        _   _         _   _
    %   d  |  dT |   d  |  dT |    d  |  dT |      dT
    %   -- |k -- | + -- |k -- | +  -- |k -- | - rc -- = 0
    %   dx |_ dx_|   dy |_ dy_|    dz |_ dz_|      dt
    %
    % La soluzione di questa equazione viene rappresentata mediante
    % una funzione che è lineare nelle coordinate spaziali all'interno di ogni
    % elemento. Quindi all'interno di un tetraedro la temperatura sarà
    % T=a+bx+cy+dz. Se Ta Tb Tc e Td sono (in un certo istante t) le temperature dei
    % vertici di tale tetraedro, allora la temperatura al suo interno risulta determinata
    % per continuità a partire solo dalle temperature dei quattro nodi. Basta
    % infatti imporre che la temperatura valga Ta Tb Tc e Td nei quattro
    % vertici e dalle 4 condizioni che ne risultano si ottengono le espressioni
    % per a, b, c, d. Una volta sostituite queste espressioni in quella per T
    % ed aver messo in evidenza la temperatura dei 4 nodi si perviene a
    %
    %      T=Na(x,y,z)Ta+Nb(x,y,z)Tb+Nc(x,y,z)Tc+Nd(x,y,z)Td
    % 
    % Le Na, Nb, Nc ed Nd sono dette funzioni di forma e dipendono, oltre che
    % dal punto (x,y,z) interno all'elemento in cui si vuole calcolare la
    % tempaeratura, anche dalle coordinate (xa,ya,za), ..., (xd,yd,zd) dei 4
    % vertici dell'elemento. Il bello delle funzioni di forma è che esistono
    % delle simpatiche formulette (valgono solo nell'ipotesi di temperatura 
    % lineare all'interno dell'elemento):
    %
    %VOLUME DELL'ELEMENTO  
    %
    %     1    | 1 xa ya za |
    % V = - det| 1 xb yb zb |
    %     6    | 1 xc yc zc |
    %          | 1 xd yd zd |
    %
    %AREA DELLA FACCIA (ad esempio tra i vertici a, b, c)
    %          _________________________________________________
    %     1   |    |ya za 1|^2      |za xa 1|^2      |xa ya 1|^2
    % A = - \ | det|yb zb 1|   + det|zb xb 1|   + det|xb yb 1|
    %     2  \|    |yc zb 1|        |zc xc 1|        |xc yc 1|
    %
    %INTEGRAZIONE IN TUTTO L'ELEMENTO DELLE FUNZIONI FORMA
    %    _
    %   |                                 p!q!r!s!6V
    %   | Na^p * Nb^q * Nc^r * Nd^s dV = ------------
    %  _|                                (p+q+r+s+3)!
    %    Tetraedro
    %
    %INTEGRAZIONE SUPERFICIALE DELLE FUNZIONI FORMA (faccia a,b,c)
    %    _
    %   |                          p!q!r!2A
    %   | Na^p * Nb^q * Nc^r dS = ----------
    %  _|                         (p+q+r+2)!
    %    Faccia
    pack;
    K=zeros(N,N);
    C=zeros(N,N);
    F=zeros(Passi,N);
       
    % -------------------------------------------------------------------------
    %                                                 dT  
    %        Formulazione delle equazioni del tipo  C -- = K T + F
    %                                                 dt
    % -------------------------------------------------------------------------
    Inizio = clock;
    fprintf('Calcolo matrici del modello...');
    for i=1:Ne %imponiamo l'annullamento su ogni elemento
        T1=T(i,1); T2=T(i,2); T3=T(i,3); T4=T(i,4);
        %Calcolo del volume dell'elemento, del suo baricentro, e dei baricentri
        %delle quattro facce
        Baricentro=[ 0.25*(x(T1)+x(T2)+x(T3)+x(T4));
                     0.25*(y(T1)+y(T2)+y(T3)+y(T4));
                     0.25*(z(T1)+z(T2)+z(T3)+z(T4)); ];
        SeiV=det([1 x(T1) y(T1) z(T1);
                  1 x(T2) y(T2) z(T2);
                  1 x(T3) y(T3) z(T3);
                  1 x(T4) y(T4) z(T4)]);
        if (SeiV<0)      
            SeiV=-SeiV;
            App=T1;
            T1=T2;
            T2=App;
        end
        BaricentroF1=[(x(T1)+x(T2)+x(T3))/3;
                      (y(T1)+y(T2)+y(T3))/3;
                      (z(T1)+z(T2)+z(T3))/3; ];
        BaricentroF2=[(x(T1)+x(T2)+x(T4))/3;
                      (y(T1)+y(T2)+y(T4))/3;
                      (z(T1)+z(T2)+z(T4))/3; ];
        BaricentroF3=[(x(T1)+x(T4)+x(T3))/3;
                      (y(T1)+y(T4)+y(T3))/3;
                      (z(T1)+z(T4)+z(T3))/3; ];
        BaricentroF4=[(x(T4)+x(T2)+x(T3))/3;
                      (y(T4)+y(T2)+y(T3))/3;
                      (z(T4)+z(T2)+z(T3))/3; ];

        %Determinazione di quali (eventuali) facce dell'elemento fanno parte
        %della superficie dell'edificio e calcolo della loro area
        [Faccia1Esterna ConvCoef1]=FacciaEsterna(BaricentroF1);
        if (Faccia1Esterna)
            DueA1=sqrt(det([y(T1) z(T1) 1; y(T2) z(T2) 1; y(T3) z(T3) 1])^2+det([z(T1) x(T1) 1; z(T2) x(T2) 1; z(T3) x(T3) 1])^2+det([x(T1) y(T1) 1; x(T2) y(T2) 1; x(T3) y(T3) 1])^2);
        end
        [Faccia2Esterna ConvCoef2]=FacciaEsterna(BaricentroF2);
        if (Faccia2Esterna)
            DueA2=sqrt(det([y(T1) z(T1) 1; y(T2) z(T2) 1; y(T4) z(T4) 1])^2+det([z(T1) x(T1) 1; z(T2) x(T2) 1; z(T4) x(T4) 1])^2+det([x(T1) y(T1) 1; x(T2) y(T2) 1; x(T4) y(T4) 1])^2);
        end
        [Faccia3Esterna ConvCoef3]=FacciaEsterna(BaricentroF3);
        if (Faccia3Esterna)
            DueA3=sqrt(det([y(T1) z(T1) 1; y(T3) z(T3) 1; y(T4) z(T4) 1])^2+det([z(T1) x(T1) 1; z(T3) x(T3) 1; z(T4) x(T4) 1])^2+det([x(T1) y(T1) 1; x(T3) y(T3) 1; x(T4) y(T4) 1])^2);
        end
        [Faccia4Esterna ConvCoef4]=FacciaEsterna(BaricentroF4);
        if (Faccia4Esterna)
            DueA4=sqrt(det([y(T2) z(T2) 1; y(T3) z(T3) 1; y(T4) z(T4) 1])^2+det([z(T2) x(T2) 1; z(T3) x(T3) 1; z(T4) x(T4) 1])^2+det([x(T2) y(T2) 1; x(T3) y(T3) 1; x(T4) y(T4) 1])^2);
        end          

        %--------- CALCOLO COEFFICIENTI FUNZIONI DI FORMA------------------
        %Calcolo esplicito coefficienti Na
        b1=(y(T2)*(z(T4)-z(T3))-y(T3)*(z(T4)-z(T2))+y(T4)*(z(T3)-z(T2)))/SeiV;
        c1=(-x(T2)*(z(T4)-z(T3))+x(T3)*(z(T4)-z(T2))-x(T4)*(z(T3)-z(T2)))/SeiV;
        d1=(x(T2)*(y(T4)-y(T3))-x(T3)*(y(T4)-y(T2))+x(T4)*(y(T3)-y(T2)))/SeiV;
        
        %Calcolo esplicito coefficienti Nb
        b2=(-y(T1)*(z(T4)-z(T3))+y(T3)*(z(T4)-z(T1))-y(T4)*(z(T3)-z(T1)))/SeiV;
        c2=(x(T1)*(z(T4)-z(T3))-x(T3)*(z(T4)-z(T1))+x(T4)*(z(T3)-z(T1)))/SeiV;
        d2=(-x(T1)*(y(T4)-y(T3))+x(T3)*(y(T4)-y(T1))-x(T4)*(y(T3)-y(T1)))/SeiV;
        
        %Calcolo esplicito coefficienti Nc
        b3=(y(T1)*(z(T4)-z(T2))-y(T2)*(z(T4)-z(T1))+y(T4)*(z(T2)-z(T1)))/SeiV;
        c3=(-x(T1)*(z(T4)-z(T2))+x(T2)*(z(T4)-z(T1))-x(T4)*(z(T2)-z(T1)))/SeiV;
        d3=(x(T1)*(y(T4)-y(T2))-x(T2)*(y(T4)-y(T1))+x(T4)*(y(T2)-y(T1)))/SeiV;
        
        %Calcolo esplicito coefficienti Nd
        b4=(-y(T1)*(z(T3)-z(T2))+y(T2)*(z(T3)-z(T1))-y(T3)*(z(T2)-z(T1)))/SeiV;
        c4=(x(T1)*(z(T3)-z(T2))-x(T2)*(z(T3)-z(T1))+x(T3)*(z(T2)-z(T1)))/SeiV;
        d4=(-x(T1)*(y(T3)-y(T2))+x(T2)*(y(T3)-y(T1))-x(T3)*(y(T2)-y(T1)))/SeiV;
        %------------------------------------------------------------------
        
        %------------- CALCOLO PROPRIETA' FISICHE DELL'ELEMENTO -----------
        [Dens Cap Cond]= ProprietaMezzo(Baricentro);
        

        %-------------------- CALCOLO MATRICE C----------------------------
        %Integrale di c*r*Na^2 all'interno del tetraedro
        I1=Cap*Dens*SeiV/60; %Formula semplificata... considero c ed r costati all'interno del tetraedro
                                                                  %e pari al valore assunto nel baricentro... da migliorare! 
                                                                  %la cosa bella è che, usando la formula semplificata, l'integrale è
                                                                  %uguale se uso Na Nb o Nc
        %Integrale di c*r*Na*Nb all'interno del tetraedro
        I2=Cap*Dens*SeiV/120;%Formula semplificata... considero c ed r costanti all'interno del tetraedro
                                                                  %e pari al valore assunto nel baricentro... da migliorare! 
                                                                  %la cosa bella è che, usando la formula semplificata, l'integrale è
                                                                  %indipendente dalle funzioni di forma coinvolte
        C(T1,T1)=C(T1,T1)+I1;
        C(T2,T2)=C(T2,T2)+I1;
        C(T3,T3)=C(T3,T3)+I1;
        C(T4,T4)=C(T4,T4)+I1;
        C(T1,T2)=C(T1,T2)+I2;
        C(T2,T1)=C(T2,T1)+I2;
        C(T1,T3)=C(T1,T3)+I2;
        C(T3,T1)=C(T3,T1)+I2;
        C(T1,T4)=C(T1,T4)+I2;
        C(T4,T1)=C(T4,T1)+I2;
        C(T2,T3)=C(T2,T3)+I2;
        C(T3,T2)=C(T3,T2)+I2;
        C(T2,T4)=C(T2,T4)+I2;
        C(T4,T2)=C(T4,T2)+I2;
        C(T3,T4)=C(T3,T4)+I2;
        C(T4,T3)=C(T4,T3)+I2;
        %------------------------------------------------------------------

        %-------------------- CALCOLO MATRICE K----------------------------    
        %Integrale di k(x,y,z) all'interno dell'elemento
        I3=Cond*SeiV/6; %Formula semplificata... considero k costante all'interno dell'elemento... da migliorare!
        
        %Integrali superficiali
        %I4 = Integrale sulla superficie dell'elemento (solo della parte che è superficie dell'edificio) di h*Ni^2
        %I5 = Integrale sulla superficie dell'elemento (solo della parte che è superficie dell'edificio) di h*Ni*Nj
        I4_1=0; I4_2=0; I4_3=0; I4_4=0;
        I5_1=0; I5_2=0; I5_3=0; I5_4=0;
        %if (i==10)
        %    ciccio=0;
        %end
        if (Faccia1Esterna)
            I4_1=ConvCoef1*DueA1/12; %Formule semplificate... considero h costante sulla faccia... da migliorare!
            I5_1=ConvCoef1*DueA1/24; % tra l'altro in questo modo sono indipendenti dalle funzioni di forma coinvolte
        end
        if (Faccia2Esterna)
            I4_2=ConvCoef2*DueA2/12; %Formule semplificate... considero h costante sulla faccia... da migliorare!
            I5_2=ConvCoef2*DueA2/24; % tra l'altro in questo modo sono indipendenti dalle funzioni di forma coinvolte
        end
        if (Faccia3Esterna)
            I4_3=ConvCoef3*DueA3/12; %Formule semplificate... considero h costante sulla faccia... da migliorare!
            I5_3=ConvCoef3*DueA3/24; % tra l'altro in questo modo sono indipendenti dalle funzioni di forma coinvolte
        end
        if (Faccia4Esterna)
            I4_4=ConvCoef4*DueA4/12; %Formule semplificate... considero h costante sulla faccia... da migliorare!
            I5_4=ConvCoef4*DueA4/24; % tra l'altro in questo modo sono indipendenti dalle funzioni di forma coinvolte
        end
        
        %        --->>  PARTE TEORICA BASTARDISSIMA!!!!! <<---
        % Se integro Na sulla superficie opposta al vertice a trovo 0
        % perché su quella faccia Na è nullo. Lo stesso vale per le altre
        % funzioni di forma. Analogamente il prodotto Na*Nb è nullo sulle
        % facce che non contengono i nodi a e b!! quindi ci vuole OCCHIO A
        % NON SOMMARLI TUTTI MA SOLO QUELLI NON NULLI!!!        
        
        %Coefficienti dal nodo i al nodo i
        K(T1,T1)=K(T1,T1)-(I4_1+I4_2+I4_3)-I3*(b1^2+c1^2+d1^2);
        K(T2,T2)=K(T2,T2)-(I4_1+I4_2+I4_4)-I3*(b2^2+c2^2+d2^2);
        K(T3,T3)=K(T3,T3)-(I4_1+I4_3+I4_4)-I3*(b3^2+c3^2+d3^2);
        K(T4,T4)=K(T4,T4)-(I4_2+I4_3+I4_4)-I3*(b4^2+c4^2+d4^2);
        %Coefficienti dal nodo i al nodo j
        K(T1,T2)=K(T1,T2)-(I5_1+I5_2)-I3*(b1*b2+c1*c2+d1*d2);
        K(T1,T3)=K(T1,T3)-(I5_1+I5_3)-I3*(b1*b3+c1*c3+d1*d3);
        K(T1,T4)=K(T1,T4)-(I5_2+I5_3)-I3*(b1*b4+c1*c4+d1*d4);
        K(T2,T1)=K(T2,T1)-(I5_1+I5_2)-I3*(b1*b2+c1*c2+d1*d2);
        K(T2,T3)=K(T2,T3)-(I5_1+I5_4)-I3*(b2*b3+c2*c3+d2*d3);
        K(T2,T4)=K(T2,T4)-(I5_2+I5_4)-I3*(b2*b4+c2*c4+d2*d4);
        K(T3,T1)=K(T3,T1)-(I5_1+I5_3)-I3*(b1*b3+c1*c3+d1*d3);
        K(T3,T4)=K(T3,T4)-(I5_3+I5_4)-I3*(b3*b4+c3*c4+d3*d4);        
        K(T3,T2)=K(T3,T2)-(I5_1+I5_4)-I3*(b2*b3+c2*c3+d2*d3);
        K(T4,T2)=K(T4,T2)-(I5_2+I5_4)-I3*(b2*b4+c2*c4+d2*d4);
        K(T4,T1)=K(T4,T1)-(I5_2+I5_3)-I3*(b1*b4+c1*c4+d1*d4);
        K(T4,T3)=K(T4,T3)-(I5_3+I5_4)-I3*(b3*b4+c3*c4+d3*d4);
        %------------------------------------------------------------------

        %----------------- CALCOLO INGRESSI (MATRICE F) ------------------- 
        if (Faccia1Esterna || Faccia2Esterna || Faccia3Esterna || Faccia4Esterna)
            I6_1=0; I6_2=0; I6_3=0; I6_4=0;
            for j=1:Passi
                % clacolo l'Integrale superficiale di Na*(-q+h*Ta)... uso il solito trucco
                %di considerare q, h e Ta costanti sulle facce                
                if (Faccia1Esterna)
                    [Text Irr]=CondizioniContorno(BaricentroF1,(j-1)*DeltaT);
                    I6_1=(DueA1/6)*(ConvCoef1*Text+Irr);
                end
                if (Faccia2Esterna)
                    [Text Irr]=CondizioniContorno(BaricentroF2,(j-1)*DeltaT);
                    I6_2=(DueA2/6)*(ConvCoef2*Text+Irr);
                end
                if (Faccia3Esterna)
                    [Text Irr]=CondizioniContorno(BaricentroF3,(j-1)*DeltaT);
                    I6_3=(DueA3/6)*(ConvCoef3*Text+Irr);
                end
                if (Faccia4Esterna)
                    [Text Irr]=CondizioniContorno(BaricentroF4,(j-1)*DeltaT);
                    I6_4=(DueA4/6)*(ConvCoef4*Text+Irr);
                end
                F(j,T1)=F(j,T1)+I6_1+I6_2+I6_3;
                F(j,T2)=F(j,T2)+I6_1+I6_2+I6_4;
                F(j,T3)=F(j,T3)+I6_1+I6_3+I6_4;
                F(j,T4)=F(j,T4)+I6_2+I6_3+I6_4;
            end
        end
        %------------------------------------------------------------------
    end
    TempoSimComplessivo=TempoSimComplessivo+etime(clock,Inizio);
    fprintf(strcat(strcat('fatto ( ',num2str(etime(clock,Inizio))),'sec)\n'));
    
    pack;
    %------- CALCOLO DELLE MATRICI DEL SISTEMA DA SIMULARE------------
    Inizio = clock;
    fprintf('Costruzione del sistema dinamico a tempo discreto...');
    B=inv(C);
    [A, B]=c2dm(B*K,B,eye(N),zeros(N,N),DeltaT,'zoh');
    TempoSimComplessivo=TempoSimComplessivo+etime(clock,Inizio);
    fprintf(strcat(strcat('fatto ( ',num2str(etime(clock,Inizio))),'sec)\n'));

    pack;
    %---------------------- SIMULAZIONE SISTEMA-----------------------
    Inizio = clock;
    fprintf('Simulazione della risposta...');
    t=DeltaT*(0:(Passi-1));
    Temp=Temp_ini';
    for i=2:Passi
        TempAtt=A*Temp(i-1,:)'+B*F(i-1,:)';
        Temp=[Temp; TempAtt'];
    end
    TempoSimComplessivo=TempoSimComplessivo+etime(clock,Inizio);
    fprintf(strcat(strcat('fatto ( ',num2str(etime(clock,Inizio))),'sec)\n'));

    %--------------- SALVATAGGIO DATI SIMULAZIONE --------------------
    if (SalvaSimulazione)
        save(FileSalvataggio, 'Temp', 't', 'x', 'y', 'z', 'H', 'nh', 'L', 'nl', 'W', 'nw', 'N', 'Passi', 'DeltaT', 'A', 'B', 'F','TempoSimComplessivo');
    end

else
    load(FileCaricamento);
end


%------------- VISUALIZZAZIONE TEMPERATURA NELLE SEZIONI -------------
if(GraficoSezioneTemp)
    handle=figure;
    v=get(0,'ScreenSize');
    set(handle,'Position',[round(v(3)/2)-700 round(v(4)/2)-200 1400 400]);
    VettAltezze=(((-(nh-1)/2):((nh-1)/2) )*(H/(nh-1)));
    [m,i]=min(abs((-(nh-1)/2:(nh-1)/2)*(H/(nh-1))-ones(1,nh)*hSezione1));
    Altezza1=VettAltezze(i);
    IndiciBuoni1=find(abs(z-Altezza1*ones(N,1))<TOLLERANZA);
    [m,i]=min(abs((-(nh-1)/2:(nh-1)/2)*(H/(nh-1))-ones(1,nh)*hSezione2));
    Altezza2=VettAltezze(i);
    IndiciBuoni2=find(abs(z-Altezza2*ones(N,1))<TOLLERANZA);
    [m,i]=min(abs((-(nh-1)/2:(nh-1)/2)*(H/(nh-1))-ones(1,nh)*hSezione3));
    Altezza3=VettAltezze(i);
    IndiciBuoni3=find(abs(z-Altezza3*ones(N,1))<TOLLERANZA);
    for j=1:(round(Passi/(TempoSezione/PassoSezione))):Passi
        AppSez1=Temp(j,IndiciBuoni1)';
        AppSez2=Temp(j,IndiciBuoni2)';
        AppSez3=Temp(j,IndiciBuoni3)';
        App=AppSez1(1:nw)';
        for i=2:nl
            App=[App; AppSez1((1+(i-1)*nw):(i*nw))'];
        end
        subplot(1,3,1);
        meshc((-(nw-1)/2:(nw-1)/2)*(W/(nw-1)), (-(nl-1)/2:(nl-1)/2)*(L/(nl-1)), App);
        title(strcat('Altezza sezione = ',num2str(Altezza1)));
        App=AppSez2(1:nw)';
        for i=2:nl
            App=[App; AppSez2((1+(i-1)*nw):(i*nw))'];
        end
        subplot(1,3,2);
        meshc((-(nw-1)/2:(nw-1)/2)*(W/(nw-1)), (-(nl-1)/2:(nl-1)/2)*(L/(nl-1)), App);
        title(strcat('Altezza sezione = ',num2str(Altezza2)));
        App=AppSez3(1:nw)';
        for i=2:nl
            App=[App; AppSez3((1+(i-1)*nw):(i*nw))'];
        end
        subplot(1,3,3);
        meshc((-(nw-1)/2:(nw-1)/2)*(W/(nw-1)), (-(nl-1)/2:(nl-1)/2)*(L/(nl-1)), App);
        title(strcat('Altezza sezione = ',num2str(Altezza3)));
        pause(PassoSezione);
    end
end













