import java.io.*;

class HT_cylinder_FD_g
{ 
  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 [] HT_cylinder( double b,double gm[],double dr, double Tr,double hr,double kr)
  {
  //b : cylinder radius 
  //gm : heat generatþon per unit volume 
  int M=(int)(b/dr); 
  int mt=M+1;   
  double A[][]=new double[mt][mt];
  double B[]=new double[mt];
  double x[]=new double[mt];
  double x1,x2,y1,y2;

  // boundary conditions
  int m=0;
  double GM=0,betaM=0,gammaM=0;
  for(m=0;m<=M;m++)
  {  GM=dr*dr*gm[m]/kr;
     betaM=1+(1+1.0/(2.0*m))*dr*hr/kr;
     gammaM=(1+1.0/(2.0*m))*dr*hr/kr*Tr;
	 if(m==0)            {B[m]+=-dr*dr*gm[0]/kr;}
     else if(m==M)       {B[m]+=-2.0*gammaM-GM;}
     else                {B[m]+=-dr*dr*gm[m]/kr;}
     if(m==0)      {A[m][m]+=-4.0;A[m][m+1]+=4.0;}
     else if(m==M) {A[m][m-1]=2.0;A[m][m]+=-2*betaM;}
     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[]=thomas(A,B); 
  return T;  
}  
  public static void main(String args[]) 
  {  int M=10;
     double b=5.0e-2;
     double dr=b/M;
     double gm[]=new double[M+1];
     double r[]=new double[M+1];
     for(int m=0;m<M+1;m++)
     {gm[m]=4.0e6;//W/m3
     r[m]=dr*m;
     }
     double Tr=20.0;//degree C
     double hr=400.0;//W/(m2K)
     double k=40.0;//W/(mK)
	 double T[]=HT_cylinder( b,gm,dr,Tr,hr,k);
	 //compare result with analytical case with qm=const.
	 double T1[]=new double[T.length];
	 for(int m=0;m<M+1;m++)
     {T1[m]=gm[0]*b*b/(4.0*k)*(1.0-(m*dr/b)*(m*dr/b))+gm[0]*b/(2.0*hr)+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();
  }
}