// pde_sphere.c  given sphere boundary values, numerically solve PDE
// r is radius,  t is theta angle from X axis,  p is phi angle from Z axis
// spherical coordinates  x = r*cos(t)*sin(p)  r = sqrt(x*x+y*y+z*z)
//                        y = r*sin(t)*sin(p)  t = atan2(y,x)
//                        z = r*cos(p)         p = atan(sqrt(x*x+y*y)/z)
// r > 0     0 <= theta < 2Pi (2Pi==0)    0 < phi < Pi  not zero or Pi
//                                        infinite del and laplacian
//
//                  ^         ^                    ^
//   del U(r,t,p) = r dU/dr + t 1/r*sin(p) dU/dt + p 1/r dU/dp
//
// laplacian del^2 U(r,t,p) = 1/(r*r) d/dr(r*r dU/dr) +
//                            1/(r*r*sin(p)*sin(p)) d^2U/dt^2
//                            1/(r*r*sin(p)) d/dp(sin(p) dU/dp) +
// simplify:                = (2/r) dU/dr + d^2U/dr^2 +
//                            (1/(r*r sin(p)*sin(p))) d^2U/dt^2 +
//                            (cos(p)/(r*r sin(p))) dU/dp + 
//                            (1/(r*r)) d^2U/dp^2

#include "deriv.h"
#include "simeq.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.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  double Pi = 3.1415926535897932384626433832795028841971;

// number of radius 
#define Nr 5
// number of theta angles
#define Nt 6
// number of phi angles                                            
#define Np 7
// all points
#define Nrtp Nr*Nt*Np
// inside boundary
#define nrow (Nr-2)*(Nt-2)*(Np-2)
#define ncol  nrow+1
#define debug 1

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 ut[nrow][ncol]; // simultaneous equation matrix
static  double Us[nrow];       // solution values inside boundary
static  double Rg[Nr]; // r grid values
static  double Tg[Nt]; // theta grid values
static  double Pg[Np]; // p grid values

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

static int S(int i, int j, int k)  // for zero based matrix i=1..Nr-2
                                   // j=1..Nt-2  k=1..Np-2
{
  if(i<1 || i>Nr-2 || j<1 || j>Nt-2 || k<1 || k>Np-2)
    printf("BAD S i=%2d, j=%2d, k=%2d \n", i, j, k);
  return (i-1)*(Nt-2)*(Np-2)+(j-1)*(Np-2)+(k-1); // index into us and ut 
} // end S

static double xval(double r, double theta, double phi)
{
  return r*cos(theta)*sin(phi);
} // end xval

static double yval(double r, double theta, double phi)
{
  return r*sin(theta)*sin(phi);
} // end yval

static double zval(double r, double theta, double phi)
{
  return r*cos(phi);
} // end zval

static double rval(double x, double y, double z)
{
  return sqrt(x*x+y*y+z*z);
} // end rval

static double thetaval(double x, double y, double z)
{
  double t = atan2(y,x);
  if(t<0.0) t = 2.0*Pi+t;
  return t; 
} // end thetaval

static double phival(double x, double y, double z)
{
  double p = atan(sqrt(x*x+y*y)/z);
  if(p<0.0) p = Pi+p;
  return p;
} // end phival

static void  checkxyzrtp(double Rg[], double Tg[], double Pg[])
{
  double x, y, z, x2, y2, z2;
  double R, T, P, r, theta, phi;
  double maxxerr, maxyerr, maxzerr, maxrerr, maxterr, maxperr;
  int i, j, k;

  printf("check xval, yval, zval, rval, thetaval, phival \n");
  maxxerr = 0.0;
  maxyerr = 0.0;
  maxzerr = 0.0;
  maxrerr = 0.0;
  maxterr = 0.0;
  maxperr = 0.0;
  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];
        // printf("r=%f, theta=%f, phi=%f \n", R, T, P);
        x = xval(R, T, P);
        y = yval(R, T, P);
	z = zval(R, T, P);
        // printf("x=%f, y=%f, z=%f \n", x, y, z);
        r = rval(x, y, z);
        theta = thetaval(x, y, z);
        phi = phival(x, y, z);
        // printf("r=%f, theta=%f, phi=%f \n", r, theta, phi);
        x2 = xval(r, theta, phi);
        y2 = yval(r, theta, phi);
        z2 = zval(r, theta, phi);
        // printf("x=%f, y=%f, z=%f \n", x2, y2, z2);
        // printf(" \n");
        maxxerr = max(maxxerr,abs(x-x2));
        maxyerr = max(maxyerr,abs(y-y2));
        maxzerr = max(maxzerr,abs(z-z2));
        maxrerr = max(maxrerr,abs(R-r));
        maxterr = max(maxterr,abs(T-theta));
        maxperr = max(maxperr,abs(P-phi));
      }
    }
  }
  printf("maxxerr=%e \n", maxxerr);
  printf("maxyerr=%e \n", maxyerr);
  printf("maxzerr=%e \n", maxzerr);
  printf("maxrerr=%e \n", maxrerr);
  printf("maxterr=%e \n", maxterr);
  printf("maxperr=%e \n", maxperr);
  printf(" \n");
} // end checkxyzrtp


static double Ub(double R, double T, double P) // boundary for testing
{
  return R*R*R + T*T*T + P*P*P;
  // future R*R*R*exp(cos(T))*exp(sin(P)); see pde_sphere2_mws.out
} // end U

static double f(double R, double T, double P) // RHS
{
  return 12.0*R + 6.0*T/(R*R*sin(P)*sin(P)) +
    3.0*cos(P)*P*P/(R*R*sin(P)) + 6.0*P/(R*R);
} // end f

static double Ur(double R, double T, double P) // analytic derivatives for test
{
  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;
} // end Ut

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

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

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

static double Uuu(double R, double T, double P) // Laplacian
{
  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

// (2/r)*ur(r,t,p) + urr(r,t,p) +
// (1/(r*r*sin(p)*sin(p)))*utt(r,t,p) +   not p==0, p==Pi
// (cos(p)/(r*r*sin(p)))*up(r,t,p) + 
// (1/(r*r))*upp(r,t,p);

static double D2(double Urv, double Urrv, double Uttv, double Upv,
		 double Uppv, double R, double P)
{  // same as  Uuu  given values of derivatives from  nuderiv
  return 2.0*Urv/R + Urrv + 
         Uttv/(R*R*sin(P)*sin(P)) +  // not p==0, p==Pi
	 cos(P)*Upv/(R*R*sin(P)) + 
         Uppv/(R*R);
} // end D2

static void check_deriv(double Rg[], double Tg[], double Pg[])
{
  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;
  int i, j, k, ii, jj, kk;

  printf("check_deriv Ur, Urr, Ut, Utt, Up, Upp, Uuu, D2 \n");
  printf("check analytic against numerically computer \n");
  printf("Uuu is Laplacian in spherical coordinates, D2 numeric \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];
        // 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);
        }
        Utv = 0.0;
        Uttv = 0.0;
        nuderiv(1, Nt, j, Tg, Ct);
        nuderiv(2, Nt, j, Tg, Ctt);
        for(jj=0; jj<Nt; jj++)
	{
          Utv = Utv + Ct[jj]*Ub(R,Tg[jj],P);
          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, abs(Err));
        Err = Ur(R,T,P)-Urv;
        Maxerrur = max(Maxerrur, abs(Err));
        Err = Urr(R,T,P)-Urrv;
        Maxerrurr = max(Maxerrurr, abs(Err));
        Err = Ut(R,T,P)-Utv;
        Maxerrut = max(Maxerrut, abs(Err));
        Err = Utt(R,T,P)-Uttv;
        Maxerrutt = max(Maxerrutt, abs(Err));
        Err = Up(R,T,P)-Upv;
        Maxerrup = max(Maxerrup, abs(Err));
        Err = Upp(R,T,P)-Uppv;
        Maxerrupp = max(Maxerrupp, abs(Err));
      }
    }
  }
  printf("Maxerr Ur=%e\n", Maxerrur);
  printf("Maxerr Urr=%e\n", Maxerrurr);
  printf("Maxerr Ut=%e\n", Maxerrut);
  printf("Maxerr Utt=%e\n", Maxerrutt);
  printf("Maxerr Up=%e\n", Maxerrup);
  printf("Maxerr Upp=%e\n", Maxerrupp );
  printf("Maxerr deriv=%e\n", Maxerr );
} // end check_deriv

// write_soln  for plot3D 
static void  write_soln(double Xa[], double Ya[], double Za[],
                        double Ua[], char file_name[])
{
  FILE * Fout;
  int i, j, k;

  Fout = fopen(file_name, "w");
  if(Fout==NULL)
  {
    printf("ERROR unable to open %s for writing \n", file_name);
    exit(1);
  }
  printf("writing %s\n", file_name);
  fprintf(Fout,"%d %d %d %d \n", 3, Nr, Nt, Np);
  for(i=0; i<Nr; i++)
  {
    for(j=0; j<Nt; j++)
    {
      for(k=0; k<Np; k++)
      {
        fprintf(Fout,"%f %f %f %f \n",
                Xa[S0(i,j,k)],Ya[S0(i,j,k)],Za[S0(i,j,k)],Ua[S0(i,j,k)]);
                // x,y,z,value 
      } // k 
    } // j 
  } // i 
  fclose(Fout);
  printf("finished writing %s\n", file_name);
}  // end write_soln 


static void print_matrix()
{
  int i, j, k, ii, jj, kk, irow, jcol, cs;
  printf("matrix = \n");
  cs = (Nr-2)*(Nt-2)*(Np-2);
  for(i=1; i<Nr-1; i++)  // inside boundary
  {
    for(j=1; j<Nt-1; j++)
    {
      for(k=1; k<Np-1; k++)
      {
        irow = S(i,j,k);
        for(ii=1; ii<Nr-1; ii++)
        {
          for(jj=1; jj<Nt-1; jj++)
          {
            for(kk=1; kk<Np-1; kk++)
            {
              jcol = S(ii,jj,kk);
              printf("Ut(%2d,%2d,%2d) = %f \n", i, j, k, ut[irow][jcol]);
	    } // kk
	  } // jj
        } // ii
        printf("RHS(%2d) = %f \n", irow, ut[irow,cs]);
      } // k      
    } // j
  } // i
  printf(" \n");
} // end print_matrix

static void init_matrix() // based on PDE 
{
  int i, j, k, m, ii, jj, kk, cs;
  int irow, jcol;
  double coefr, coeft, coeftt, coefp, coefpp;
  double R, T, P;

  cs = nrow; // constant column subscript 
  printf("Nr=%d, Nt=%d, Np=%d, nrow=%d, ncol=%d \n", Nr, Nt, Np, nrow, ncol);
  printf("S(1,1,1)=%d, S(Nr-2,Nt-2,Np-2)=%d \n", S(1,1,1), S(Nr-2,Nt-2,Np-2));
  // zero internal cells (zero based row and column, S(i,j) )
  for(i=1; i<Nr-1; i++)
  {
    for(j=1; j<Nt-1; j++)
    {
      for(k=1; k<Np-1; k++)
      {
        irow = S(i,j,k);
        for(ii=1; ii<Nr-1; ii++)
        {
          for(jj=1; jj<Nt-1; jj++)
	  {
            for(kk=1; kk<Np-1; kk++)
	    {
              jcol = S(ii,jj,kk);
              ut[irow][jcol] = 0.0;
	    } // kk
	  } // jj
        } // ii
        ut[irow][cs] = 0.0;
      } // k
    } // j
  } // i
  printf("matrix zeroed \n");

  for(i=1; i<Nr-1; i++) // inside boundary 
  {
    R = Rg[i];
    for(j=1; j<Nt-1; j++)                 
    {
      T = Rg[j];
      for(k=1; k<Np-1; k++)                 
      {
	P = Pg[k];
        irow = S(i,j,k); // row of matrix 
        // for each term 
        // (2/r) dU/dr 
        nuderiv(1,Nr,i,Rg,Cr);
        coefr = 2.0/R;
        for(m=0; m<Nr; m++)
        {
	  if(m==0 || m==Nr-1) ut[irow][cs] = ut[irow][cs] - 
                              coefr*Cr[m]*Ub(Rg[m],Tg[j],Pg[k]);
	  else  ut[irow][S(m,j,k)] = ut[irow][S(m,j,k)] + coefr*Cr[m];
        }

        // d^2U/dr^2
        nuderiv(2,Nr,i,Rg,Crr);
        for(m=0; m<Nr; m++)
        {
	  if(m==0 || m==Nr-1) ut[irow][cs] = 
			      ut[irow][cs] - Crr[m]*Ub(Rg[m],Tg[j],Pg[k]);
	  else  ut[irow][S(m,j,k)] = ut[irow][S(m,j,k)] + Crr[m];
        }

        // (1/(r*r sin(p)*sin(p))) d^2U/dt^2 
        nuderiv(2,Nt,j,Tg,Ctt);
        coeftt = 1.0/(R*R*sin(P)*sin(P));
        for(m=0; m<Nt; m++)
        {
	  if(m==0 || m==Nt-1) ut[irow][cs] =  ut[irow][cs] - 
                              coeftt*Ctt[m]*Ub(Rg[i],Tg[m],Pg[k]);
	  else  ut[irow][S(i,m,k)] = ut[irow][S(i,m,k)] + coeftt*Ctt[m];
        }

        // (cos(p)/(r*r sin(p))) dU/dp 
        nuderiv(1,Np,k,Pg,Cp);
	coefp = cos(P)/(R*R*sin(P));
        for(m=0; m<Np; m++)
        {
	  if(m==0 || m==Np-1) ut[irow][cs] = ut[irow][cs] -
                              coefp*Cp[m]*Ub(Rg[i],Tg[j],Pg[m]);
	  else  ut[irow][S(i,j,m)] = ut[irow][S(i,j,m)] + coefp*Cp[m];
        }

        // (1/(r*r)) d^2U/dp^2
        nuderiv(2,Np,k,Pg,Cpp);
	coefpp = 1.0/(R*R);
        for(m=0; m<Np; m++)
        {
	  if(m==0 || m==Np-1) ut[irow][cs] = ut[irow][cs] - 
                              coefpp*Cpp[m]*Ub(Rg[i],Tg[j],Pg[m]);
	  else  ut[irow][S(i,j,m)] = ut[irow][S(i,j,m)] + coefpp*Cpp[m];
        }

      // RHS constant 
	ut[irow][cs] = ut[irow][cs] + f(Rg[i],Tg[j],Pg[k]);
      // end terms
      } // k 
    } // j
  } // i
} // end init_matrix 

int main(int argc, char * argv[])
{
  double Ua[Nrtp];
  double Xa[Nrtp]; // for x,y,z plotting
  double Ya[Nrtp];
  double Za[Nrtp];
  int i, j, k, ii, jj, kk;
  double R, T, P;

  printf("pde_sphere.c running \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[%1d]=%f\n", i, Rg[i]);
  printf("theta gird values, not equally spaced\n");
  Tg[0] = Pi*  5.0/180.0; // middle number is degrees
  Tg[1] = Pi* 60.0/180.0;
  Tg[2] = Pi*125.0/180.0;
  Tg[3] = Pi*180.0/180.0;
  Tg[4] = Pi*265.0/180.0;
  Tg[5] = Pi*330.0/180.0;
  for(j=0; j<Nt; j++)
    printf("Tg[%1d]=%f\n", j, Tg[j]);
  printf("phi gird values, not equally spaced\n");
  Pg[0] = Pi*15.0/180.0;
  Pg[1] = Pi*40.0/180.0;
  Pg[2] = Pi*60.0/180.0;
  Pg[3] = Pi*90.0/180.0;
  Pg[4] = Pi*110.0/180.0;
  Pg[5] = Pi*140.0/180.0;
  Pg[6] = Pi*165.0/180.0;
  for(k=0; k<Np; k++)
    printf("Pg[%1d]=%f\n", k, Pg[k]);

  if(debug) checkxyzrtp(Rg, Tg, Pg);
  if(debug) check_deriv(Rg, Tg, Pg);

  printf("S0(i,j,k) 0,0,0=%d, max=%d \n", S0(0,0,0), S0(Nr-1,Nt-1,Np-1));
  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[S0(i,j,k)] = Uuu(R, T, P);
        Xa[S0(i,j,k)] = xval(R, T, P);
        Ya[S0(i,j,k)] = yval(R, T, P);
        Za[S0(i,j,k)] = zval(R, T, P);
      }
    }
  }
  printf("Ua has analytic solution \n");
  write_soln(Xa, Ya, Za, Ua, "pde_sphere.dat");

  init_matrix();
  // print_matrix();  // debug print, big
  simeq2b(nrow, ut, Us);
  printf("                      analytic \n");
  for(i=0; i<Nr; i++)
  {
    for(j=0; j<Nt; j++)
    {
      for(k=0; k<Np; k++)
      {
        if(i<1 || i>Nr-2 || j<1 || j>Nt-2 || k<1 || k>Np-2)
          printf("Ub(%d,%d,%d) =              %4.4f \n",
                 i, j, k, Ub(Rg[i],Tg[j],Pg[k]));
        else
          printf(" U(%d,%d,%d) = %4.4f     %4.4f \n",
                 i, j, k, Us[S(i,j,k)], Ub(Rg[i],Tg[j],Pg[k]));
      } // k
      printf(" \n");
    } // j
    printf(" \n");
  } // i

  printf("pde_sphere.c finished \n");
  return 0;
} // end pde_sphere.c