//IRAKASGAIA: ESPERIMENTAZIOA INGENIERITZA KIMIKOAN I IRAKASGAIKO 
// OCW-EHU
// http://ocw.ehu.es/irakaskuntza-teknikoak/zenbakizko-kalkulua-ingeniari-kimikoentzat/Course_listing
// LANA 2: BERO TRANSMISIOA BI NORABIDETAN
//ARIKETA MOTA: EKUAZIO DIFERENTZIAL PARTZIALAK -LAPLACE MOTAKOA- 
//ALGORITMOA: DIFERENTZIA FINITOAK ETA GAUSS-SEIDEL
//EGILEA: ASIER ARANZABAL
//DATA: 2009/7/9


 
function [T,errorea,iter,a,l,dx]=tenperatura_banaketa4(Tezker,Teskuin,Tgoi,Tbehe,A,L) 

//NODOAK SORTZEN.....HASIERAKO KALKULUA (5 X 5 MATRIZEA)
tarte=4
tol=0.1
dx=L/tarte;
dy=A/tarte;
l=tarte+1;
a=tarte+1;
T0=(Tezker+Teskuin+Tgoi+Tbehe)/4; //ustezko balioa
Tzahar=ones(a,l)*T0;


//KALKULUA 
iter=0
errorea=1000000
    while errorea>tol
  //MATRIZE BERRIA DEFINITZEN
  //-------------------------------------------------------------- 
   tarte=tarte*2;
   l=tarte+1;
   a=tarte+1;
   dx=L/tarte;
   dy=A/tarte;
   iter=iter+1
   T=ones(a,l)*T0;
      for i=2:a-1
      T(i,1)=Tezker;
      T(i,l)=Teskuin;
      end
     for j=2:l-1
      T(1,j)=Tgoi;
      T(a,j)=Tbehe;
    end
    T(1,1)=%nan; T(a,1)=%nan; T(1,l)=%nan; T(a,l)=%nan;

  //Gauss-Siedel algoritmoa diferentzia finitoen metodoaz
  //---------------------------------------------------------------
  err_gauss_seidel=1000000
       while err_gauss_seidel>tol/10
          Tz=T;
          for i=2:a-1
             for j=2:l-1
              T(i,j)=(T(i-1,j)+T(i+1,j)+T(i,j-1)+T(i,j+1))/4;
            end
        end
         err_gauss_seidel=norm(T(2:a-1,2:l-1)-Tz(2:a-1,2:l-1));
     end
  
  //----------------------------------------------------------------   
  [nv,mv]=size(Tzahar);
   for n=1:nv;
       for m=1:mv;
         Tmurriztua(n,m)=T(2*n-1,2*m-1);
       end
   end
   errorea=abs(norm(Tzahar(2:nv-1,2:mv-1))-norm(Tmurriztua(2:nv-1,2:mv-1)));
   Tzahar=T;
   disp(errorea) 
   end 

 T
 contour2d([1:l],[1:a],T',20)
 endfunction
//-------------------------------------------------------------------------------

//XAFLAREN PROPIETATEAK JASOTZEN
L = input("sartu piezaren luzera (zenbaki osoa):");
A = input("sartu piezaren zabalera (zenbaki osoa):");

Tezker = input("plakaren ezker aldeko tenperatura sartu(K):");
Teskuin = input("plakaren eskuin aldeko tenperatura sartu(K):");
Tgoi = input("plakaren goi aldeko tenperatura sartu(K):");
Tbehe = input("plakaren behe aldeko tenperatura sartu(K):");
[T,errorea,iter,a,l,dx]=tenperatura_banaketa3(Tezker,Teskuin,Tgoi,Tbehe,L,A)