class fx extends f_x
{ double beta,b;
  public fx(double betai,double bi)
  {beta=betai;
   b=bi;
  }
  
  public void setbeta(double betai)
  {beta=betai;}
  
 public double func(double r)
 {double y=r*Mathd.J(0,beta*r)*20.0*Math.cos(Math.PI/2.0*r/b);
  return y;
 }
}

public class HT_EX3_1
{   double x0,x1,x2;
    int ix;
    double beta[];
    double C[];
    double b,H;
    double alpha;
    public HT_EX3_1(double bi,double Hi,double alphai)
	{ix=0;	
     x1=0.0;
	 b=bi;
	 H=Hi;
	 alpha=alphai;
	 x0=0;
	 beta=new double[10];
     C=new double[10];
     for(int n=1;n<=10;n++)
     {beta[n-1]=beta();
     C[n-1]=Cm(n-1);
     }
   }	
    
public  double fx(double beta)
{return beta*Mathd.J(1,beta*b)-H*Mathd.J(0,beta*b);}

public  double dfx(double x,double h)
{ //Türev hesabı
  return (-fx(x+2.0*h)+8.0*fx(x+h)-8.0*fx(x-h)+fx(x-2.0*h))/(12.0*h);
}
public void genislet()
{
	//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);
	f2=fx(x2);
	for (j=1;j<=NTRY;j++) {
		if (f1*f2 < 0.0) {ix++;break;}
		else 
		{f2=fx(x2 += FACTOR*(x2-x1));}
	}
	System.out.println("x0="+x0+"x1="+x1+"x2="+x2);
}
//  x1 ve x2 nin yeni değerleri geldi


public double beta()
{
//lineer olmayan denklemin kökü : ikiyebölme-Newton Raphson metodu
genislet();         //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);
fh=fx(x2);
//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);   //fonksiyonun r noktasındaki değeri
dfr=dfx(r,h); //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);
dfr=dfx(r,h);
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 N(int i)
{  double J0=Mathd.bessj0(beta[i]*b);
   double H1=(H*H+beta[i]*beta[i]);
   double bb=2.0*beta[i]*beta[i]/(J0*J0*b*b)/H1;
   System.out.println("J0="+J0+"H1="+H1+"bb="+bb);
   return bb;
}

public static double[][] gauss_legendre_coefficients(double x1,double x2,int n)
{
	//calculates legendre gauss-coefficients as coefficients of the integral
	//for n terms
	double EPS=3.0e-15;
	int m,j,i;
	double z1,z,xm,xl,pp,p3,p2,p1;
    //double x[]=new double[n];
    //double w[]=new double[n];
    double a[][]=new double[2][n];//a[0][i]=x[i]   a[1][i]=w[i]
	m=(n+1)/2;
	xm=0.5*(x2+x1);
	xl=0.5*(x2-x1);
	for (i=1;i<=m;i++)  {
		z=Math.cos(Math.PI*((i-0.25)/(n+0.5)));
		do {
			p1=1.0;
			p2=0.0;
			for (j=1;j<=n;j++) {
				p3=p2;
				p2=p1;
				p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
			}
			pp=n*(z*p1-p2)/(z*z-1.0);
			z1=z;
			z=z1-p1/pp;
		} while (Math.abs(z-z1) > EPS);
		a[0][i-1]=xm-xl*z;
		a[0][n-i]=xm+xl*z;
		a[1][i-1]=2.0*xl/((1.0-z*z)*pp*pp);
		a[1][n-i]=a[1][i-1];
			}
return a;
}

public static double integral(f_x f_xnt,double x1,double x2,int n)
{
//n : number of integral coefficients
// this routine first generates gauss legendre coefficients
// for [x1,x2] band 
// then calculates gauss legendre integral
double a[][]=new double[2][n];
a=gauss_legendre_coefficients(x1,x2,n);
//Text.print(Text.T(a));
//System.out.println(Matrix.toString(a));
double z=0;
for(int i=0;i<n;i++)
{ z+=a[1][i]*f_xnt.func(a[0][i]);}
return z;
}

public double Cm(int i)
{   fx f=new fx(b,beta[i]);
System.out.println("i="+i+"ni"+N(i)+"integ="+integral(f,0,b,50));
	double cc=N(i)*integral(f,0,b,50);
return cc;
}

public double T(double r,double t)
{
double T=0;
for(int n=1;n<=10;n++)
{
 T+=C[n-1]*Math.exp(-alpha*beta[n-1]*beta[n-1]*t)*Mathd.J(0,beta[n-1]*r);
}
return T;	
}
public static void main(String arg[])
{ HT_EX3_1 w=new HT_EX3_1(0.1,0.68,2.034e-5);
  for(int n=0;n<10;n++)
  {System.out.println("beta["+n+"] = "+w.beta[n]+"Cm="+w.C[n]);}
}
}