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

class f_duvar extends f_xt
{
double alpha; // thermal diffusivity
double h;     // thermal convectivity
double k;     // thermal conductivity
double Cp;    // Specific heat
double ro;    // density
double Ti;    // initial temperature
double Tinf;  // environmrntal temperature
double r0;    // radius of the cylinder
double Bi;    // Biot number
double C[];   // Coefficient of wall equation
double ln[];  // 
double x0,x1,x2;
int ix;
public f_duvar(double hi,double ki,double Cpi,double roi,double r0i,double Tii,double Tinfi)
{   r0=r0i;
	h=hi;
    k=ki;
    Cp=Cpi;
    ro=roi;
    Ti=Tii;
    Tinf=Tinfi;
	Bi=h*r0/k;
    alpha=k/(ro*Cp);
    ix=0;	
    x1=0.0;
    System.out.println("Bi="+Bi+"alpha="+alpha);
    C=new double[20];
    ln=new double[20];
    for(int n=1;n<20;n++)
    { ln[n]=lambda(Bi);
      double ffx=fx(ln[n],Bi);
      if(Math.abs(ffx)<1.0e-4)
      { C[n]=Cn(ln[n]);
        System.out.println("n ="+n+"  ln="+ln[n]+"C="+C[n]);
      } 
     }
}
public  static double fx(double lambda,double Bi)
{return lambda*Mathd.J(1,lambda)/Mathd.J(0,lambda)-Bi;}

public  static double dfx(double x,double h,double Bi)
{ //Türev hesabı
  return (-fx((x+2.0*h),Bi)+8.0*fx((x+h),Bi)-8.0*fx((x-h),Bi)+fx((x-2.0*h),Bi))/(12.0*h);
}

public  double Cn(double lambda)
{
 double J1=Mathd.J(1,lambda);
 double J0=Mathd.J(0,lambda);	
 return 2.0/lambda*J1/(J0*J0+J1*J1);
}

public void genislet(double Bi)
{
	//kök olan bir bölge bulana kadar verilen sınırı genişletir
	x0=x1;
	x2=x1+0.00001;
	int NTRY=200;
	double a[]=new double[2];
	double FACTOR=1.01;
	int j;
	double f1,f2;
	if (x1 == x2) System.out.println("Giriş değerleri yanlış verildi");
	f1=fx(x1,Bi);
	f2=fx(x2,Bi);
	for (j=1;j<=NTRY;j++) {
		if (f1*f2 < 0.0) {ix++;break;}
		else 
		{f2=fx(x2 += FACTOR*(x2-x1),Bi);}
	}
}
//  x1 ve x2 nin yeni değerleri geldi


public double lambda(double Bi)
{
//lineer olmayan denklemin kökü : ikiyebölme-Newton Raphson metodu
genislet(Bi);         //kök sınırlarını saptar
double xacc=1.0e-10; //iterasyonlar arasındaki fark, bu rakamdan 
// küçük olduğunda çözüme ulaştığımızı kabul edeceğiz
double EPS=1.0e-10;
double h=0.0001; //türev step büyüklüğü
int MAXIT=100;  //Maksimum iterasyon sayısı
int j;
double dfr,dx,dxold,fr,fh,fl;
double temp,b,a,r;
fl=fx(x1,Bi);
fh=fx(x2,Bi);
//System.out.println("x1 ="+x1+"fl="+fl+"x2="+x2+"fh="+fh);
//Bölgemizde kök değeri varmı kontrol et
if ((fl > 0.0 && fh > 0.0) || (fl < 0.0 && fh < 0.0))
System.out.println("n = "+ix+"x1="+x0+"x2="+x2+"verilen sınırlar içinde kök değeri yok");
if (fl == 0.0) return x1;
if (fh == 0.0) return x2;
if (fl < 0.0)   { a=x1;b=x2;} 
else            { b=x1;a=x2;}
r=0.5*(x1+x2); //orta nokta değeri 
dxold=Math.abs(x2-x1); 
dx=dxold; 
fr=fx(r,Bi);   //fonksiyonun r noktasındaki değeri
dfr=dfx(r,h,Bi); //fonksiyonun türevinin r noktasındaki değeri
for (j=1;j<=MAXIT;j++) 
{ 
if ((((r-b)*dfr-fr)*((r-a)*dfr-fr) > 0.0) || (Math.abs(2.0*fr) > Math.abs(dxold*dfr))) 
    { //ikiye bölme 
	  dxold=dx;
      dx=0.5*(b-a);
      r=a+dx;
      if ( a == r && fr>EPS) {x1=r+0.001;return r;} //çözüm!!!!
    } 
else{ //Newton-Raphson
	  dxold=dx;
      dx=fr/dfr;
      temp=r;
      r -= dx;
      if ((temp == r) && fr<EPS) {x1=r+0.001;return r;}  //çözüm!!!!!
    }
if ((Math.abs(dx) < xacc)&& fr<EPS) {x1=r+0.001;return r;} //çözüm!!!!!
fr=fx(r,Bi);
dfr=dfx(r,h,Bi);
if (fr < 0.0) 
a=r;
else
b=r;
}
System.out.println("Maksimum iterasyon sayısı aşıldı");
return 0.0; //program buraya gelmemeli 
}

//***********


public double T(double t,double Bi,double r)
{ double Fo=alpha*t/(r0*r0);
  double tt=0;
  for(int n=1;n<20;n++)
  {   
	 tt+=C[n]*Math.exp(-ln[n]*ln[n]*Fo)*Mathd.J(0,ln[n]*r/r0);
  }
  return Tinf+tt*(Ti-Tinf);	
}

public double func(double x,double t)
{
return T(t,Bi,x);
}
}

public class timeplottest4
{
     public static void main (String args[])  
     { //Material : Carbon steel
	   double h=50;//W/m^2K
       double k=60.5; //W/mK
       double Cp=434.0;//J/kgK
       double ro=7854.0;//kg/m^3
       double L=0.05;
       double Ti=100.0;//C
       double Tinf=20.0;//C
	   f_duvar ff=new f_duvar(h,k,Cp,ro,L,Ti,Tinf);
       timeplot p=new timeplot(ff,0.0,L,0.001,0.0,500.0,10.0,0.1);
       p.start();
     } 
}