// nabla.c  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 *drr, double da[], double dda[]);
//
// method for basic nabla_n: 
// use iterative function, n>3, for computing nabla_n of U(r,a[])
// using d for partial derivative symbol and U for U(r,a[]) and
// adjusting subscripts for a[0],...,a[n-2]
//
// before reduction for numerical computation:
// nabla_n = 1/r^n-1 d(r^n-1 dU/dr)/dr + 1/r^2 L2(n)
//
// after reduction:
// nabla_n = (n-1)/r dU/dr + d^2U/dr^2 + 1/r^2 L2(n, r, a)
//
// in code:
// nabla_n = ((n-1)/r)*dr + ddr + (1.0/(r*r))*L2(n,r,a);
//
// before reduction and adjusting subscripts of a_i code a[]:
// L2(n,r,a) = sum(i=2,n){prod(j=i+1,n){1/sin^2(a_j)} *
//                        1/sin(a_i)^(i-2)*d(sin(a_i)^(i-2)*dU/da_i)/da_i}
//
// after reduction:
// L2(n,r,a) = sum(i=2,n){{prod(j=i+1,n){1/sin^2(a_j)} *
//             (i-2)*cos(a_i)/sin(a_i) *dU/da_i + d^2U/da_i^2}
//
// after adjusting subscripts:
// L2(n,r,a) = sum(i=2,n){{prod(j=i+1,n){1/sin^2(a_j-2)} *
//             ((i-2)*cos(a_i-2)/sin(a_i-2)) * dU/da_i-2 + d^2U/(da_i-2)^2}
//
// in code (n>3):
// tmp = 0.0;
// for(i=2; i<=n; i++)
// {
//   ptmp = 1.0;
//   for(j=i+1; j<=n; j++)
//   {
//     ptmp = ptmp * (1.0/sin(a[j-2])*sin(a[j-2]));
//   }
//   ptmp = ptmp * ((i-2.0)*cos(s[i-2])/sin(a[i-2])) * da[i-2] * dda[i-2];
//   tmp = tmp + ptmp;
// }
// L2(n,r,a) = tmp;

#include <math.h>
#include "nabla.h"
#define Pi 3.1415926535897932384626433832795028841971

double nabla(int n, double r, double a[], double (*U)())
{
  double dr, da[n];
  double ddr, dda[n];
  sphereDeriv(n, U, r, a, &dr, &ddr, da, dda);
  return  nablapd(n, r, a, dr, ddr, da, dda);
} // end nabla

double nablapd(int n, double r, double a[], 
               double dr, double ddr, double da[], double dda[])
{
  double tmp, ptmp;
  int i, j;

  if(n<=3) return 0.0;
  tmp = 0.0;
  for(i=2; i<=n; i++)
  {
    ptmp = 1.0;
    for(j=i+1; j<=n; j++)
    {
      ptmp = ptmp * (1.0/sin(a[j-2])*sin(a[j-2]));
    }
    ptmp = ptmp * ((i-2.0)*cos(a[i-2])/sin(a[i-2])) * da[i-2] * dda[i-2];
    tmp = tmp + ptmp;
  }
  tmp = ((n-1.0)/r)*dr + ddr + (1.0/(r*r))*tmp; // tmp is L(n,r,a)
  return tmp;
} // end nablapd

// example 8 Dimensional sphere, n=8 symmetry with zero based indexing 
// x0 = r cos(a0) 
// x1 = r sin(a0) cos(a1) 
// x2 = r sin(a0) sin(a1) cos(a2) 
// x3 = r sin(a0) sin(a1) sin(a2) cos(a3)
// x4 = r sin(a0) sin(a1) sin(a2) sin(a3) cos(a4) 
// x5 = r sin(a0) sin(a1) sin(a2) sin(a3) sin(a4) cos(a5)  
// x6 = r sin(a0) sin(a1) sin(a2) sin(a3) sin(a4) sin(a5) cos(a6) 
// x7 = r sin(a0) sin(a1) sin(a2) sin(a3) sin(a4) sin(a5) sin(a6) 
//
// r = sqrt(x0^2 + x1^2 + x2^2 + x3^2 + x4^2 + x5^2 + x6^2 + x7^2)
// a0 = acos(x0/r)
// a1 = acos(x1/sqrt(x1^2 + x2^2 + x3^2 + x4^2 + x5^2 + x6^2 + x7^2))
// a2 = acos(x2/sqrt(x2^2 + x3^2 + x4^2 + x5^2 + x6^2 + x7^2))
// a3 = acos(x3/sqrt(x3^2 + x4^2 + x5^2 + x6^2 + x7^2))
// a4 = acos(x4/sqrt(x4^2 + x5^2 + x6^2 + x7^2))
// a5 = acos(x5/sqrt(x5^2 + x6^2 + x7^2))
// a6 = acos(x6/sqrt(x6^2 + x7^2))      if x7>=0
// a6 = 2Pi-acos(x6/sqrt(x6^2 + x6^2))  if x7<0

void toCartesian(int n, double r, double a[], double x[])
{
  int i, j;
  for(i=0; i<n; i++) x[i] = r;
  for(j=0; j<n; j++)
  {
    if(j<n-1) x[j] = x[j] * cos(a[j]);
    for(i=0; i<j; i++) x[j] = x[j] * sin(a[i]);
  }
  
} // end toCartesian

void toSpherical(int n, double x[], double *r, double a[])
{
  int i;
  double rr = 0.0;

  for(i=0; i<n; i++) rr = rr + x[i]*x[i];
  *r = sqrt(rr);
  for(i=0; i<n-1; i++)
  {
    a[i] = acos(x[i]/sqrt(rr));
    rr = rr - x[i]*x[i];
  }
  if(x[n-1]<0.0) a[n-2] = 2.0*Pi - a[n-2]; 
} // end toSpherical

void sphereDeriv(int n, double (*U)(), double r, double a[],
                 double *dr, double *ddr, double da[], double dda[])
{
  int i, j;
  double h = 0.01;
  double ta[n];

  *dr = (1.0/(60.0*h))*(-U(r-3.0*h, a) + 9.0*U(r-2.0*h, a) 
                        -45.0*U(r-h, a) + 45.0*U(r+h, a) 
                        -9.0*U(r+2.0*h, a) + U(r+3.0*h, a) );
  *ddr = (1.0/(180.0*h*h))*(2.0*U(r-3.0*h, a) - 27.0*U(r-2.0*h, a)
			    +270.0*U(r-h, a) -490.0*U(r, a)
                            +270.0*U(r+h, a) - 27.0*U(r+2.0*h, a)
                            +2.0*U(r+3.0*h, a)    );
  for(i=0; i<n-1; i++)
  {
    for(j=0; j<n-1; j++) ta[j] = a[j];
    ta[i] = ta[i]-3.0*h;
    da[i] = -U(r, ta);
    dda[i] = dda[i] + 2.0*U(r, ta);
    ta[i] = ta[i] + h; // now -2.0*h
    da[i] = da[i] + 9.0*U(r, ta); 
    dda[i] = dda[i] - 27.0*U(r, ta);
    ta[i] = ta[i] + h; // now -h
    da[i] = da[i] - 45.0*U(r, ta);
    dda[i] = dda[i] + 270.0*U(r, ta);
    ta[i] = ta[i] + h; // now original, no da term
    dda[i] = dda[i] - 490.0*U(r, ta);
    ta[i] = ta[i] + h; // now +h
    da[i] = da[i] + 45.0*U(r, ta); 
    dda[i] = dda[i] + 270.0*U(r, ta); 
    ta[i] = ta[i] + h; // now +2.0h
    da[i] = da[i] - 9.0*U(r, ta);
    dda[i] = dda[i] - 27.0*U(r, ta);
    ta[i] = ta[i] + h; // now +3.0h
    da[i] = da[i] + U(r, ta);
    dda[i] = dda[i] + 2.0*U(r, ta);
    da[i] = da[i] * (1.0/(60.0*h));
    dda[i] = dda[i] * (1.0/(180.0*h*h));
  }
} // end sphereDeriv

// just a test function
double U4(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];
}

// just a test call
void test()
{
  int n = 4;
  double val;
  double r = 1.5;
  double a[3] = {2.6, 2.8, 3.1};
  val =  nabla(n, r, a, U4);
} // end test