// check_cylinder_deriv.java  cylindrical coordinates
// x = r cos theta    0 <= theta < 2 Pi
// y = r sin theta    r < 0  restriction for derivative
// z = z              z double
// r = sqrt(x^2+y^2)
// theta = arctan(y,x)  if theta<0  theta = 2 Pi + arctan(y,x)
// z = z
//
// compute del U(radius,theta,z) r=radius, t=theta angle, z=height
// del U  =  (d/dr U, 1/r d/dt U, d/dz U)
//
// compute del^2 U(radius,theta,z) r=radius, t=theta angle, z=height
// del^2 U  =  d/dr(r d/dr U) + 1/r^2 d^2/dt^2 U + d^2/dz^2 U
//          =  2/r d/dr U + d^2/dr^2 U + 1/r^2 d^2/dt^2 U + d^2/dz^2 U

// needs nuderiv.java


public class check_cylinder_deriv
{
  final int Nr = 5; // number of radius 
  final int Nt = 7; // number of angles
  final int Nz = 5; // number of Z coordinates

  public check_cylinder_deriv()
  {
  double Pi = 3.1415926535897932384626433832795028841971;
  double Rmin = 0.2;
  double Rmax = 1.0;
  double Hr = (Rmax-Rmin)/(double)(Nr-1);
  double Tmin = 0.0;
  double Tmax = 300.0*Pi/180.0;
  // Ht not applicable
  double Zmin = -1.0;
  double Zmax =  1.0;
  double Hz = (Zmax-Zmin)/(double)(Nz-1);
  double Rg[] = new double[Nr]; // r grid values 0..Nr-1
  double Tg[] = new double[Nt]; // theta grid values
  double Zg[] = new double[Nz]; // z grid values
  final int Nrtz = Nr*Nt*Nz;
  double Ua[] = new double[Nrtz]; // U is solution
                                  // Ua solution actual for checking
  double Cr[] = new double[Nr];  // coefficients of first deriv wrt r
  double Crr[] = new double[Nr]; // coefficients of second deriv wrt r
  double Ctt[] = new double[Nt]; // coefficients of second deriv wrt t
  double Czz[] = new double[Nz]; // coefficients of second deriv wrt z
  double R, T, Z, Ana, Cmp, Err;
  double Urv, Urrv, Uttv, Uzzv; // computed derivatives
  double Maxerrur = 0.0, Maxerrurr = 0.0, Maxerrutt = 0.0;
  double Maxerruzz = 0.0, Maxerr = 0.0;
  double X, Y, Tmp;

  System.out.println("check_cylinder_deriv.java starting");
  System.out.println("r gird values, not necessarily equally spaced");
  for(int i=0; i<Nr; i++)
  {
    Rg[i] = Rmin + Hr*(double)(i);
    System.out.println("Rg["+i+"]= "+Rg[i]);
  }
  System.out.println("theta gird values, not equally spaced");
  Tg[0] = 0.0;
  Tg[1] = 40.0*Pi/180.0;
  Tg[2] = 85.0*Pi/180.0;
  Tg[3] = 130.0*Pi/180.0;
  Tg[4] = 180.0*Pi/180.0;
  Tg[5] = 240.0*Pi/180.0;
  Tg[6] = 300.0*Pi/180.0;
  for(int j=0; j<Nt; j++)
  {
    System.out.println("Tg["+j+"]= "+Tg[j]);
  }
  System.out.println("z gird values, not necessarily equally spaced");
  for(int k=0; k<Nz; k++)
  {
    Zg[k] = Zmin+ Hz*(double)(k);
    System.out.println("Tg["+k+"]= "+ Zg[k]);
  }

  System.out.println("analytic partial second derivative");
  for(int i=0; i<Nr; i++)
  {
    R = Rg[i];
    for(int j=0; j<Nt; j++)
    {
      T = Tg[j];
      for(int k=0; k<Nz; k++)
      {
        Z = Zg[k];
        Ua[S(i,j,k)] = U(Rg[i],Tg[j],Zg[k]);
        System.out.println("Ua["+i+","+j+","+k+"]="+
	       Ua[S(i,j,k)]+", D2= "+Uuu(Rg[i],Tg[j],Zg[k]));
        X = R*Math.cos(T);
        Y = R*Math.sin(T);
	if(Math.abs(Math.sqrt(X*X+Y*Y)-R)>0.0001)
           System.out.println("R error");
        
        Tmp = Math.atan2(Y,X);
        if(Tmp<0.0)  Tmp = 2.0*Pi + Tmp;
        if(Math.abs(Tmp-T)>0.0001)
           System.out.println("Theta error");
        
        // Z = Z
      } // k
    } // j
  } // i

  for(int i=0; i<Nr; i++)
  {
    R = Rg[i];
    for(int j=0; j<Nt; j++)
    {
      T = Tg[j];
      for(int k=0; k<Nz; k++)
      {
        Z = Zg[k];
        // compute Ur, Urr, Utt, Uzz and call D2
        new nuderiv(1, Nr, i, Rg, Cr);
        new nuderiv(2, Nr, i, Rg, Crr);
        Urv = 0.0;
        Urrv = 0.0;
        for(int ii=0; ii<Nr; ii++)
	{
           Urv = Urv + Cr[ii]*U(Rg[ii],T,Z);
           Urrv = Urrv + Crr[ii]*U(Rg[ii],T,Z);
	}
        Uttv = 0.0;
        new nuderiv(2, Nt, j, Tg, Ctt);
        for(int jj=0; jj<Nt; jj++)
	{
           Uttv = Uttv + Ctt[jj]*U(R,Tg[jj],Z);
        }
        Uzzv = 0.0;
        new nuderiv(2, Nz, k, Zg, Czz);
        for(int kk=0; kk<Nz; kk++)
	{
           Uzzv = Uzzv + Czz[kk]*U(R,T,Zg[kk]);
        }

        Cmp = D2(Urv, Urrv, Uttv, Uzzv, R);
        Ana = Uuu(R,T,Z);
        Err = Ana - Cmp;
        Maxerr = Math.max(Maxerr, Math.abs(Err));
        System.out.println("i="+i+", j="+j+", k="+k);
        System.out.println("Ana="+Ana+", Cmp="+Cmp+", Err="+Err);
        Err = Ur(R,T,Z)-Urv;
        Maxerrur = Math.max(Maxerrur, Math.abs(Err));
        System.out.println("Ur="+Ur(R,T,Z)+", Uvr="+Urv+", Err="+Err);
        Err = Urr(R,T,Z)-Urrv;
        Maxerrurr = Math.max(Maxerrurr, Math.abs(Err));
        System.out.println("Urr="+Urr(R,T,Z)+", Urrv="+Urrv+", Err="+Err);
        Err = Utt(R,T,Z)-Uttv;
        Maxerrutt = Math.max(Maxerrutt, Math.abs(Err));
        System.out.println("Utt= "+Utt(R,T,Z)+", Uttv= "+
                           Uttv+", Err= "+Err);
        Err = Uzz(R,T,Z)-Uzzv;
        Maxerruzz = Math.max(Maxerruzz, Math.abs(Err));
        System.out.println("Uzz= "+Uzz(R,T,Z)+", Uzzv= "+
                            Uzzv+", Err= "+Err);
        System.out.println(" ");
      } // k
    } // j
  } // i
  System.out.println("Maxerr Ur = "+ Maxerrur);
  System.out.println("Maxerr Urr= "+ Maxerrurr);
  System.out.println("Maxerr Utt= "+ Maxerrutt);
  System.out.println("Maxerr Uzz= "+ Maxerruzz);
  System.out.println("Maxerr deriv= "+ Maxerr);
  System.out.println("check_cylinder_deriv.java finishing");
  } // end cylinder_deriv

  // solution for boundary and checking
  double U(double R, double T, double Z)
  {
    return R*R*R + T*T*T + Z*Z*Z;
  } // end U

  // analytic derivatives of solution, for checking
  double Ur(double R, double T, double Z)
  {
    return 3.0*R*R;
  } // end Ur

  double Urr(double R, double T, double Z)
  {
    return 6.0*R;
  } // end Urr

  double Ut(double R, double T, double Z)
  {
    return 3.0*T*T;
  } // end Ut

  double Utt(double R, double T, double Z)
  {
    return 6.0*T;
  } // end Utt

  double Uz(double R, double T, double Z)
  {
    return 3.0*Z*Z;
  } // end Uz

  double Uzz(double R, double T, double Z)
  {
    return 6.0*Z;
  } // end Uzz

  double Uuu(double R, double T, double Z)
  {
    return 2.0*Ur(R,T,Z)/R + Urr(R,T,Z) + Utt(R,T,Z)/(R*R) + Uzz(R,T,Z);
  } // end Uuu

  double D2(double Urv, double Urrv, double Uttv, double Uzzv, double R)
  {
    return 2.0*Urv/R + Urrv + Uttv/(R*R)+ Uzzv;
  } // end D2

  int S(int i, int j, int k)
  { // index of each grid point
    return i*Nt*Nz + j*Nz + k;
  } // end S

  public static void main (String[] args) // "main" required
  {
    new check_cylinder_deriv();
  }
} // end class and check_cylinder_deriv.java