import java.io.*;
import visad.*;
import visad.java3d.DisplayImplJ3D;
import visad.java2d.DisplayImplJ2D;
import java.rmi.RemoteException;

class HT_xyt_FD
{   public static double[] pivotlugauss(double a[][],double b[])
  { //kısmi pivotlu gauss eleme yöntemi
  int n=b.length;	  
  double x[]=new double[n];
  double carpan=0;
  double toplam=0;
  double buyuk;
  double dummy=0;
  //gauss eleme
  int i,j,k,p,ii,jj;
  //System.out.println("Orijinal Matris :\n"+Matrix.toString(a));
  for(k=0;k<(n-1);k++)
  { //pivotlama
	p=k;
	buyuk=Math.abs(a[k][k]);
	for(ii=k+1;ii<n;ii++)
	{  dummy=Math.abs(a[ii][k]);
	   if(dummy > buyuk) {buyuk=dummy;p=ii;}
	}
	if(p!=k)
	{ for(jj=k;jj<n;jj++)
	  { dummy=a[p][jj];
	    a[p][jj]=a[k][jj];
	    a[k][jj]=dummy;
	  }
	  dummy=b[p];
	  b[p]=b[k];
	  b[k]=dummy;
	}
	  //System.out.println("pivotlu Matris : k="+k+"\n"+Matrix.toString(a));
	//gauss elemeyi çözme  
    for(i=k+1;i<n;i++)
    {  carpan=a[i][k]/a[k][k];
       a[i][k]=0;
       for(j=k+1;j<n;j++)
       {   a[i][j]-=carpan*a[k][j]; }
           b[i]   =b[i]   -carpan*b[k];
	}
  }
  //geriye doğru yerine koyma
  x[n-1]=b[n-1]/a[n-1][n-1];
  for(i=n-2;i>=0;i--)
  {
    toplam=0;
    for(j=i+1;j<n;j++)
    {  toplam+=a[i][j]*x[j];}
    x[i]=(b[i]-toplam)/a[i][i];
  }
  return x;
  }
	
	public static double [] thomas(double a[][],double b[])
  {
  //3 band matrix
  int n=a.length;
  double f[]=new double[n];
  double e[]=new double[n];
  double g[]=new double[n];
  double x[]=new double[n];
  for(int i=0;i<n;i++)     {f[i]=a[i][i];}
  for(int i=0;i<(n-1);i++) {g[i]=a[i][i+1];}
  for(int i=0;i<(n-1);i++) {e[i+1]=a[i+1][i];}
  for(int k=1;k<n;k++)
     {e[k]=e[k]/f[k-1];
      f[k]=f[k]-e[k]*g[k-1];
     }
  for(int k=1;k<n;k++)
     {b[k]=b[k]-e[k]*b[k-1];
     }
  x[n-1]=b[n-1]/f[n-1];   
  for(int k=(n-2);k>=0;k--)
      {x[k]=(b[k]-g[k]*x[k+1])/f[k];}
  return x;      
  }
  public static double [][] q(double k, double dx,double T[][])
  {  //heat transfer
	  int nx=T[0].length;
   int nt=T.length;
   double qi[][]=new double[nt][nx-1];
     for(int p=0;p<nt;p++)
     {  for(int m=1;m<nx;m++)
        { qi[p][m-1]=-k*(T[p][m]-T[p][m-1])/dx;}
     }
     return qi;
  }
 	       	
  public static double [][][] HT_xytcoz(double fi, double L,double W,double t, double dt,double T0[][],double alpha,double Bix0,double BixM,double Biy0,double BiyN,double Tx0inf,double TxMinf,double Ty0inf,double TyNinf)
  {
  //  L duvar kalınlığı (x boyutu)
  //  W duvar yüksekliği (y boyutu)
  //  t maksimum zaman
  //  dt zaman stebi   
  //  T1 x=0 t>=0 boundary temperature
  //  T2 x=L t>=0 boundary temperature 
  int M=T0[0].length; //x step sayısı
  int N=T0.length;
  int Mm1=M-1;
  int Nm1=N-1;
  //System.out.println(nx="+nx)
  int K=M*N;
  double A[][]=new double[K][K];
  double B[]=new double[K];
  double dx=L/(M-1);
  //it is assumed that dx=dy is provided
  int P=(int)(t/dt+0.001)+1;
  double Fo=alpha*dt/(dx*dx);
  double T[][][]=new double[P][N][M];
  for(int n=0;n<N;n++) 
  {  for(int m=0;m<M;m++)
       {T[0][n][m]=T0[n][m];}
  } 
  int k=0,kp1=0,km1=0;
  for(int p=0;p<P-1;p++)
  { //t>0 wall temperatures
    for(int n=0;n<N;n++)
    {
	for(int m=0;m<M;m++)
    { k=N*n+m;
      kp1=N*(n+1)+m;
      km1=N*(n-1)+m;
	  if(m==0 && n==0) 
	  {B[k]+=2.0*Fo*fi*(T[p][n][m+1]+T[p][n+1][m])+
	  (1.0-4.0*Fo*fi-2.0*Fo*fi*Bix0-2.0*Fo*fi*Biy0)*T[p][n][m]+
	  (2.0*Fo*Bix0*Tx0inf+2.0*Fo*Biy0*Ty0inf);
	  A[k][k]+=(1.0+4.0*(1-fi)*Fo+2.0*Fo*(1-fi)*Bix0+2.0*Fo*(1-fi)*Biy0);
	  A[k][k+1]+=-2.0*(1-fi)*Fo;
	  A[k][kp1]+=-2.0*(1-fi)*Fo;
	  }
	  else if(m==0 && n==Nm1) 
	  {//System.out.println("n="+n+"m="+m+"p="+p);
		  B[k]+=2.0*Fo*fi*(T[p][n][m+1]+T[p][n-1][m])+
		  (1.0-4.0*Fo*fi-2.0*Fo*fi*Bix0-2.0*Fo*fi*BiyN)*T[p][n][m]+
	      (2.0*Fo*Bix0*Tx0inf+2.0*Fo*BiyN*TyNinf);
	  A[k][k]+=(1.0+4.0*(1-fi)*Fo+2.0*Fo*(1-fi)*Bix0+2.0*Fo*(1-fi)*BiyN);
	  A[k][k+1]+=-2.0*(1-fi)*Fo;
	  A[k][km1]+=-2.0*(1-fi)*Fo;
	  }
	  else if(m==Mm1 && n==0) {B[k]+=2.0*Fo*fi*(T[p][n][m-1]+T[p][n+1][m])+
	  (1.0-4.0*Fo*fi-2.0*Fo*fi*BixM-2.0*Fo*fi*Biy0)*T[p][n][m]+
	  (2.0*Fo*BixM*TxMinf+2.0*Fo*Biy0*Ty0inf);
	  A[k][k]+=(1.0+4.0*(1-fi)*Fo+2.0*Fo*(1-fi)*BixM+2.0*Fo*(1-fi)*Biy0);
	  A[k][k-1]+=-2.0*(1-fi)*Fo;
	  A[k][kp1]+=-2.0*(1-fi)*Fo;
	  }
	  else if(m==Mm1 && n==Nm1) {
	  //System.out.println("MM1 k="+k+"Nm1 n="+n+"m="+m+"p="+p);
	  B[k]+=2.0*Fo*fi*(T[p][n][m-1]+T[p][n-1][m])+
	  (1.0-4.0*Fo*fi-2.0*Fo*fi*BixM-2.0*Fo*fi*BiyN)*T[p][n][m]+
	  (2.0*Fo*BixM*TxMinf+2.0*Fo*BiyN*TyNinf);
	  A[k][k]+=(1.0+4.0*(1-fi)*Fo+2.0*Fo*(1-fi)*BixM+2.0*Fo*(1-fi)*BiyN);
	  A[k][k-1]+=-2.0*(1-fi)*Fo;
	  A[k][km1]+=-2.0*(1-fi)*Fo;
	  }
	  else if(m==0) {B[k]+=Fo*fi*(2.0*T[p][n][m+1]+T[p][n+1][m]+T[p][n-1][m])+
	  (1.0-4.0*Fo*fi-2.0*Fo*fi*Bix0)*T[p][n][m]+
	  +2.0*Fo*Bix0*Tx0inf;
	  A[k][k]+=(1.0+4.0*(1-fi)*Fo+2.0*Fo*(1-fi)*Bix0);
	  A[k][kp1]+=-(1-fi)*Fo;
	  A[k][km1]+=-(1-fi)*Fo;
	  A[k][k+1]+=-2.0*(1-fi)*Fo;
	  }
	  else if(m==Mm1) {B[k]+=Fo*fi*(2.0*T[p][n][m-1]+T[p][n+1][m]+T[p][n-1][m])+
	  (1.0-4.0*Fo*fi-2.0*Fo*fi*BixM)*T[p][n][m]+
	  +2.0*Fo*BixM*TxMinf;
	  A[k][k]+=(1.0+4.0*(1-fi)*Fo+2.0*Fo*(1-fi)*BixM);
	  A[k][kp1]+=-(1-fi)*Fo;
	  A[k][km1]+=-(1-fi)*Fo;
	  A[k][k-1]+=-2.0*(1-fi)*Fo;
	  }
	  else if(n==0) {B[k]+=Fo*fi*(T[p][n][m+1]+T[p][n][m-1]+2.0*T[p][n+1][m])+
	  (1.0-4.0*Fo*fi-2.0*Fo*fi*Biy0)*T[p][n][m]+
	  +2.0*Fo*Biy0*Ty0inf;
	  A[k][k]+=(1.0+4.0*(1-fi)*Fo+2.0*Fo*(1-fi)*Biy0);
	  A[k][k+1]+=-(1-fi)*Fo;
	  A[k][k-1]+=-(1-fi)*Fo;
	  A[k][kp1]+=-2.0*(1-fi)*Fo;
	  }
	  else if(n==Nm1) {//System.out.println("Nm1 n="+n+"m="+m+"p="+p);
		  B[k]+=Fo*fi*(T[p][n][m+1]+T[p][n][m-1]+2.0*T[p][n-1][m])+
	  (1.0-4.0*Fo*fi-2.0*Fo*fi*BiyN)*T[p][n][m]+
	  +2.0*Fo*BiyN*TyNinf;
	  A[k][k]+=(1.0+4.0*(1-fi)*Fo+2.0*Fo*(1-fi)*BiyN);
	  A[k][k+1]+=-(1-fi)*Fo;
	  A[k][k-1]+=-(1-fi)*Fo;
	  A[k][km1]+=-2.0*(1-fi)*Fo;
	  }
	  else	  { //internal region
		       B[k]+=Fo*fi*(T[p][n][m-1]+T[p][n][m+1]+T[p][n-1][m]+T[p][n+1][m])+
	           (1-4.0*fi*Fo)*T[p][n][m];
	           A[k][k]+=1.0+4.0*(1-fi)*Fo;
	           A[k][km1]+=-(1-fi)*Fo;
	           A[k][kp1]+=-(1-fi)*Fo;
	           A[k][k-1]+=-(1-fi)*Fo;
	           A[k][k+1]+=-(1-fi)*Fo;	           
	          }
    }
    }
    //Text.print(A);
    //Text.printT(B);
    double aa[]=pivotlugauss(A,B);
    for(int n=0;n<N;n++)
    { for(int m=0;m<M;m++)
       {k=N*n+m;
        //System.out.println("TT aa.length="+aa.length+"P="+P+"N="+N+"M="+M+"k="+k+"p="+p+"n="+n+"m="+m);
        T[p+1][n][m]=aa[k];
       }
    }
  }
  return T;
  }
    
  public static void main(String args[]) throws IOException,RemoteException, VisADException
  {
	  double T0[][]=new double[11][11];
	  double L=0.1;
	  double W=0.1;
	  // t=0 daki T0(x) sıcaklık profili 
	  int N=T0.length;
	  int M=T0[0].length;
	  for(int n=0;n<N;n++) 
	  {  for(int m=0;m<N;m++)
	     {T0[n][m]=20.0;}
	  }
	  Text.print(T0,"T0 main");
	  double Tx0inf=50.0;  //x=0 da t>0 zamanı için duvar sıcaklığı sınır şartı
	  double TxMinf=50.0; 
	  double Ty0inf=500.0;  //x=0 da t>0 zamanı için duvar sıcaklığı sınır şartı
	  double TyNinf=100.0; 
	   //x=L de t>0 zamanı için duvar sıcaklığı sınır şartı
	  double hx0=250.0;  //x=0 da t>0 zamanı için duvar sıcaklığı sınır şartı
	  double hxN=250.0;
	  double hy0=250.0;  //x=0 da t>0 zamanı için duvar sıcaklığı sınır şartı
	  double hyN=250.0;  
	  double t=6000.0; //sn toplam zaman
	  double dt=1.0; //sn gerçek zaman adımı
	  double dtplot=0.3; //sn her bir plotta bekleme zamanı
	  double alpha=8.4e-6;
	  double x0=0.0;
	  double xn=L;
	  double y0=0.0;
	  double yn=W;
	  double fi=0.0;
	  double k=300.0;
	  double dx=L/T0.length;
	  double h=250.0;
	  double Bix0=hx0*dx/k;
	  double BixM=hxN*dx/k;
	  double Biy0=hy0*dx/k;
	  double BiyN=hyN*dx/k;
	   double Fo=alpha*dt/dx;
	  System.out.println("Fo ="+Fo+" Bix0="+Bix0+" BixM="+BixM+" Biy0= "+Biy0+"BiyN= "+BiyN);
	  double T[][][]=HT_xytcoz(fi,L,W,t,dt,T0,alpha,Bix0,BixM,Biy0, BiyN,Tx0inf,TxMinf,Ty0inf, TyNinf);
	  //Text.print(T[1],"1 result");
	  //Text.print(T[10],"200 result");
	  Text.print(T[199],"499 result");
	  //double b[][]=Matrix.T(T);
	  //double q[][]=q(k,dx,T);
	  //Plot3D.data(T,x0,xn,0.0,t,20,"x","t","T=f(x,t)");
	  timeContourPlot pp=new timeContourPlot(T,x0,xn,y0,yn,dt,dt,20,"contour","x","T",true);
	  //timeplot p=new timeplot(T,x0,xn,dtplot,dt);
	  //timeplot p=new timeplot(q,x0,xn,dtplot,dt);
     pp.start();
  }
}