/* pde_spherical_eq.c  sphereical coordinates
 * solve del^2 U(radius,theta,phi) + U(radius,theta,phi) = f(r,t,p)
 * given Dirchlet boundary values Ub(r,t,p), and f(r,t,p)
 * much checking in this code
 * analytical solution is U(r,t,p)=r^3 + t^3/216 + p^3/27
 *
 * 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 + arctan(y,x)
 * phi = arctan(sqrt(x^2+y^2),z)
 *
 * use existing partial derivative computation in Cartesian coordinates to
 * compute partial derivatives in spherical coordinates
 *
 * 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^2 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(p)) d/dp U + (1/r^2) d^2/dp^2 U
 *
 * f(r,t,p) := 12*r + (1/36)*t/(r^2*sin(p)^2) + 
 *             (1/9)*cos(p)*p^2/(r^2*sin(p)) + (2/9)*p/r^2 +
 *             r^3 + (1/216)*t^3 + (1/27)*p^3
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "simeq.h"
#include "nuderiv.h"

#undef max
#define max(a,b) ((a)>(b)?(a):(b))
#undef abs
#define abs(a) ((a)<0.0?(-(a)):(a))

#define Nr 5
#define Nt 7
#define Np 5
#define Nrtp  (Nr-2)*(Nt-2)*(Np-2) /* DOF */
static double Rg[Nr]; /* r grid values */
static double Tg[Nt]; /* theta grid values */
static double Pg[Np]; /* p grid values */
static double Ua[Nr*Nt*Np]; /* boundary and analytic solution for check */
static double Ud2[Nr*Nt*Np]; /* del^2 for check */
static double Cr[Nr]; /* first deriv wrt r */
static double Crr[Nr]; /* second deriv wrt r */
static double Ct[Nt]; /* first deriv wrt t */
static double Ctt[Nt]; /* second deriv wrt t */
static double Cp[Np]; /* first deriv wrt p */
static double Cpp[Np]; /* second deriv wrt p */
static double A[Nrtp][Nrtp+1]; /* no bound, has RHS */
static double X[Nrtp];

static double Ub(double R, double T, double P); /* boundary values */
static double Ur(double R, double T, double P);
static double Urr(double R, double T, double P);
static double Ut(double R, double T, double P);
static double Utt(double R, double T, double P);
static double Up(double R, double T, double P);
static double Upp(double R, double T, double P);
static double Uuu(double R, double T, double P);
static double D2(double Urv, double Urrv, double Uttv, double Upv,
                 double Uppv, double R, double P);
static int S(int I, int J, int K);
static int Sd(int I, int J, int K);
static double f(double r, double t, double p); /* RHS */
static void build_A();

int main(int argc, char *argv[])
{
  double Pi = 3.1415926535897932384626433832795028841971;
  double Rmin = 0.3;
  double Rmax = 1.0;
  double Tmin = 0.1;
  double Tmax = 300.0*Pi/180.0;
  double Pmin = 0.2;
  double Pmax = Pi-0.5;
  double R, T, P, Ana, Cmp, Err;
  double Urv, Urrv, Uttv, Upv, Uppv; /* computed derivatives */
  double Maxerrur=0.0, Maxerrurr=0.0, Maxerrutt=0.0;
  double Maxerrup=0.0, Maxerrupp=0.0, Maxerr=0.0;
  int I, J, K, Ii, Jj, Kk;

  printf("pde_spherical_eq.c starting\n");
  printf("r gird values, not equally spaced\n");
  Rg[0] = 0.3;
  Rg[1] = 0.5;
  Rg[2] = 0.7;
  Rg[3] = 0.9;
  Rg[4] = 1.0;
  for(I=0; I<Nr; I++)
    printf("Rg[%d]=%g\n", I, Rg[I]);
  printf("theta gird values, not equally spaced\n");
  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(J=0; J<Nt; J++)
    printf("Tg[%d]=%g\n", J, Tg[J]);
  printf("phi gird values, not equally spaced\n");
  Pg[0] = 0.2;
  Pg[1] = 0.9;
  Pg[2] = 1.7;
  Pg[3] = 2.4;
  Pg[4] = Pi-0.5;
  for(K=0; K<Np; K++)
    printf("Pg[%d]=%g\n", K, Pg[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<Np; K++)
      {
        P = Pg[K];
        Ua[S(I,J,K)] = Ub(R,T,P);
        Ud2[S(I,J,K)] = D2(Ur(R,T,P),Urr(R,T,P),Utt(R,T,P),
                           Up(R,T,P),Upp(R,T,P),R,P);
        printf("[%1d,%1d,%1d] Ua=%g, Ud2=%g\n",
               I, J, K, Ua[S(I,J,K)], Ud2[S(I,J,K)]);
        printf("Uuu+U=%g, f=%g\n", Uuu(R,T,P)+Ub(R,T,P),f(R,T,P));
      }
    }
  }
  for(I=0; I<Nr; I++)
  {
    R = Rg[I];
    for(J=0; J<Nt; J++)
    {
      T = Tg[J];
      for(K=0; K<Np; K++)
      {
        P = Pg[K];
        /* compute Ur, Urr, Utt, Upp 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]*Ub(Rg[Ii],T,P);
          Urrv = Urrv + Crr[Ii]*Ub(Rg[Ii],T,P);
        }
        Uttv = 0.0;
        nuderiv(2, Nt, J, Tg, Ctt);
        for(Jj=0; Jj<Nt; Jj++)
	{
          Uttv = Uttv + Ctt[Jj]*Ub(R,Tg[Jj],P);
        }
        Upv = 0.0;
        Uppv = 0.0;
        nuderiv(1, Np, K, Pg, Cp);
        nuderiv(2, Np, K, Pg, Cpp);
        for(Kk=0; Kk<Np; Kk++)
	{
          Upv = Upv + Cp[Kk]*Ub(R,T,Pg[Kk]);
          Uppv = Uppv + Cpp[Kk]*Ub(R,T,Pg[Kk]);
        }

        Cmp = D2(Urv, Urrv, Uttv, Upv, Uppv, R, P);
        Ana = Uuu(R,T,P);
        Err = Ana - Cmp;
        Maxerr = max(Maxerr, Err);
        Err = Ur(R,T,P)-Urv;
        Maxerrur = max(Maxerrur, Err);
        Err = Urr(R,T,P)-Urrv;
        Maxerrurr = max(Maxerrurr, Err);
        Err = Utt(R,T,P)-Uttv;
        Maxerrutt = max(Maxerrutt, Err);
        Err = Up(R,T,P)-Upv;
        Maxerrup = max(Maxerrup, Err);
        Err = Upp(R,T,P)-Uppv;
        Maxerrupp = max(Maxerrupp, Err);
      }
    }
  }
  printf("Maxerr Ur=%g\n", Maxerrur);
  printf("Maxerr Urr=%g\n", Maxerrurr);
  printf("Maxerr Utt=%g\n", Maxerrutt);
  printf("Maxerr Utp=%g\n", Maxerrup);
  printf("Maxerr Upp=%g\n", Maxerrupp);
  printf("Maxerr deriv=%g\n", Maxerr);

  build_A();
  simeq2b(Nrtp, A, X); /* X is solution */
  printf("check against analytic solution\n");
  Maxerr = 0.0;
  for(I=1; I<Nr-1; I++)
  {
    for(J=1; J<Nt-1; J++)
    {
      for(K=1; K<Np-1; K++)
      {
        Err = X[Sd(I, J, K)]-Ua[S(I,J,K)];
        Maxerr = max(Maxerr, abs(Err));
        printf("[%1d,%1d,%1d] X=%g, ana=%g, Err=%g\n",
               I, J, K, X[Sd(I, J, K)], Ua[S(I,J,K)], Err);
      } /* end K */
    } /* end J */
  } /* end I */
  printf("Maxerr soln=%g\n", Maxerr);

  printf("pde_spherical_eq.c finishing\n");
} /* end pde_spherical_eq */

static double Ub(double R, double T, double P) /* boundary values */
{
  return R*R*R + T*T*T/216.0 + P*P*P/27.0;
} /* end Ub */

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

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

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

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

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

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

static 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*sin(P)*sin(P)) +
         cos(P)*Up(R,T,P)/(R*R*sin(P)) + 
         Upp(R,T,P)/(R*R);
} /* end Uuu */

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

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

static int Sd(int I, int J, int K)
{ /* index of each grid point inside boundary, A and X */
  return (I-1)*(Nt-2)*(Np-2) + (J-1)*(Np-2) + (K-1);
} /* end S */

static double f(double r, double t, double p) /* RHS */
{
  double sp = sin(p);
  return 12.0*r + t/(36.0*r*r*sp*sp) +
         cos(p)*p*p/(9.0*r*r*sp) + 
         2.0*p/(9.0*r*r) +
         r*r*r + t*t*t/216.0 + p*p*p/27.0;
} /* end f  RHS */

static void build_A() /* A is discrete system of equations, including Y */
{              /* solution is X, from A*X=Y */
  int row, idof, cs;
  double Urrs, Utts, Upps, Ud2ck, Err;
  double Err1, Err2, Err3;
  double R, T, P;
  int I, J, K, Ii, Jj, Kk;

  /* equation Urrs(r,t,p)+Utts(r,t,p)+Upps(r,t,p)+U(r,t,p) = f(r,t,p) */
  /* for every (r,t,p) inside boundary */
  printf("building discrete system of %d equations\n", Nrtp);
  cs = Nrtp; /* RHS */
  /* initialize A and X */
  for(I=1; I<Nr-1; I++)
  {
    for(J=1; J<Nt-1; J++)
    {
      for(K=1; K<Np-1; K++)
      {
        row = Sd(I, J, K); /* equation index */
        for(Ii=1; Ii<Nr-1; Ii++)
        {
          for(Jj=1; Jj<Nt-1; Jj++)
          {
            for(Kk=1; Kk<Np-1; Kk++)
            {
              A[row][Sd(Ii,Jj,Kk)] = 0.0;
	    } /* end Kk */
          } /* end Jj */
        } /* end Ii */
        A[row][cs] = f(Rg[I],Tg[J],Pg[K]); /* RHS */
        X[row] = 0.0;
      } /* end K */
    } /* end J */
  } /* end I */
  /* now build A inside boundary (could use sparse) */
  for(I=1; I<Nr-1; I++)
  {
    R = Rg[I];
    for(J=1; J<Nt-1; J++)
    {
      T = Tg[J];
      for(K=1; K<Np-1; K++)
      {
        P = Pg[K];
        row = Sd(I, J, K); /* equation index */
        Ud2ck = 0.0; /* only for development checking */

        /* equation Urrs(r,t,p)+Utts(r,t,p)+Upps(r,t,p)+U(r,t,p) = f(r,t,p) */
        /* Urrs term */
        nuderiv(1, Nr, I, Rg, Cr);
        nuderiv(2, Nr, I, Rg, Crr);
        for(Ii=0; Ii<Nr; Ii++)
	{
          Urrs = 2.0*Cr[Ii]/R + Crr[Ii];
          Ud2ck += Urrs;

	  if(Ii==0 || Ii==Nr-1) /* boundary, negative on RHS */
	  {
	    A[row][cs] -= Urrs*Ub(Rg[Ii],T,P);
	  } 
          else
	  {
            idof = Sd(Ii,J,K);
            A[row][idof] += Urrs;
	  } /* end if */
        } /* end Ii */

          /* Utts term */ 
          nuderiv(2, Nt, J, Tg, Ctt);
          for(Jj=0; Jj<Nt; Jj++)
	  {
	    Utts = Ctt[Jj]/(R*R*sin(P)*sin(P));
            Ud2ck += Utts;
	    if(Jj==0 || Jj==Nt-1) /* boundary, negative on RHS */
	    {
	       A[row][cs] -= Utts*Ub(R,Tg[Jj],P);
	    } 
            else
	    {
              idof = Sd(I,Jj,K);
              A[row][idof] += Utts;
	    } /* end if */
	  } /* end Jj */

          /* Upps term */ 
          nuderiv(1, Np, K, Pg, Cp);
          nuderiv(2, Np, K, Pg, Cpp);
          for(Kk=0; Kk<Np; Kk++)
	  {
	    Upps = cos(P)*Cp[Kk]/(R*R*sin(P)) + Cpp[Kk]/(R*R);
            Ud2ck += Upps;
	    if(Kk==0 || Kk==Np-1) /* boundary, negative on RHS */
	    {
	       A[row][cs] -= Upps*Ub(R,T,Pg[Kk]);
	    } 
            else
	    {
              idof = Sd(I,J,Kk);
              A[row][idof] += Upps;
	    } /* end if */
	  } /* end Kk */

	  /* U term */
	  A[row][row] += 1.0;

          /* check */
          Err = abs(Ud2ck);
          if(Err> 0.001)
	    printf("[%1d,%1d,%1d] Ud2ck=%g, Err=%g\n",
                   I, J, K, Ud2ck, Err);
          if(row==2) /* debug print */
	  {
            Err1 = 0.0;
            Err2 = 0.0;
            Err3 = 0.0;
	    for(Ii=1; Ii<Nr-1; Ii++)
	      for(Jj=1; Jj<Nt-1; Jj++) 
		for(Kk=1; Kk<Np-1; Kk++)
		{
                  printf("[%1d,%1d,%1d] A[0]=%g, A[1]=%g, A[2]=%g\n",
			 Ii, Jj, Kk, A[0][Sd(Ii,Jj,Kk)],
                         A[1][Sd(Ii,Jj,Kk)],
                         A[2][Sd(Ii,Jj,Kk)]);
                  Err1 += A[0][Sd(Ii,Jj,Kk)]*Ub(Rg[Ii],Tg[Jj],Pg[Kk]);
                  Err2 += A[1][Sd(Ii,Jj,Kk)]*Ub(Rg[Ii],Tg[Jj],Pg[Kk]);
                  Err3 += A[2][Sd(Ii,Jj,Kk)]*Ub(Rg[Ii],Tg[Jj],Pg[Kk]);
		} /* Kk, Jj, Ii */
	  printf("RHS A[0][cs]=%g,  A[1][cs]=%g, A[2][cs]=%g\n",
		  A[0][cs],A[1][cs],A[2][cs]);
          Err1 -= A[0][cs];
          Err2 -= A[1][cs];
          Err3 -= A[2][cs];
          printf("row check\n");
          printf("Err1=%E,  Err2=%E, Err3=%E \n", Err1,Err2,Err3);
        } /* end if  debug print */
      } /* end K */
    } /* end J */
  } /* end I */
} /* end build_A */
/* end pde_spherical_eq.c */