// test_nabla8.c  test all functions in  nabla.h  nabla.c
//               just n=8 dimensions, in this file
// del^2 = nabla = Laplacian  n dimensional space of function U(r,a)
//        n-dimensional sphere, radius r
//        a[0],...,a[n-2] are angles in radians (User can rename in their code)
//        compute Laplacian of U(r,a) at given angles
// double nabla(int n, double r, double a[], double (*U()));
//
// alternate call if partial derivative values of U are known
//        da[0],...,da[n-2]   are first partial derivatives of users U(r,a)
//        dda[0],...,dda[n-2] are second partial derivatives of users U(r,a)
//                            with respect to a[0],...,a[n-2] evaluated
//                            at a[0],...,a[n-1]
// double nablapd(int n, double r, double a[],
//                double dr, double ddr,  double da[], double dda[]);
//
// utility function for n-dimensional cartesian to spherical coordinates
// n>3 x[0..n-1]  r, a[0..n-2]
// void toSpherical(int n, double x[], double *r, double a[]);
//
// utility function for n-dimensional spherical to cartesian coordinates 
// n>3 r, a[0..n-2] x[0..n-1]
// void toCartesian(int n, double r, double a[], double x[]);
//
// utility function for n-dimensional U(r, a) at r, a[]
// to derivatives da[] and dda[]
// void sphereDeriv(int n, double (*U)(), double r, double a[],
//                  double *dr, double *ddr, double da[], double dda[]);
//


#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "nabla.h"

// just a test function and its derivatives
double U8f(double r, double a[])
{
  return r*r + 1.1*a[0]*a[0] + 1.2*a[1]*a[1] + 1.3*a[2]*a[2] +
         1.4*a[3]*a[3]+ 1.5*a[4]*a[4] + 1.6*a[5]*a[5] +
         1.7*a[6]*a[6] +
         2.0*r*a[0] + 2.1*r*a[1] + 2.2*r*a[2] + 2.3*r*a[3] +
         2.4*r*a[4] + 2.5*r*a[5] + 2.6*r*a[6] +
         3.0*a[0]*a[1] + 3.1*a[0]*a[2] + 3.2*a[0]*a[3] +
         3.3*a[0]*a[4] + 3.4*a[0]*a[5] + 3.5*a[0]*a[6] + 
         3.7*a[1]*a[2] + 3.8*a[1]*a[3] + 3.9*a[1]*a[4] + 
         4.0*a[1]*a[5] + 4.1*a[1]*a[6] + 4.3*a[2]*a[3] +
         4.4*a[2]*a[4] + 4.5*a[2]*a[5] + 4.6*a[2]*a[6] + 
         4.8*a[3]*a[4] + 4.9*a[3]*a[5] + 5.0*a[3]*a[6] +
         5.2*a[4]*a[5] + 5.3*a[4]*a[6] + 5.5*a[5]*a[6];
}

double U8r(double r, double a[])
{
  return 2.0*r +  2.0*a[0] + 2.1*a[1] + 2.2*a[2] + 2.3*a[3] +
         2.4*a[4] + 2.5*a[5] + 2.6*a[6];
}

double U8rr(double r, double a[])
{
  return 2.0;
}

double U8a0(double r, double a[])
{
  return 2.2*a[0] + 2.0*r    + 3.0*a[1] + 3.1*a[2] + 3.2*a[3] +
         3.3*a[4] + 3.4*a[5] + 3.5*a[6];
}

double U8a0a0(double r, double a[])
{
  return 2.2;
}

double U8a1(double r, double a[])
{
  return 2.4*a[1] + 2.1*r +    3.0*a[0] + 3.7*a[2] + 3.8*a[3] +
         3.9*a[4] + 4.0*a[5] + 4.1*a[6];
}

double U8a1a1(double r, double a[])
{
  return 2.4;
}

double U8a2(double r, double a[])
{
  return 2.6*a[2] + 2.2*r    + 3.1*a[0] + 3.7*a[1] + 4.3*a[3] +
         4.4*a[4] + 4.5*a[5] + 4.6*a[6];
}

double U8a2a2(double r, double a[])
{
  return 2.6;
}

double U8a3(double r, double a[])
{
  return 2.8*a[3] + 2.3*r    + 3.2*a[0] + 3.8*a[1] + 4.3*a[2] +
         4.8*a[4] + 4.9*a[5] + 5.0*a[6];
}

double U8a3a3(double r, double a[])
{
  return 2.8;
}

double U8a4(double r, double a[])
{
  return 3.0*a[4] + 2.4*r    + 3.3*a[0] + 3.9*a[1] + 4.4*a[2] +
         4.8*a[3] + 5.2*a[5] + 5.3*a[6];
}

double U8a4a4(double r, double a[])
{
  return 3.0;
}

double U8a5(double r, double a[])
{
  return 3.2*a[5] + 2.5*r    + 3.4*a[0] + 4.0*a[1] + 4.5*a[2] +
         4.9*a[3] + 5.2*a[4] + 5.5*a[6];
}

double U8a5a5(double r, double a[])
{
  return 3.2;
}

double U8a6(double r, double a[])
{
  return 3.4*a[6] + 2.6*r    + 3.5*a[0] + 4.1*a[1] + 4.6*a[2] +
         5.0*a[3] + 5.3*a[4] + 5.5*a[5];
}

double U8a6a6(double r, double a[])
{
  return 3.4;
}

int main(int argc, char* argv[])
{
  int n = 8;
  double val;
  double r = 0.9;
  double a[7] = {1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8};
  double x[8];
  double rc; // computed for check
  double ac[7]; // computed for check
  double xc[8], xe[8];
  double dr, ddr, da[7], dda[7];
  double drc, ddrc, dac[7], ddac[7];

  printf("test_nabla.c running \n");
  printf("n=%d, r=%f, a[0]=%f, a[1]=%f, a[2]=%f \n", n, r, a[0], a[1], a[2]);
  printf("a[3]=%f, a[4]=%f, a[5]=%f, a[6]=%f \n", a[3], a[4], a[5], a[6]);
  printf(" \n");

  sphereDeriv(n, U8f, r, a, &dr, &ddr, da, dda);
  printf("dr=%f, da0=%f, da1=%f, da2=%f \n", dr, da[0], da[1], da[2]);
  printf("da3=%f, da4=%f, da5=%f, da6=%f \n", da[3], da[4], da[5], da[6]);
  printf("   %f,     %f,     %f,     %f \n",
         U8r(r,a), U8a0(r,a), U8a1(r,a), U8a2(r,a));
  printf("    %f,     %f,     %f,     %f \n",
         U8a3(r,a), U8a4(r,a), U8a5(r,a), U8a6(r,a));
  printf(" \n");
  printf("ddr=%f, dda0=%f, dda1=%f, dda2=%f \n", ddr, dda[0], dda[1], dda[2]);
  printf("dda3=%f, dda4=%f, dda5=%f, dda6=%f \n", dda[3], dda[4], dda[5], dda[6]);
  printf("    %f,      %f,      %f,      %f \n",
         U8rr(r,a), U8a0a0(r,a), U8a1a1(r,a), U8a2a2(r,a));
  printf("     %f,      %f,      %f,      %f \n",
         U8a3a3(r,a), U8a4a4(r,a), U8a5a5(r,a), U8a6a6(r,a));
  printf(" \n");

  toCartesian(n, r, a, x);
  xe[0] = r*cos(a[0]);
  xe[1] = r*sin(a[0])*cos(a[1]);
  xe[2] = r*sin(a[0])*sin(a[1])*cos(a[2]);
  xe[3] = r*sin(a[0])*sin(a[1])*sin(a[2])*cos(a[3]);
  xe[4] = r*sin(a[0])*sin(a[1])*sin(a[2])*sin(a[3])*cos(a[4]);
  xe[5] = r*sin(a[0])*sin(a[1])*sin(a[2])*sin(a[3])*sin(a[4])*cos(a[5]);
  xe[6] = r*sin(a[0])*sin(a[1])*sin(a[2])*sin(a[3])*sin(a[4])*sin(a[5])*cos(a[6]);
  xe[7] = r*sin(a[0])*sin(a[1])*sin(a[2])*sin(a[3])*sin(a[4])*sin(a[5])*sin(a[6]);
  printf("x0 =%f, x1 =%f, x2 =%f, x3 =%f \n", x[0], x[1], x[2], x[3]);
  printf("x4 =%f, x5 =%f, x6 =%f, x7 =%f \n", x[4], x[5], x[6], x[7]);
  printf("xe0=%f, xe1=%f, xe2=%f, xe3=%f \n", xe[0], xe[1], xe[2], xe[3]);
  printf("xe4=%f, xe5=%f, xe6=%f, xe7=%f \n", xe[4], xe[5], xe[6], xe[7]);
  printf(" \n");

  toSpherical(n, x, &rc, ac);
  toCartesian(n, rc, ac, xc);

  printf("r =%f, a0 =%f, a1 =%f, a2 =%f \n", r, a[0], a[1], a[2]);
  printf("a3 =%f, a4 =%f, a5 =%f, a6=%f \n", a[3], a[4], a[5], a[6]);
  printf("rc=%f, ac0=%f, ac1=%f, ac2=%f \n", rc, ac[0], ac[1], ac[2]);
  printf("ac3=%f, ac4=%f, ac5=%f, ac6=%f \n", ac[3], ac[4], ac[5], ac[6]);
  printf(" \n");
  printf("x0 =%f, x1 =%f, x2 =%f, x3 =%f \n", x[0], x[1], x[2], x[3]);
  printf("x4 =%f, x5 =%f, x6 =%f, x7 =%f \n", x[4], x[5], x[6], x[7]);
  printf("xc0=%f, xc1=%f, xc2=%f, xc3=%f \n", xc[0], xc[1], xc[2], xc[3]);
  printf("xc4=%f, xc5=%f, xc6=%f, xc7=%f \n", xc[4], xc[5], xc[6], xc[7]);
  printf(" \n");

  val = nablapd(n, r, a, dr, ddr, da, dda);
  printf("val = nablapd(n, r, a, dr, ddr, da, daa) = %f \n", val);

  val = nabla(n, r, a, U8f);
  printf("val = nabla(n, r, a, U8f) = %f \n", val);

  printf("test_nabla8.c finished \n");
  return 0.0;
} // end test_nabla8.c