import java.io.*;

class HT_cylinder_FD_g1
{ 

  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];
	}
  }
  //back substitution
  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_cylinder( double b,double gm[],double dr, double Tr,double k)
  {
  //b : cylinder radius 
  //gm : heat generatşon per unit volume 
  int M=(int)(b/dr+0.01);    
  double A[][]=new double[M][M];
  double B[]=new double[M];
  double x[]=new double[M+1];
  double x1,x2,y1,y2;

  // boundary conditions
  int m=0;
  for(m=0;m<M;m++)
  { 
	 if(m==0)            {B[m]+=-dr*dr*gm[m]/k;}
     else if(m==(M-1))   {B[m]+=-dr*dr*gm[m]/k-(1+1/(2.0*m))*Tr;}
     else                {B[m]+=-dr*dr*gm[m]/k;}
     if(m==0)      {A[m][m]+=-4.0;A[m][m+1]+=4.0;}
     else if(m==(M-1)) {A[m][m-1]+=1.0-1/(2.0*m);A[m][m]+=-2.0;}
     else          {A[m][m-1]+=1.0-1/(2.0*m);A[m][m]+=-2.0;A[m][m+1]+=1.0+1/(2.0*m);}
  }
  Text.print(A,"A");
  Text.printT(B,"B");
  //double T[]=conjugate_gradient1.conjugate_gradient(A,B);     
  double T[]=pivotlugauss(A,B); 
  return T;  
}  
  public static void main(String args[]) 
  {  int M=10;
     double b=1.0e-2;
     double dr=b/M;     
     double gm[]=new double[M+1];
     double r[]=new double[M+1];
     for(int m=0;m<=M;m++)
     {gm[m]=5.0e8;//W/m3
      r[m]=dr*m;
     }
     double Tr=100.0;//degree C
     double k=25.0;//W/(mK)
     System.out.println("dr="+dr+"gg="+gm[0]*dr*dr/k);
	 double a[]=HT_cylinder( b,gm,dr,Tr,k);
	 double T[]=new double[a.length+1];
	 //exact solution
	 double T1[]=new double[a.length+1];
	 T[M]=Tr;
	 T1[M]=Tr;	 
	 for(int m=0;m<M;m++)
     {T[m]=a[m];     
      T1[m]=gm[m]*b*b/(4.0*k)*(1.0-(r[m]/b)*(r[m]/b))+Tr;
     }
     double TT[][]={T,T1};
     Text.print(TT);
	 Plot pp=new Plot(r,T);
	 pp.addData(r,T1);
	 pp.setPlabel("heat transfer with heat generation");
	 pp.setXYlabel("r (m)","T degree c");
	 pp.plot();
  }
}