import java.io.*;
/*
import visad.*;
import visad.java3d.DisplayImplJ3D;
import visad.java2d.DisplayImplJ2D;
import java.rmi.RemoteException;
*/
class HT_plane_FD1
{ 
  public static double [] gauss_seidel(double a[][],double b[],double x[],double lambda)
  {
  //denklem sistemi çözümü b iteratif gauss_seidel metodu metodu
  //yakınsama faktörü lambda ile birlikte verilmiş formül	  
  int n=a.length;
  int imax=1000;//maksimum iterasyon sayısı 
  double es=1.0e-10;
  double dummy=0;
  double ea=1e-1;
  double sum=0;
  double old=0;
  int iter=0;
  int sentinel=0;
  //System.out.print("iter = "+(iter+1)+" ");
  for(int i=0;i<n;i++)
  {  dummy=a[i][i];
     for(int j=0;j<n;j++)
       {a[i][j]/=dummy;}
     b[i]/=dummy;
  }
  for(int i=0;i<n;i++)
  {  sum=b[i];
     for(int j=0;j<n;j++)
       {if(i!=j)sum-=a[i][j]*x[j];}
     x[i]=sum;
     //System.out.print(x[i]+ "  ");  
  }
   //System.out.println();  
  iter++;
  do
  {
  sentinel=1;
  //System.out.print("iter = "+(iter+1)+" ");  
  for(int i=0;i<n;i++)
     { 
	 old=x[i];
	 sum=b[i];
	 for(int j=0;j<n;j++)
	    {if(i!=j) sum-=a[i][j]*x[j];}
	    x[i]=lambda*sum+(1.0-lambda)*old;
	    if(sentinel==1 && x[i]!=0.0 ) {ea=Math.abs((x[i]-old)/x[i])*100.0;}
	    if(ea<es) sentinel=0;
	  //System.out.print(x[i]+ "  ");     
     }
  //System.out.println();  
  iter++;
  }while(sentinel==1 && iter<imax);  
  if(iter>=imax) System.out.println("Maksimum iterasyon sayısı aşıldı sonuç geçerli olmayabilir iter = "+iter);
  return x;
  }
  
  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 [][] HT_planecoz( double L,double W,double dx, double Tx0,double Txn,double Ty0,double Tyn,double hx0,double hxn,double hy0,double hyn,double kc)
  {
  //  Assume dx=dy=const;
  //  Temperature of 4 sides are known
  int nx=(int)(L/dx+0.01)+1;
  int ny=(int)(W/dx+0.01)+1;
  System.out.println("nx="+nx+"ny="+ny);
  //  teta = (t-Ty0)/(Tyn-Ty0)
  double tetay0=0.0;
  double tetayn=1.0;
  double tetax0=(Tx0-Ty0)/(Tyn-Ty0);
  double tetaxn=(Txn-Ty0)/(Tyn-Ty0);
  int mx=nx;
  int my=ny;
  int mt=mx*my;    
  double A[][]=new double[mt][mt];
  double b[]=new double[mt];
  double x[]=new double[mt];
  double x1,x2,y1,y2;
  // boundary conditions
  int i=0,j=0;
  for(i=0;i<my;i++)
  {   
     for(j=0;j<mx;j++)
     {
	 x1=dx*j;
	 x2=(L-x1);
	 y1=dx*i;
	 y2=(W-y1);    
	 int k=i*my+j;
	 x[k]=((x1*tetax0+x2*tetaxn)/L+(y1*tetay0+y2*tetayn)/W)/2.0;	 
	int kp1=(i+1)*my+j;
	int km1=(i-1)*my+j;
	 // 1
     if((i==0)&&(j==0)) {A[k][k]+=-(hy0*dx/kc+hx0*dx/kc+2.0);A[kp1][k]+=1.0;A[k][k+1]+=1.0;b[k]+=-hy0*dx/kc*tetay0-hx0*dx/kc*tetax0;}
     // 4
     else if((i==0)&&(j==(mx-1))) {A[k][k]+=-(hxn*dx/kc+hy0*dx/kc+2.0);A[kp1][k]+=1.0;A[k][k-1]+=1.0;b[k]+=-hxn*dx/kc*tetaxn-hy0*dx/kc*tetay0;}
     // 2
     else if((i==(my-1))&&(j==0)) {A[k][k]+=-(hx0*dx/kc+hyn*dx/kc+2.0);A[km1][k]+=1.0;A[k][k+1]+=1.0;b[k]+=-hyn*dx/kc*tetayn-hx0*dx/kc*tetax0;}
     // 3
     else if((i==(my-1))&&(j==(mx-1))) {A[k][k]+=-(hxn*dx/kc+hyn*dx/kc+2.0);A[km1][k]+=1.0;A[k][k-1]+=1.0;b[k]+=-hyn*dx/kc*tetayn-hxn*dx/kc*tetaxn;}
     // 7
     else if(j==0) {A[k][k]+=-(2.0*hx0*dx/kc+4.0);A[kp1][k]+=1.0;A[km1][k]+=1.0;A[k][k+1]+=2.0;b[k]+=-2.0*hx0*dx/kc*tetax0;}
     // 5
     else if(i==0) {A[k][k]+=-(2.0*hy0*dx/kc+4.0);A[kp1][k]+=2.0;A[k][k-1]+=1.0;A[k][k+1]+=1.0;b[k]+=-2.0*hy0*dx/kc*tetay0;}
     // 8
     else if(j==(mx-1)) {A[k][k]+=-(2.0*hxn*dx/kc+4.0);A[kp1][k]+=1.0;A[km1][k]+=1.0;A[k][k-1]+=2.0;b[k]+=-2.0*hxn*dx/kc*tetaxn;}
     // 6
     else if(i==(my-1)) {A[k][k]+=-(2.0*hyn*dx/kc+4.0);A[k][k-1]+=1.0;A[k][k+1]+=1.0;A[km1][k]+=2.0;b[k]+=-2.0*hyn*dx/kc*tetayn;}
     // iç bölge
     else {A[k][k]+=-4.0;A[km1][k]+=1.0;A[kp1][k]+=1.0;A[k][k-1]+=1.0;A[k][k+1]+=1.0;}
     }
  }
  //Text.print(A);
  double T[][]=new double[ny][nx];   
  //double a[]=conjugate_gradient1.conjugate_gradient(x,A,b);
  double a[]=gauss_seidel(A,b,x,1.2);
  //double a[]=pivotlugauss(A,b);
    for(i=0;i<a.length;i++) {a[i]=a[i]*(Tyn-Ty0)+Ty0;}
  for(i=0;i<ny;i++)
  { for(j=0;j<nx;j++)
   { 
    {int k=i*my+j;
     T[i][j]=a[k];
    }
   }
  }

  return T;
  }
       
  public static void main(String args[]) //throws IOException,RemoteException, VisADException
  {   //HT_planecoz(L,W,dx,Tx0,Txn,Ty0,Tyn,hx0,hxn,hy0,hyn,kc)
	  double a[][]=HT_planecoz(0.4,0.4,0.01,20.0,20.0,20.0,20.0,200.0,200.0,200.0,200.0,200.0);
	  Text.print(a);
	  ContourPlot pp=new ContourPlot(a,0.0,0.4,0.0,0.4,40, false,"x","y","T=f(x,y)");
	  pp.plot();
	  //Plot3D.data(a,0.0,0.4,0.0,0.4,400,"x","y","T=f(x,y)");
  }
}