//======================================================
// Numerical Analysis package in java
// Gauss integration of four dimensional function 
// Hamza TAŞ
// =====================================================
import java.io.*;
abstract class fdi4
{
public double a,b;
public fdi4(double ai,double bi)
{a=ai;b=bi;}
public fdi4()
{a=0;b=1;}
public void setab(double ai,double bi)
{a=ai;b=bi;}	
abstract public  double c(double x);
abstract public  double d(double x);
abstract public  double alfa(double x[]);
abstract public  double beta(double x[]);
abstract public  double gama(double x[]);
abstract public  double teta(double x[]);
abstract public  double f(double x[]); 
}

class f1 extends fdi4 
{ double h,R1,R2;
  public f1(double hi,double R1i, double R2i)
  {h=hi;R1=R1i;R2=R2i;a=0;b=R1;}
  
  public double c(double x)
  {return 0;}
  
  public double d(double x)
  {return R2;}
  
   public double alfa(double x[]) 
  {return 0;}
  
  public double beta(double x[]) 
  {return 2.0*Math.PI;}
  
  public double gama(double x[]) 
  {return 0.0;}
  
  public double teta(double x[]) 
  {return 2.0*Math.PI;}
  
  public double f(double x[]) 
  {
  //çözümü istenen fonksiyon	
  double ff;
  double z=(h*h+x[0]*x[0]+x[1]*x[1]-2.0*x[0]*x[1]*Math.cos(x[2]-x[3]));
  ff=1.0/(Math.PI*Math.PI*R1*R1)*h*h*x[0]*x[1]/(z*z);
  return ff; 
  }
  public double  F12(double h,double R1,double R2)
  {double X=(h*h+R1*R1+R2*R2)/(2.0*R1*R2);
   double F=R2/R1*(X-Math.sqrt(X*X-1));
   return F;
  }
}  

class NA78A
{	
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 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 I(fdi4 f,int n,int m,int p,int s) 
{
 int i,j,k,h;
 double rn[],cn[];
 double rm[],cm[];
 double rp[],cp[];
 double rs[],cs[];
 rn=new double[n];
 rm=new double[m];
 rp=new double[p];
 rs=new double[s];
 cn=new double[n];
 cm=new double[m];
 cp=new double[p];
 cs=new double[s];
 double a[][]=new double[2][n];
 a=gauss_legendre_coefficients(-1.0,1.0,n);
 rn=a[0];
 cn=a[1];
 a=gauss_legendre_coefficients(-1.0,1.0,m);
 rm=a[0];
 cm=a[1];
 a=gauss_legendre_coefficients(-1.0,1.0,p);
 rp=a[0];
 cp=a[1];
 a=gauss_legendre_coefficients(-1.0,1.0,s);
 rs=a[0];
 cs=a[1];
double x[],h1,h2,J,a1,b1,c1,d1,k1,k2,JX,JY,JZ,Q,l1,l2,alfa1,beta1,gama1,teta1,m1,m2;
x=new double[4];
a1=f.a;
b1=f.b;
h1=(b1-a1)/2.0;
h2=(b1+a1)/2.0;
J=0;
for(i=0;i<m;i++)
{
JX=0;
x[0]=h1*rm[i]+h2;
d1=f.d(x[0]);
c1=f.c(x[0]);
k1=(d1-c1)/2.0;
k2=(d1+c1)/2.0;
      for(j=0;j<n;j++)
      {   
	  JY=0;
      x[1]=k1*rn[j]+k2;
      beta1=f.beta(x);
      alfa1=f.alfa(x);
      l1=(beta1-alfa1)/2.0;
      l2=(beta1+alfa1)/2.0;
            for(k=0;k<p;k++)
            {
	        JZ=0;
	        x[2]=l1*rp[k]+l2;
	        teta1=f.teta(x);
	        gama1=f.gama(x);
	        m1=(teta1-gama1)/2.0;
	        m2=(teta1+gama1)/2.0;
	               for (h=0;h<s;h++)
	               {
		           x[3]=m1*rs[h]+m2;
	               Q=f.f(x);
	               JZ+=cs[h]*Q;
                   }
          JY+=cp[k]*m1*JZ;
          }
     JX+=cn[j]*l1*JY;
      }		
J=J+cm[i]*k1*JX;	
}
J=h1*J;
return J;
}
public static double I(fdi4 f,int n)
{return I(f,n,n,n,n);}
public static double I(fdi4 f)
{return I(f,60,60,60,60);}
  public static void main(String args[]) throws IOException
  {
  double h=1.0;
  double R1=1.0;
  double R2=1.0;
  f1 ff=new f1(h,R1,R2);	  
  System.out.println("FA1-A2= "+I(ff,100)+" numeric four dimensional gauss-legendre integral");
  System.out.println("FA1-A2= "+ff.F12(h,R1,R2)+ " analytic solution ");
  }
}