// check_sphere_deriv.java  sphereical coordinates
// x = r cos(theta) sin(phi)    0 < theta < 2 Pi
// y = r sin(theta) sin(phi)    0 < phi < Pi restricted for derivatives
// z = r cos(phi)               r < 0        restricted for derivatives
// r = sqrt(x^2+y^2+z^2)
// theta = arctan(y,x) if theta<0 then theta = 2 Pi - theta
// phi = arctan(sqrt(x^2+y^2),z)
//
// compute del U(radius,theta,phi) r=radius, t=theta angle, p=phi angle
// del U  =  (d/dr U, 1/(r*sin(p)) d/dt U, (1/r) d/dp U)
//
// compute del^2 U(radius,theta,phi) r=radius, t=theta angle, p=phi angle
// del^2 U  = (1/r^2) d/dr(r^2 d/dr U) + 1/(r^2 sin(p)^2) d^2/dt^2 U +
//            (1/(r^2 sin(p) d/dp(sin(p) d/dp U)
//          = (2/r) d/dr U + d^2/dr^2 U + (1/(r^2 sin(p)^2) d^2/dt^2 U +
//            (cos(p)/(r^2 sin(t)) d/dp U + (1/(r^2)) d^2/dp^2 U

import java.text.*;
import java.io.*;

// uses nuderiv.java

public class check_sphere_deriv
{
  double Pi = 3.1415926535897932384626433832795028841971;
  final int Nr = 5;
  double Rmin = 0.3;
  double Rmax = 1.0;
  final int Nt = 7;
  double Tmin = 0.1;
  double Tmax = 300.0*Pi/180.0;
  final int Np = 5;
  double Pmin = 0.2;
  double Pmax = Pi-0.5;
  double Rg[] = new double[Nr]; // r grid values
  double Tg[] = new double[Nt]; // theta grid values
  double Pg[] = new double[Np]; // p grid values
  final int Nrtp = Nr*Nt*Np;
  double Ua[] = new double[Nrtp];
  double Cr[] = new double[Nr]; // first deriv wrt r
  double Crr[] = new double[Nr]; // second deriv wrt r
  double Ct[] = new double[Nt]; // first deriv wrt t
  double Ctt[] = new double[Nt]; // second deriv wrt t
  double Cp[] = new double[Np]; // first deriv wrt p
  double Cpp[] = new double[Np]; // second deriv wrt p
  double R, T, P, Ana, Cmp, Err;
  double Urv, Urrv, Utv, Uttv, Upv, Uppv; // computed derivatives
  double x, y, z, rtr, ttr, ptr;
  double Maxerrur=0.0, Maxerrurr=0.0, Maxerrut=0.0, Maxerrutt=0.0;
  double Maxerrup=0.0, Maxerrupp=0.0, Maxerr=0.0, Maxerrlap=0.0;

  check_sphere_deriv()
  {
    System.out.println("check_sphere_deriv.java starting");
    System.out.println("r gird values, not equally spaced");
    Rg[0] = 0.3;
    Rg[1] = 0.5;
    Rg[2] = 0.7;
    Rg[3] = 0.9;
    Rg[4] = 1.0;
    for(int I=0; I<Nr; I++)
       System.out.println("Rg["+(I)+"]="+(Rg[I]));
    System.out.println("theta gird values, not equally spaced");
    Tg[0] = 0.1;
    Tg[1] = 40.0*Pi/180.0;
    Tg[2] = 85.0*Pi/180.0;
    Tg[3] = 130.0*Pi/180.0;
    Tg[4] = 185.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("phi gird values, not equally spaced");
    Pg[0] = 0.2;
    Pg[1] = 0.9;
    Pg[2] = 1.7;
    Pg[3] = 2.4;
    Pg[4] = Pi-0.5;
    for(int K=0; K<Np; K++)
       System.out.println("Pg["+(K)+"]="+(Pg[K]));

    System.out.println("analytic partial second derivative");
    for(int I=0; I<Nr; I++)
      for(int J=0; J<Nt; J++)
	for(int K=0; K<Np; K++)
	{
          Ua[S(I,J,K)] = U(Rg[I],Tg[J],Pg[K]);
          System.out.println("Ua["+(I)+","+(J)+","+(K)+
                             "]="+(Ua[S(I,J,K)])+
                             ", D2="+(Uuu(Rg[I],Tg[J],Pg[K])));
	  x = Rg[I]*Math.cos(Tg[J])*Math.sin(Pg[K]);
	  y = Rg[I]*Math.sin(Tg[J])*Math.sin(Pg[K]);
          z = Rg[I]*Math.cos(Pg[K]);
          if(I==0 && J==0 && K==0) 
            System.out.println("x="+x+", y="+y+", z="+z);
          rtr = Math.sqrt(x*x+y*y+z*z);
          ttr = Math.atan2(y, x);
          if(ttr<0.0) ttr = 2.0*Pi + ttr; // atan2 -Pi to Pi fix
          ptr = Math.atan2(Math.sqrt(x*x+y*y),z);
          Err = Math.abs(rtr-Rg[I]);
          if(Err>0.001) 
            System.out.println("R="+Rg[I]+", rtr="+rtr+", err="+Err);
          Err = Math.abs(ttr-Tg[J]);
          if(Err>0.001) 
	    {System.out.println("T="+Tg[J]+", ttr="+ttr+", err="+Err);
	     System.out.println("x="+x+", y="+y+",z="+z);}
          Err = Math.abs(ptr-Pg[K]);
          if(Err>0.001) 
            System.out.println("P="+Pg[K]+", ptr="+ptr+", err="+Err);
          Err = Math.abs(U(Rg[I],Tg[J],Pg[K])-Uxyz(x,y,z));
          if(Err>0.001) 
	    {System.out.println("U="+U(Rg[I],Tg[J],Pg[K])+
                                ", Uxyz="+Uxyz(x,y,z)+", err="+Err);
	     System.out.println("x="+x+", y="+y+",z="+z);}

	}

    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<Np; K++)
	{
          P = Pg[K];
          // compute Ur, Urr, Utt, Upp 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,P);
            Urrv = Urrv + Crr[Ii]*U(Rg[Ii],T,P);
          }
          Utv = 0.0;
          Uttv = 0.0;
          new nuderiv(1, Nt, J, Tg, Ct);
          new nuderiv(2, Nt, J, Tg, Ctt);
          for(int Jj=0; Jj<Nt; Jj++)
	  {
            Utv = Utv + Ct[Jj]*U(R,Tg[Jj],P);
            Uttv = Uttv + Ctt[Jj]*U(R,Tg[Jj],P);
          }
          Upv = 0.0;
          Uppv = 0.0;
          new nuderiv(1, Np, K, Pg, Cp);
          new nuderiv(2, Np, K, Pg, Cpp);
          for(int Kk=0; Kk<Np; Kk++)
	  {
            Upv = Upv + Cp[Kk]*U(R,T,Pg[Kk]);
            Uppv = Uppv + Cpp[Kk]*U(R,T,Pg[Kk]);
          }

          Cmp = D2(Urv, Urrv, Uttv, Upv, Uppv, R);
          Ana = Uuu(R,T,P);
          Err = Ana - Cmp;
          Maxerr = Math.max(Maxerr, Math.abs(Err));
          System.out.println("i="+(I)+
                  ", j="+(J)+
                  ", k="+(K)+
                  ", Ana="+(Ana)+
                  ", Cmp="+(Cmp)+
                  ", err="+(Err));
          Err = Ur(R,T,P)-Urv;
          Maxerrur = Math.max(Maxerrur, Math.abs(Err));
          System.out.println("i="+(I)+", j="+(J)+", k="+(K)+
                  ", Ur="+(Ur(R,T,P))+
                  ", Urv="+(Urv)+
                  ", err="+(Err));
          Err = Urr(R,T,P)-Urrv;
          Maxerrurr = Math.max(Maxerrurr, Math.abs(Err));
          System.out.println("i="+(I)+", j="+(J)+", k="+(K)+
                  ", Urr="+(Urr(R,T,P))+
                  ", Urrv="+(Urrv)+
                  ", err="+(Err));
          Err = Ut(R,T,P)-Utv;
          Maxerrut = Math.max(Maxerrut, Math.abs(Err));
          System.out.println("i="+(I)+", j="+(J)+", k="+(K)+
                  ", Ut="+(Ut(R,T,P))+
                  ", Utv="+(Utv)+
                  ", err="+(Err));
          Err = Utt(R,T,P)-Uttv;
          Maxerrutt = Math.max(Maxerrutt, Math.abs(Err));
          System.out.println("i="+(I)+", j="+(J)+", k="+(K)+
                  ", Utt="+(Utt(R,T,P))+
                  ", Uttv="+(Uttv)+
                  ", err="+(Err));
          Err = Up(R,T,P)-Upv;
          Maxerrupp = Math.max(Maxerrupp, Math.abs(Err));
          System.out.println("i="+(I)+", j="+(J)+", k="+(K)+
                  ", Up="+(Upp(R,T,P))+
                  ", Upv="+(Upv)+
                  ", err="+(Err));
          System.out.println(" ");
          Err = Upp(R,T,P)-Uppv;
          Maxerrupp = Math.max(Maxerrupp, Math.abs(Err));
          System.out.println("i="+(I)+", j="+(J)+", k="+(K)+
                  ", Upp="+(Upp(R,T,P))+
                  ", Uppv="+(Uppv)+
                  ", err="+(Err));
          System.out.println(" ");
        }
      }
    }
    System.out.println("Maxerr Ur="+(Maxerrur));
    System.out.println("Maxerr Urr="+(Maxerrurr));
    System.out.println("Maxerr Ut="+(Maxerrut));
    System.out.println("Maxerr Utt="+(Maxerrutt));
    System.out.println("Maxerr Up="+(Maxerrup));
    System.out.println("Maxerr Upp="+(Maxerrupp));
    System.out.println("Maxerr deriv="+(Maxerr));

    System.out.println("check_sphere_deriv.java finishing");
  } //end check_sphere_deriv

  double Uxyz(double x, double y, double z)
  {
    double rtr, ttr, ptr;
    rtr = Math.sqrt(x*x+y*y+z*z);
    ttr = Math.atan2(y, x);
    if(ttr<0.0) ttr = 2.0*Pi + ttr; // atan2 -Pi to Pi fix
    ptr = Math.atan2(Math.sqrt(x*x+y*y),z);
    return rtr*rtr*rtr + ttr*ttr*ttr + ptr*ptr*ptr;
  } // end Uxyz

  double U(double R, double T, double P)
  {
    return R*R*R + T*T*T + P*P*P;
  } // end U

  double Ur(double R, double T, double P)
  {
    return 3.0*R*R;
  } // end Ur

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

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

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

  double Up(double R, double T, double P)
  {
    return 3.0*P*P;
  } // end Up

  double Upp(double R, double T, double P)
  {
    return 6.0*P;
  } // end Upp

  double Uuu(double R, double T, double P)
  {
    return 2.0*Ur(R,T,P)/R + Urr(R,T,P) + 
	   Utt(R,T,P)/(R*R*Math.sin(P)*Math.sin(P)) +
           Math.cos(P)*Up(R,T,P)/(R*R*Math.sin(P)) + 
	   Upp(R,T,P)/(R*R);
  } // end Uuu

  double D2(double Urv, double Urrv, double Uttv, double Upv,
            double Uppv, double R)
  {
    return 2.0*Urv/R + Urrv + 
	   Uttv/(R*R*Math.sin(P)*Math.sin(P)) + 
	   Math.cos(P)*Upv/(R*R*Math.sin(P)) + 
           Uppv/(R*R);
  } // end D2

  int S(int I, int J, int K)
  { // index of each grid point
    return I*Nt*Np + J*Np + K;
  } // end S

  public static void main (String[] args)
  {
    new check_sphere_deriv();
  }
} // end check_sphere_deriv.java