import java.util.*;
import java.awt.*;
import javax.swing.*;
import java.awt.event.*;

class f1 extends f_x
{ 
  double a[][];
  double T;
  double xlimit;
  double sigma=5.670503208914327e-8;
  double C1=3.74177489e8;
  double C2=1.438769e4;   
  public f1(double x[][],double xlimiti,double Ti)
  {   T=Ti;  
	  for(int i=0;i<x[0].length;i++)
        {x[0][i]=T*x[0][i];}
	    a=cubic_spline(x[0],x[1],0.0,0.0);T=Ti;xlimit=C2/xlimiti;
  }
  
  public f1(String name,double xlimiti,double Ti)
  {double x[][]=Text.readDoubleT(name);T=Ti;xlimit=C2/xlimiti;
  T=Ti;
  for(int i=0;i<x[0].length;i++)
  {x[0][i]=T*x[0][i];
   //x[1][i]=x[1][i]/100.0;
  }   
  Text.print(x);
   a=cubic_spline(x[0],x[1],0.0,0.0);
    Text.printT(a);
  }  
  public f1(double xlimiti,double Ti)
  { T=Ti;
    double x[][]=Text.readDoubleT("a.txt");T=Ti;xlimit=C2/xlimiti;
    for(int i=0;i<x[0].length;i++)
    { x[0][i]=T*x[0][i];}  
    Text.printT(x);
    a=cubic_spline(x[0],x[1],0.0,0.0);
    Text.printT(a);
  }
  
  public double func(double y)
  {
   double x=C2/y;  
   double eps=0;
   if(x>xlimit) eps=0.0;
   else eps=funcSpline(a,x);
   if(eps<0) eps=0;
   return eps/sigma*C1/(C2*C2*C2*C2)*y*y*y/(Math.exp(y)-1.0);
  }
  
  public static double [] thomas(double a[][],double r[])
  {
  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++)
     {r[k]=r[k]-e[k]*r[k-1];
     }
  x[n-1]=r[n-1]/f[n-1];   
  for(int k=(n-2);k>=0;k--)
     {x[k]=(r[k]-g[k]*x[k+1])/f[k];}
  return x;      
  }
  
  public static double [] thomas(double f[],double e[],double g[],double r[])
  {
  int n=f.length;
  double x[]=new double[n];
  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++)
     {r[k]=r[k]-e[k]*r[k-1];
     }
  x[n-1]=r[n-1]/f[n-1];   
  for(int k=(n-2);k>=0;k--)
     {x[k]=(r[k]-g[k]*x[k+1])/f[k];}
  return x;      
  }
    
  public static double [][] cubic_spline(double xi[],double yi[],double c0,double cn)
  {
  int n=xi.length;
  double h[]=new double[n];
  double w[]=new double[n];
  double f[]=new double[n];
  double e[]=new double[n];
  double g[]=new double[n];
  double d[]=new double[n];
  double x[]=new double[n];
  double S[][]=new double[4][n];
  int k;	  
  for(k=0;k<(n-1);k++)
     {h[k]=xi[k+1]-xi[k];
      w[k]=(yi[k+1]-yi[k])/h[k];      
	 }
  d[0]=c0;
  d[n-1]=cn;
  for(k=1;k<(n-1);k++)
     {d[k]=6.0*(w[k]-w[k-1]);}
  f[0]=1.0;
  f[n-1]=1.0;
  g[0]=0.0;
  g[n-1]=0.0;
  e[0]=0.0;
  e[n-1]=0.0;
  for(k=1;k<(n-1);k++)
  {f[k]=2.0*(h[k]+h[k-1]);e[k]=h[k-1];g[k]=h[k];}
  S[2]=thomas(f,e,g,d);
  S[3]=xi;
  for(k=0;k<(n-1);k++)
    {S[0][k]=(6.*yi[k+1]-h[k]*h[k]*S[2][k+1])/(6.0*h[k]);
     S[1][k]=(6.*yi[k]-h[k]*h[k]*S[2][k])/(6.0*h[k]);
    } 
  return S;    	 	     	 
  }
  
  public static double funcSpline(double S[][],double x)
  {
  int n=S[0].length;
  double xx1=0;
  double xx2=0;
  double y=0;
  double hk=0;
  for(int k=0;k<(n-1);k++)
     {if(S[3][k]<=x && x<=S[3][k+1]) 
      {hk=(S[3][k+1]-S[3][k]);
       xx1=(x-S[3][k]);
       xx2=(S[3][k+1]-x);
       y=S[0][k]*xx1+S[1][k]*xx2+(xx1*xx1*xx1*S[2][k+1]+xx2*xx2*xx2*S[2][k])/(6.0*hk);
       break;
       }
     }
       if(y==0 && S[3][n-2]<=x )
       {    
	   int k=n-2;
	   hk=(S[3][k+1]-S[3][k]);
       xx1=(x-S[3][k]);
       xx2=(S[3][k+1]-x);
       y=S[0][k]*xx1+S[1][k]*xx2+(xx1*xx1*xx1*S[2][k+1]+xx2*xx2*xx2*S[2][k])/(6.0*hk);
       }
  return y;
  }    
   
 public static double[][]  funcSpline(double xi[],double yi[],int aradegersayisi)
 {
	//aradegersayisi:  x--o--o--x--o--o--x zincirinde x deneysel noktalar ise
	// ara değer sayısı 2 dir
	int n=xi.length;
	int nn=(n-1)*(aradegersayisi+1)+1;
	double z[][]=new double[2][nn];
    double S[][]=cubic_spline(xi,yi,0,0);	
    double dx=0;
    int k=0;
    int i;
    for(i=0;i<(n-1);i++)
       {  z[0][k]=xi[i];z[1][k]=funcSpline(S,z[0][k]);k++;
          for(int j=0;j<aradegersayisi;j++) 
             {dx=(xi[i+1]-xi[i])/((double)aradegersayisi+1.0);
          z[0][k]=z[0][k-1]+dx;z[1][k]=funcSpline(S,z[0][k]);k++;}
       }
       z[0][k]=xi[i];z[1][k]=funcSpline(S,z[0][k]);
    return z;
 }
 public static void plotSpline(double xi[],double yi[],int aradegersayisi)
 {	//aradegersayisi:  x--o--o--x--o--o--x zincirinde x deneysel noktalar ise
	// ara değer sayısı 2 dir
	int n=xi.length;
	int nn=(n-1)*(aradegersayisi+1)+1;
	double z[][]=new double[2][nn];
    double S[][]=cubic_spline(xi,yi,0,0);	
    double dx=0;
    int k=0;
    int i;
    for(i=0;i<(n-1);i++)
       {  z[0][k]=xi[i];z[1][k]=funcSpline(S,z[0][k]);k++;
          for(int j=0;j<aradegersayisi;j++) 
             {dx=(xi[i+1]-xi[i])/((double)aradegersayisi+1.0);
          z[0][k]=z[0][k-1]+dx;z[1][k]=funcSpline(S,z[0][k]);k++;}
       }
       z[0][k]=xi[i];z[1][k]=funcSpline(S,z[0][k]);
       Plot pp=new Plot(z[0],z[1]);
       pp.setPlotType(0,0); 
       pp.addData(xi,yi);
       pp.setPlotType(1,24);
       //pp.setPlabel("R134a");
       //pp.setXlabel("T derece C");
       //pp.setYlabel("h ");
       pp.plot();  
       }
 // Türev formülleri
// =================       
  public static double dfSpline(double xi[],double yi[],double c0,double cn,double x)
  { //kübik şerit türev formülü 
	double S[][]=cubic_spline(xi,yi,c0,cn);	
    return dfSpline(S,x);
  }
  

  
  public static double dfSpline(double S[][],double x)
  {
    //kübik şerit türev formülü 
  int n=S[0].length;
  double xx1=0;
  double xx2=0;
  double y=0;
  double hk=0;
  for(int k=0;k<(n-1);k++)
     {if(S[3][k]<=x && x<=S[3][k+1]) 
      {hk=(S[3][k+1]-S[3][k]);
       xx1=(x-S[3][k]);
       xx2=(S[3][k+1]-x);
       y=S[0][k]-S[1][k]+(xx1*xx1*S[2][k+1]-xx2*xx2*S[2][k])/(2.0*hk);
       break;
       }
     }
       if(y==0 && S[3][n-2]<=x )
       {    
	   int k=n-2;
	   hk=(S[3][k+1]-S[3][k]);
       xx1=(x-S[3][k]);
       xx2=(S[3][k+1]-x);
       y=S[0][k]*xx1+S[1][k]*xx2+(xx1*xx1*xx1*S[2][k+1]+xx2*xx2*xx2*S[2][k])/(6.0*hk);
       }
  return y;
  } 
       
 public static double  dfSpline(double xi[],double yi[],double x)
 {
    //doğal kübik şerit türev formülü 
	double S[][]=cubic_spline(xi,yi,0,0);	
    return dfSpline(S,x);
 }
  
  public static double intSpline(double S[][],double a,double b)
  {
    //kübik şerit türev formülü 
  int n=S[0].length;
  double xx1=0;
  double xx2=0;
  double y1,y2;
  double hk=0;
  double toplam=0;
  y1=0;y2=0;
  for(int k=0;k<(n-1);k++)
     { hk=(S[3][k+1]-S[3][k]);
       //System.out.println("x1="+S[3][k]+"x2="+S[3][k+1]+"a="+a+"b="+b);
      if(a>S[3][k+1])
      {toplam=0;//System.out.println("şart 1 toplam = "+toplam);
      }//şart1
	  else if(S[3][k]<=a && a<=S[3][k+1] && b>S[3][k+1])//şart 2 
       {
       xx1=(a-S[3][k]);
       xx2=(S[3][k+1]-a);
       y1=S[0][k]*xx1*xx1/2.0-S[1][k]*xx2*xx2/2.0+(xx1*xx1*xx1*xx1*S[2][k+1]-xx2*xx2*xx2*xx2*S[2][k])/(24.0*hk);
       xx1=hk;
       xx2=0;
       y2=S[0][k]*xx1*xx1/2.0-S[1][k]*xx2*xx2/2.0+(xx1*xx1*xx1*xx1*S[2][k+1]-xx2*xx2*xx2*xx2*S[2][k])/(24.0*hk);
       toplam+=(y2-y1);
       //System.out.println("şart 2toplam = "+toplam+"y1="+y1+"y2="+y2);
       }
       else if(S[3][k]<=a && a<=S[3][k+1] && S[3][k]<=b && b<=S[3][k+1])//şart 3 
       {
       xx1=(a-S[3][k]);
       xx2=(S[3][k+1]-a);
       y1=S[0][k]*xx1*xx1/2.0-S[1][k]*xx2*xx2/2.0+(xx1*xx1*xx1*xx1*S[2][k+1]-xx2*xx2*xx2*xx2*S[2][k])/(24.0*hk);
       xx1=(b-S[3][k]);
       xx2=(S[3][k+1]-b);
       y2=S[0][k]*xx1*xx1/2.0-S[1][k]*xx2*xx2/2.0+(xx1*xx1*xx1*xx1*S[2][k+1]-xx2*xx2*xx2*xx2*S[2][k])/(24.0*hk);
       toplam+=(y2-y1);
       //System.out.println("şart 3 toplam = "+toplam+"y1="+y1+"y2="+y2);
       }
       else if(a<S[3][k] && b>=S[3][k+1]) //şart 4
       {
       xx1=0;
       xx2=hk;
       y1=S[0][k]*xx1*xx1/2.0-S[1][k]*xx2*xx2/2.0+(xx1*xx1*xx1*xx1*S[2][k+1]-xx2*xx2*xx2*xx2*S[2][k])/(24.0*hk);
       xx1=hk;
       xx2=0;
       y2=S[0][k]*xx1*xx1/2.0-S[1][k]*xx2*xx2/2.0+(xx1*xx1*xx1*xx1*S[2][k+1]-xx2*xx2*xx2*xx2*S[2][k])/(24.0*hk);
       toplam+=(y2-y1);
       //System.out.println("şart 4 toplam = "+toplam+"y1="+y1+"y2="+y2);
       }
       else if(S[3][k]<=b && b<=S[3][k+1])//şart 5
       {xx1=0;
        xx2=hk;
        y1=S[0][k]*xx1*xx1/2.0-S[1][k]*xx2*xx2/2.0+(xx1*xx1*xx1*xx1*S[2][k+1]-xx2*xx2*xx2*xx2*S[2][k])/(24.0*hk);
        xx1=(b-S[3][k]);
        xx2=(S[3][k+1]-b);
        y2=S[0][k]*xx1*xx1/2.0-S[1][k]*xx2*xx2/2.0+(xx1*xx1*xx1*xx1*S[2][k+1]-xx2*xx2*xx2*xx2*S[2][k])/(24.0*hk);
        toplam+=(y2-y1);
        //System.out.println("şart 5 toplam = "+toplam+"y1="+y1+"y2="+y2);
        }
        else break;
     }
  return toplam;
  }
    
  public static double intSpline(double xi[],double yi[],double c0,double cn,double a,double b)
  {double S[][]=cubic_spline(xi,yi,c0,cn);	
   return intSpline(S,a,b);
  }  
}

public class emissivity 
{
public static double emissivity(double xlimiti,double Ti)
{   
	int i;
	double error=0.001;
	int MAX_ITERATION=100;
	f1 ff=new f1(xlimiti,Ti);
	double a=1.0;
	double b=xlimiti;
    double h=(b-a);
    double h2=h/4.0;
	double x[]=new double[5];
    double f[]=new double[5];
    double w1[]=new double[3];
    double z1[]=new double[3]; 
    int j;
    for(i=0;i<5;i+=2) {j=i/2;x[i]=a+i*h2;f[i]=ff.func(x[i]);w1[j]=x[i];z1[j]=f[i];}	
    double sum=h*(f[0]+4.0*f[2]+f[4])/6.0; 
	double Q=0;
	Plot pp=new Plot(ff,a,b,100);
	pp.setPlotType(0,0);// 	 
	pp.addData(w1,z1,28);
	int k=0;
    Q=simpson1_3(ff,x,f,sum,error,k,MAX_ITERATION,Q,pp); 
    pp.plot();
    return Q;	
}
	
public static double simpson1_3(f_x ff,double x[],double f[],double sum,double eps,int k,int MAX_ITERATION,double Q,Plot pp)
{
        //adaptive Simpson 1/3 rule
      
        double h=(x[4]-x[0]);
        double h1=h/2.0;
        double h2=h/4.0;        
        int p,i;
        double x1[]=new double[5];
        double f1[]=new double[5];
        double x2[]=new double[5];
        double f2[]=new double[5];
        double w1[]=new double[2];
        double z1[]=new double[2]; 
        int j;
        for(i=1;i<5;i+=2) {j=i/2;x[i]=x[0]+i*h2;w1[j]=x[i];f[i]=ff.func(x[i]);z1[j]=f[i];}           
        for(i=0;i<5;i++) {System.out.println(i+" "+x[i]+" "+f[i]);}
        {System.out.println("==================");}
        double sum1=h1*(f[0]+4.0*f[1]+f[2])/6.0;
        double sum2=h1*(f[2]+4.0*f[3]+f[4])/6.0; 
        if(1.0/15.0*Math.abs(sum1+sum2-sum)<eps || k>MAX_ITERATION) Q+=sum1+sum2;        
        else 
        { x1[0]=x[0];f1[0]=f[0];
          x1[2]=x[1];f1[2]=f[1];
          x1[4]=x[2];f1[4]=f[2];
          x2[0]=x[2];f2[0]=f[2];
          x2[2]=x[3];f2[2]=f[3];
          x2[4]=x[4];f2[4]=f[4]; 
	      k++;
	      Q+=simpson1_3(ff,x1,f1,sum1,eps,k,MAX_ITERATION,Q,pp)+
             simpson1_3(ff,x2,f2,sum2,eps,k,MAX_ITERATION,Q,pp);
        }
        pp.addData(w1,z1,28);
        return Q;
}
  public static double F(double x)
  {//Black body shape function
  double C2=1.438769e4;
  double y=C2/x;
  double total=0;
  for(int n=1;n<40;n++)
  {total+=Math.exp(-n*y)/n*(y*y*y+3.0*y*y/n+6.0*y/(n*n)+6.0/(n*n*n));}  
  double pi=Math.PI;
  total=15.0/(pi*pi*pi*pi)*total;
  return total;
  }
  
  public double total_emmissivity(double xlimit, double Ti)
  { f1 b1=new f1(xlimit,Ti);
    int n=1000;
	double x=0; //lambda*T
    double C2=1.438769e4;  
    double y=0; 
    double fo=0.0;
    double epso=b1.funcSpline(b1.a,0.0);
    double fn=0;
    double eps=0;
    double total=0;
    double dx=xlimit*Ti/n;
    for(int i=0;i<n;i++) 
    { x+=dx;y=C2/x;fn=F(x);
      if(x>xlimit) eps=0;
      else if(x<0) eps=0;
      else eps=b1.funcSpline(b1.a,x);
      total+=eps*(fn-fo);
      fo=fn;epso=eps;
      System.out.println("i="+i+"x="+x+"fn="+fn+"fn-1="+fo+"eps="+eps+"total="+total);
    }
    return total;
  }

      public static void main(String[] args)
      { double x=2.0*1600.0;
	    f1 b1=new f1("a.txt",x,1600.0);
	   for(double y=5;y<x;y+=1000.0)
	   {System.out.println(y+" "+b1.func(y));}
       Plot pp=new Plot(b1,0.0,x,8);
	   pp.plot();
	   emissivity e=new emissivity();
       System.out.println("F(900)="+F(900.0));
       System.out.println("F(5590)="+F(5000.0));
       double et=0.4*F(4110.0)+0.8*( F(5590.0) - F(4110.0) );
       System.out.println("total emissivity ="+et);
       System.out.println("total emissivity(integral) ="+e.total_emmissivity(1000.0,400.0));
      }
}
