/* spherical_deriv.c   check numerical calculation of derivatives   */
/*                     in spherical coordinates.                    */
/* data structure  rtp[][][] field values at radius, theta and phi  */
/*                 rval[]    values of r at index 0, 1, ..., nr-1   */
/*                 tval[]    values of theta -Pi, -Pi+dt, ... Pi-dt */
/*                 pval[]    values of phi   =pi/2, -Pi/2+dp, ,,,   */
/*                                                                  */
/* method, use nuderiv to compute numerical derivatives             */
/* u(r,t,p) = */
/* du/dr =  analytic */
/* du/dt =  analytic */
/* du/dp =  analytic */

#include <stdio.h>
#include <math.h>
#include "nuderiv.h"
#define Pi     3.14159265358979323846

#define nr 7
#define nt 16
#define np 12
static double rtp[nr][nt][np];
static double rval[nr] = { 0.25, 0.5, 0.75, 1.0 , 1.25, 1.5, 1.75};
static double tval[nt]; /* computed below */
static double pval[np];
static double rcoef[nr];
static double tcoef[nt];
static double pcoef[np];

double u(double r, double t, double p)
{
  return 1.0/((r*r+0.1) *
              ((t-Pi/4.0)*(t-Pi/4.0)+0.2) *
              ((p-Pi/8.0)*(p-Pi/8.0)+0.3) );
}

double dudr(double r, double t, double p)
{
  return -(2.0*r)/((r*r+0.1)*(r*r+0.1) *
                   ((t-Pi/4.0)*(t-Pi/4.0)+0.2) *
                   ((p-Pi/8.0)*(p-Pi/8.0)+0.3) );
 
}

double dudt(double r, double t, double p)
{
  return -(2.0*t-Pi/2.0)/((r*r+0.1) *
                  ((t-Pi/4.0)*(t-Pi/4.0)+0.2) * ((t-Pi/4.0)*(t-Pi/4.0)+0.2) *
                  ((p-Pi/8.0)*(p-Pi/8.0)+0.3));

}

double dudp(double r, double t, double p)
{
  return -(2.0*p-Pi/4.0)/((r*r+0.1) * 
                  ((t-Pi/4.0)*(t-Pi/4.0)+0.2) *
                  ((p-Pi/8.0)*(p-Pi/8.0)+0.3) * ((p-Pi/8.0)*(p-Pi/8.0)+0.3));

}

int main(int argc, char * argv[])
{
  int i, j, k, m;
  int mm, mb;
  double rd, td, pd;
  double drc, dtc, dpc;

  printf("spherical_deriv.c checking numerical derivatives \n");
  for(i=0; i<nr; i++) printf("rval[%d]=%f \n", i, rval[i]);
  printf("\n");

  tval[0] = -Pi;
  td = 2.0*Pi/(double)(nt);
  for(j=1; j<nt; j++)
  {
    tval[j] = tval[j-1]+td;
  }
  for(j=0; j<nt; j++) printf("tval[%d]=%f \n", j, tval[j]);
  printf("\n");

  pval[0] = -Pi/2;
  pd = Pi/(double)(np);
  for(k=1; k<np; k++)
  {
    pval[k] = pval[k-1]+pd;
  }
  for(j=0; j<np; j++) printf("pval[%d]=%f \n", j, pval[j]);
  printf("\n");

  for(i=0; i<nr; i++)
  {
    for(j=0; j<nt; j++)
    {
      for(k=0; k<np; k++)
      {
        rtp[i][j][k] = u(rval[i], tval[j], pval[k]);
      }
    }
  }

  for(i=0; i<nr; i+=nr-1)
  {
    for(j=0; j<nt; j+=nt-1)
    {
      for(k=0; k<np; k+=np-1)
      {
        printf("rtp[%d][%d][%d]=%f \n", i, j, k, rtp[i][j][k]);
      }
    }
    printf("\n");
  }

  for(i=0; i<nr; i+=nr-1)
  {
    for(j=0; j<nt; j+=nt-1)
    {
      for(k=0; k<np; k+=np-1)
      {
        printf("at %d,%d,%d r=%f, t=%f, p=%f \n", 
               i, j, k, rval[i], tval[j], pval[k]);
        printf("  analytic    dudr=%f,  dudt=%f,  dudp=%f \n",
               dudr(rval[i], tval[j], pval[k]),
               dudt(rval[i], tval[j], pval[k]),
               dudp(rval[i], tval[j], pval[k]));
        nuderiv(1, nr, i, rval, rcoef);
        drc = 0.0;
        for(m=0; m<nr; m++) drc = drc + rcoef[m]*rtp[m][j][k];
        nuderiv(1, nt, j, tval, tcoef);
        dtc = 0.0;
        for(m=0; m<nt; m++) dtc = dtc + tcoef[m]*rtp[i][m][k];
        nuderiv(1, np, k, pval, pcoef);
        dpc = 0.0;
        for(m=0; m<np; m++) dpc = dpc + pcoef[m]*rtp[i][j][m];
        printf("  computed    dudr=%f,  dudt=%f,  dudp=%f \n\n",
               drc, dtc, dpc);
      }
    }
    printf("\n");
  }
  for(m=0; m<nr; m++) printf("rcoef[%d]=%f \n", m, rcoef[m]);
  for(m=0; m<nt; m++) printf("tcoef[%d]=%f \n", m, tcoef[m]);
  for(m=0; m<np; m++) printf("pcoef[%d]=%f \n", m, pcoef[m]);
  printf("/n");
  i = nr/2;
  j = nt/4;
  k = np/4;
  mm = 5;
  mb = 2;
        printf("mm=%d, mb=%d \n", mm, mb);
        printf("at %d,%d,%d r=%f, t=%f, p=%f \n",
               i, j, k,  rval[i], tval[j], pval[k]);
        printf("  analytic    dudr=%f,  dudt=%f,  dudp=%f \n",
               dudr(rval[i], tval[j], pval[k]),
               dudt(rval[i], tval[j], pval[k]),
               dudp(rval[i], tval[j], pval[k]));
        nuderiv(1, mm, mb, &rval[i-mb], rcoef);
        drc = 0.0;
        for(m=0; m<mm; m++) drc = drc + rcoef[m]*rtp[i-mb+m][j][k];
        nuderiv(1, mm, mb, &tval[j-mb], tcoef);
        dtc = 0.0;
        for(m=0; m<mm; m++) dtc = dtc + tcoef[m]*rtp[i][j-mb+m][k];
        nuderiv(1, mm, mb, &pval[k-mb], pcoef);
        dpc = 0.0;
        for(m=0; m<mm; m++) dpc = dpc + pcoef[m]*rtp[i][j][k-mb+m];
        printf("  computed    dudr=%f,  dudt=%f,  dudp=%f \n\n",
               drc, dtc, dpc);
        for(m=0; m<mm; m++) printf("rcoef[%d]=%f \n", m, rcoef[m]);
        for(m=0; m<mm; m++) printf("tcoef[%d]=%f \n", m, tcoef[m]);
        for(m=0; m<mm; m++) printf("pcoef[%d]=%f \n", m, pcoef[m]);
	printf("/n");
  mm = 7;
  mb = 3;
        printf("mm=%d, mb=%d \n", mm, mb);
        printf("at %d,%d,%d r=%f, t=%f, p=%f \n",
               i, j, k,  rval[i], tval[j], pval[k]);
        printf("  analytic    dudr=%f,  dudt=%f,  dudp=%f \n",
               dudr(rval[i], tval[j], pval[k]),
               dudt(rval[i], tval[j], pval[k]),
               dudp(rval[i], tval[j], pval[k]));
        nuderiv(1, mm, mb, &rval[i-mb], rcoef);
        drc = 0.0;
        for(m=0; m<mm; m++) drc = drc + rcoef[m]*rtp[i-mb+m][j][k];
        nuderiv(1, mm, mb, &tval[j-mb], tcoef);
        dtc = 0.0;
        for(m=0; m<mm; m++) dtc = dtc + tcoef[m]*rtp[i][j-mb+m][k];
        nuderiv(1, mm, mb, &pval[k-mb], pcoef);
        dpc = 0.0;
        for(m=0; m<mm; m++) dpc = dpc + pcoef[m]*rtp[i][j][k-mb+m];
        printf("  computed    dudr=%f,  dudt=%f,  dudp=%f \n\n",
               drc, dtc, dpc);
        for(m=0; m<mm; m++) printf("rcoef[%d]=%f \n", m, rcoef[m]);
        for(m=0; m<mm; m++) printf("tcoef[%d]=%f \n", m, tcoef[m]);
        for(m=0; m<mm; m++) printf("pcoef[%d]=%f \n", m, pcoef[m]);

  printf("\nend spherical_deriv check \n");
  return 0;
} /* end main of spherical_deriv.c */