// check_cylinder_deriv.c  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

#include <stdio.h>
#include <math.h>
#include "nuderiv.h"
#undef  abs
#define abs(a) ((a)<0.0?(-(a)):(a))
#undef  min
#define min(a,b) ((a)<(b)?(a):(b))
#undef  max
#define max(a,b) ((a)>(b)?(a):(b))

static const int Nr = 5; // number of radius 
static const int Nt = 7; // number of angles
static const int Nz = 5; // number of Z coordinates

  // solution for boundary and checking
static 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
  static double Ur(double R, double T, double Z)
  {
    return 3.0*R*R;
  } // end Ur

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

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

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

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

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

  static 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

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

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

int main(int argc, char * argv[]) // 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[Nr]; // r grid values 0..Nr-1
  double Tg[Nt]; // theta grid values
  double Zg[Nz]; // z grid values
  const int Nrtz = Nr*Nt*Nz;
  double Ua[Nrtz]; // U is solution, Ua solution actual for checking
  double Cr[Nr];  // coefficients of first deriv wrt r
  double Crr[Nr]; // coefficients of second deriv wrt r
  double Ctt[Nt]; // coefficients of second deriv wrt t
  double Czz[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;
  int i, j, k, ii, jj, kk;

  printf("check_cylinder_deriv.c starting\n");
  printf("r gird values, not necessarily equally spaced\n");
  for(i=0; i<Nr; i++)
  {
    Rg[i] = Rmin + Hr*(double)(i);
    printf("Rg[%1d]=%f\n", i, Rg[i]);
  }
  printf("theta gird values, not equally spaced\n");
  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(j=0; j<Nt; j++)
  {
    printf("Tg[%1d]=%f\n", j, Tg[j]);
  }
  printf("z gird values, not necessarily equally spaced\n");
  for(k=0; k<Nz; k++)
  {
    Zg[k] = Zmin+ Hz*(double)(k);
    printf("Tg[%1d]=%f\n", k, Zg[k]);
  }

  printf("analytic partial second derivative\n");
  for(i=0; i<Nr; i++)
  {
    R = Rg[i];
    for(j=0; j<Nt; j++)
    {
      T = Tg[j];
      for(k=0; k<Nz; k++)
      {
        Z = Zg[k];
        Ua[S(i,j,k)] = U(Rg[i],Tg[j],Zg[k]);
        printf("Ua[%1d,%1d,%1d]=%f, D2=%e\n",
	       i, j, k, Ua[S(i,j,k)], Uuu(Rg[i],Tg[j],Zg[k]));
        X = R*cos(T);
        Y = R*sin(T);
	if(abs(sqrt(X*X+Y*Y)-R)>0.0001)
           printf("R error\n");
        
        Tmp = atan2(Y,X);
        if(Tmp<0.0)  Tmp = 2.0*Pi + Tmp;
        if(abs(Tmp-T)>0.0001)
           printf("Theta error\n");
        
        // Z = Z
      } // k
    } // j
  } // i

  for(i=0; i<Nr; i++)
  {
    R = Rg[i];
    for(j=0; j<Nt; j++)
    {
      T = Tg[j];
      for(k=0; k<Nz; k++)
      {
        Z = Zg[k];
        // compute Ur, Urr, Utt, Uzz and call D2
        nuderiv(1, Nr, i, Rg, Cr);
        nuderiv(2, Nr, i, Rg, Crr);
        Urv = 0.0;
        Urrv = 0.0;
        for(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;
        nuderiv(2, Nt, j, Tg, Ctt);
        for(jj=0; jj<Nt; jj++)
	{
           Uttv = Uttv + Ctt[jj]*U(R,Tg[jj],Z);
        }
        Uzzv = 0.0;
        nuderiv(2, Nz, k, Zg, Czz);
        for(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 = max(Maxerr, abs(Err));
        printf("i=%1d, j=%1d, k=%1d \n", i, j, k);
        printf("Ana=%e, Cmp=%e, Err=%e \n",
               Ana, Cmp, Err);
        Err = Ur(R,T,Z)-Urv;
        Maxerrur = max(Maxerrur, abs(Err));
        printf("Ur=%e, Urv=%e, Err=%e \n",
	       Ur(R,T,Z), Urv, Err);
        Err = Urr(R,T,Z)-Urrv;
        Maxerrurr = max(Maxerrurr, abs(Err));
        printf("Urr=%e, Urrv=%e, Err=%e \n",
               Urr(R,T,Z), Urrv, Err);
        Err = Utt(R,T,Z)-Uttv;
        Maxerrutt = max(Maxerrutt, abs(Err));
        printf("Utt=%e, Uttv=%e, Err=%e \n",
	       Utt(R,T,Z), Uttv, Err);
        Err = Uzz(R,T,Z)-Uzzv;
        Maxerruzz = max(Maxerruzz, abs(Err));
        printf("Uzz=%e, Uzzv=%e, Err=%e \n",
	       Uzz(R,T,Z), Uzzv, Err);
        printf(" \n");
      } // k
    } // j
  } // i
  printf("Maxerr Ur =%e\n", Maxerrur);
  printf("Maxerr Urr=%e\n", Maxerrurr);
  printf("Maxerr Utt=%e\n", Maxerrutt);
  printf("Maxerr Uzz=%e\n", Maxerruzz);
  printf("Maxerr deriv=%e\n", Maxerr);
  printf("check_cylinder_deriv.c finishing\n");
  return 0;
} // end check_cylinder_deriv.c