public class radiation_exchange
{
double T[];
double q[];
double eps[];
double F[][];
double J[];
double sigma=5.67e-8; 
public radiation_exchange(double Ti[],double qi[],double epsi[],double Fi[][])
{
int k=Ti.length;
int w=qi.length;
int n=k+w;
T=new double[k];
eps=new double[k];
q=new double[w];
F=new double[n][n];
J=new double[n];
for(int i=0;i<k;i++)
  {T[i]=Ti[i];
   eps[i]=epsi[i];
   for(int j=0;j<n;j++)
    {F[i][j]=Fi[i][j];}
  }
for(int i=k;i<n;i++)
  {q[i-k]=qi[i-k];
   for(int j=0;j<n;j++)
    {F[i][j]=Fi[i][j];}
  }  
calculate();	
}
  public static double[] gauss(double a[][],double b[])
  { //gauss elimination with partial pivoting
  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("Original matrix :\n"+Matrix.toString(a));
  for(k=0;k<(n-1);k++)
  { //pivoting
	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("pivoted matrix : k="+k+"\n"+Matrix.toString(a));
	//Solving gauss elimination
    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];
	}
  }
  //System.out.println("final upper triangular matrix : \n"+Matrix.toString(a));
  //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 int dirac_delta(int i, int j)
  {
   int k=0;
   if(i==j) k=1;
   return k;	  
  }
  public void calculate()
  {   int n=F.length; 
      int k=T.length;
      int w=q.length;
	  double M[][]=new double[n][n];
	  double S[]=new double[n];
	  for(int i=0;i<k;i++)
      {   for(int j=0;j<n;j++)
           { M[i][j]=(dirac_delta(i,j)-(1.0-eps[i])*F[i][j])/eps[i];}
           S[i]=sigma*T[i]*T[i]*T[i]*T[i];
      }
      for(int i=k;i<n;i++)
      {   for(int j=0;j<n;j++)
           {  M[i][j]=dirac_delta(i,j)-F[i][j];}
           S[i]=q[i-k];
      }     
      J=gauss(M,S);
  }
  public double[] q()
  {int n=F.length;
   int k=T.length; 
   double qi[]=new double[n];
   for(int i=0;i<k;i++)
      {qi[i]=eps[i]/(1.0-eps[i])*(sigma*T[i]*T[i]*T[i]*T[i]-J[i]);}
   for(int i=k;i<n;i++)
      {   System.out.println("i="+i);
	      qi[i]=q[i-k];}
   return qi;
  }
}