// test_nabla.c  test all functions in  nabla.h  nabla.c
//               just n=4 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 U4f(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] +
         2.0*r*a[0] + 2.1*r*a[1] + 2.2*r*a[2] +
         3.0*a[0]*a[1] + 3.1*a[0]*a[2] + 3.2*a[1]*a[2];
}

double U4r(double r, double a[])
{
  return 2.0*r + 2.0*a[0] + 2.1*a[1] + 2.2*a[2];
}

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

double U4a0(double r, double a[])
{
  return 2.2*a[0] + 2.0*r + 3.0*a[1] + 3.1*a[2];
}

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

double U4a1(double r, double a[])
{
  return 2.4*a[1] + 2.1*r + 3.0*a[0] + 3.2*a[2];
}

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

double U4a2(double r, double a[])
{
  return 2.6*a[2] + 2.2*r + 3.1*a[0] + 3.2*a[1];
}

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

int main(int argc, char* argv[])
{
  int n = 4;
  double val;
  double r = 1.5;
  double a[3] = {1.9, 2.1, 2.3};
  double x[4];
  double rc; // computed for check
  double ac[3]; // computed for check
  double xc[4], xe[4];
  double dr, ddr, da[3], dda[3];
  double drc, ddrc, dac[3], ddac[3];

  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(" \n");

  sphereDeriv(n, U4f, r, a, &dr, &ddr, da, dda);
  printf("dr=%f, da0=%f, da1=%f, da2=%f \n", dr, da[0], da[1], da[2]);
  printf("   %f,     %f,     %f,     %f \n",
         U4r(r,a), U4a0(r,a), U4a1(r,a), U4a2(r,a));
  printf(" \n");
  printf("ddr=%f, dda0=%f, dda1=%f, dda2=%f \n", ddr, dda[0], dda[1], dda[2]);
  printf("    %f,      %f,      %f,      %f \n", 2.0, 2.2, 2.4, 2.6);
  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]);
  printf("x0 =%f, x1 =%f, x2 =%f, x3 =%f \n", x[0], x[1], x[2], x[3]);
  printf("xe0=%f, xe1=%f, xe2=%f, xe3=%f \n", xe[0], xe[1], xe[2], xe[3]);
  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("rc=%f, ac0=%f, ac1=%f, ac2=%f \n", rc, ac[0], ac[1], ac[2]);
  printf(" \n");
  printf("x0 =%f, x1 =%f, x2 =%f, x3 =%f \n", x[0], x[1], x[2], x[3]);
  printf("xc0=%f, xc1=%f, xc2=%f, xc3=%f \n", xc[0], xc[1], xc[2], xc[3]);
  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, U4f);
  printf("val = nabla(n, r, a, U4f) = %f \n", val);

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