#ifndef BASISFUNCTION
#define BASISFUNCTION

#include <math.h>
#include <iostream.h>

//#define E2 7.389056098930650227230427460575
#define EE16 2.1653645317858030703

class BasisFunction{
  
public:
  float centerX, centerY, centerZ;
  float cXX, cYY, cZZ, cXY, cYZ, cZX, cX, cY, cZ, c0;
  float support;
  int level;
  
  double value(float x, float y, float z){
    float vx = x - centerX;
    float vy = y - centerY;
    float vz = z - centerZ;
    
    double d = vx*vx+vy*vy+vz*vz;
    if(d > support*support)
      return 0;
    double w = weight(sqrt(d), support);
    
    return w*poly(vx, vy, vz);
  }
  
  void gradient(double g[3], float x, float y, float z){
    float vx = x - centerX;
    float vy = y - centerY;
    float vz = z - centerZ;
    
    double d = sqrt(vx*vx+vy*vy+vz*vz);
    if(d > support || d == 0){
      g[0] = g[1] = g[2] = 0;
      return;
    }
    double w = weight(d, support);
    double gw[3];
    weightG(gw, d, support, vx, vy, vz);
    
    double f = poly(vx, vy, vz);
    double gf[3];
    polyG(gf, vx, vy, vz);
    
    g[0] = gw[0]*f + w*gf[0];
    g[1] = gw[1]*f + w*gf[1];
    g[2] = gw[2]*f + w*gf[2];
  }
  
  void gradient2(double gg[3][3], float x, float y, float z){
    float vx = x - centerX;
    float vy = y - centerY;
    float vz = z - centerZ;
    
    double d = sqrt(vx*vx+vy*vy+vz*vz);
    if(d > support || d == 0){
      gg[0][0] = gg[0][1] = gg[0][2] =
      gg[1][0] = gg[1][1] = gg[1][2] =
      gg[2][0] = gg[2][1] = gg[2][2] = 0;
      return;
    }
    
    double w = weight(d, support);
    double gw[3];
    weightG(gw, d, support, vx, vy, vz);
    double ggw[3][3];
    weightGG(ggw, d, support, vx, vy, vz);
    
    double f = poly(vx, vy, vz);
    double gf[3];
    polyG(gf, vx, vy, vz);
    double ggf[3][3];
    polyGG(ggf, vx, vy, vz);
    
    gg[0][0] = ggf[0][0]*w + 2.0*gf[0]*gw[0] + f*ggw[0][0];
    gg[0][1] = ggf[0][1]*w + gf[0]*gw[1] + gf[1]*gw[0] + f*ggw[0][1];
    gg[0][2] = ggf[0][2]*w + gf[0]*gw[2] + gf[2]*gw[0] + f*ggw[0][2];
    
    gg[1][0] = gg[0][1];
    gg[1][1] = ggf[1][1]*w + 2.0*gf[1]*gw[1] + f*ggw[1][1];
    gg[2][1] = ggf[1][2]*w + gf[1]*gw[2] + gf[2]*gw[1] + f*ggw[1][2];
    
    gg[2][0] = gg[0][2];
    gg[2][1] = gg[1][2];
    gg[2][2] = ggf[2][2]*w + 2.0*gf[2]*gw[2] + f*ggw[2][2];
  }
  
  double valueAndGradient(double g[3], float x, float y, float z){
    float vx = x - centerX;
    float vy = y - centerY;
    float vz = z - centerZ;
    
    double d = sqrt(vx*vx+vy*vy+vz*vz);
    if(d > support){
      g[0] = g[1] = g[2] = 0;
      return 0;
    }
    else if(d == 0){
      g[0] = g[1] = g[2] = 0;
      return weight(0, support)*poly(vx, vy, vz);
    }
    double w = weight(d, support);
    double gw[3];
    weightG(gw, d, support, vx, vy, vz);
    
    double f = poly(vx, vy, vz);
    double gf[3];
    polyG(gf, vx, vy, vz);
    
    g[0] = gw[0]*f + w*gf[0];
    g[1] = gw[1]*f + w*gf[1];
    g[2] = gw[2]*f + w*gf[2];
    
    return w*f;
  }
  
  double valueAndGradient12(double g[3], double gg[3][3],
                            float x, float y, float z){
    float vx = x - centerX;
    float vy = y - centerY;
    float vz = z - centerZ;
    
    double d = sqrt(vx*vx+vy*vy+vz*vz);
    if(d > support){
      g[0] = g[1] = g[2] = 
      gg[0][0] = gg[0][1] = gg[0][2] =
      gg[1][0] = gg[1][1] = gg[1][2] =
      gg[2][0] = gg[2][1] = gg[2][2] = 0;
      return 0;
    }
    else if(d == 0){
      g[0] = g[1] = g[2] = 
      gg[0][0] = gg[0][1] = gg[0][2] =
      gg[1][0] = gg[1][1] = gg[1][2] =
      gg[2][0] = gg[2][1] = gg[2][2] = 0;
      return weight(0, support)*poly(vx, vy, vz);
    }
    
    double w = weight(d, support);
    double gw[3], ggw[3][3];
    weightG12(gw, ggw, d, support, vx, vy, vz);
    
    double f = poly(vx, vy, vz);
    double gf[3];
    polyG(gf, vx, vy, vz);
    double ggf[3][3];
    polyGG(ggf, vx, vy, vz);
    
    g[0] = gw[0]*f + w*gf[0];
    g[1] = gw[1]*f + w*gf[1];
    g[2] = gw[2]*f + w*gf[2];
    
    gg[0][0] = ggf[0][0]*w + 2.0*gf[0]*gw[0] + f*ggw[0][0];
    gg[0][1] = ggf[0][1]*w + gf[0]*gw[1] + gf[1]*gw[0] + f*ggw[0][1];
    gg[0][2] = ggf[0][2]*w + gf[0]*gw[2] + gf[2]*gw[0] + f*ggw[0][2];
    
    gg[1][1] = ggf[1][1]*w + 2.0*gf[1]*gw[1] + f*ggw[1][1];
    gg[1][2] = ggf[1][2]*w + gf[1]*gw[2] + gf[2]*gw[1] + f*ggw[1][2];
    
    gg[2][2] = ggf[2][2]*w + 2.0*gf[2]*gw[2] + f*ggw[2][2];
    
    return w*f;
  }
  
  double valueAndGradient123(double g[3], double gg[3][3], double ggg[3][3][3],
                             float x, float y, float z){
    float vx = x - centerX;
    float vy = y - centerY;
    float vz = z - centerZ;
    
    double d = sqrt(vx*vx+vy*vy+vz*vz);
    if(d > support){
      g[0] = g[1] = g[2] = 
      gg[0][0] = gg[0][1] = gg[0][2] =
      gg[1][1] = gg[1][2] = gg[2][2] =
      ggg[0][0][0] = ggg[0][0][1] = ggg[0][0][2] =
      ggg[0][1][1] = ggg[0][1][2] = ggg[0][2][2] =
      ggg[1][1][1] = ggg[1][1][2] = ggg[1][2][2] =
      ggg[2][2][2] = 0;
      
      return 0;
    }
    else if(d == 0){
      g[0] = g[1] = g[2] = 
      gg[0][0] = gg[0][1] = gg[0][2] =
      gg[1][1] = gg[1][2] = gg[2][2] =
      ggg[0][0][0] = ggg[0][0][1] = ggg[0][0][2] =
      ggg[0][1][1] = ggg[0][1][2] = ggg[0][2][2] =
      ggg[1][1][1] = ggg[1][1][2] = ggg[1][2][2] =
      ggg[2][2][2] = 0;
      
      return weight(0, support)*poly(vx, vy, vz);
    }
    
    double w = weight(d, support);
    double gw[3], ggw[3][3], gggw[3][3][3];
    weightG123(gw, ggw, gggw, d, support, vx, vy, vz);
    
    double f = poly(vx, vy, vz);
    double gf[3];
    polyG(gf, vx, vy, vz);
    double ggf[3][3];
    polyGG(ggf, vx, vy, vz);
    
    g[0] = gw[0]*f + w*gf[0];
    g[1] = gw[1]*f + w*gf[1];
    g[2] = gw[2]*f + w*gf[2];
    
    gg[0][0] = ggf[0][0]*w + 2.0*gf[0]*gw[0] + f*ggw[0][0];
    gg[0][1] = ggf[0][1]*w + gf[0]*gw[1] + gf[1]*gw[0] + f*ggw[0][1];
    gg[0][2] = ggf[0][2]*w + gf[0]*gw[2] + gf[2]*gw[0] + f*ggw[0][2];
    gg[1][1] = ggf[1][1]*w + 2.0*gf[1]*gw[1] + f*ggw[1][1];
    gg[1][2] = ggf[1][2]*w + gf[1]*gw[2] + gf[2]*gw[1] + f*ggw[1][2];
    gg[2][2] = ggf[2][2]*w + 2.0*gf[2]*gw[2] + f*ggw[2][2];
    
    ggg[0][0][0] = 3.0*ggf[0][0]*gw[0] + 3.0*gf[0]*ggw[0][0] + f*gggw[0][0][0];
    ggg[1][1][1] = 3.0*ggf[1][1]*gw[1] + 3.0*gf[1]*ggw[1][1] + f*gggw[1][1][1];
    ggg[2][2][2] = 3.0*ggf[2][2]*gw[2] + 3.0*gf[2]*ggw[2][2] + f*gggw[2][2][2];
    ggg[0][0][1] = 2.0*ggf[0][1]*gw[0] + ggf[0][0]*gw[1] + 2.0*gf[0]*ggw[0][1] + gf[1]*ggw[0][0] + f*gggw[0][0][1];
    ggg[0][0][2] = 2.0*ggf[0][2]*gw[0] + ggf[0][0]*gw[2] + 2.0*gf[0]*ggw[0][2] + gf[2]*ggw[0][0] + f*gggw[0][0][2];
    ggg[0][1][1] = 2.0*ggf[0][1]*gw[1] + ggf[1][1]*gw[0] + 2.0*gf[1]*ggw[0][1] + gf[0]*ggw[1][1] + f*gggw[0][1][1];
    ggg[0][2][2] = 2.0*ggf[0][2]*gw[2] + ggf[2][2]*gw[0] + 2.0*gf[2]*ggw[0][2] + gf[0]*ggw[2][2] + f*gggw[0][2][2];
    ggg[1][2][2] = 2.0*ggf[1][2]*gw[2] + ggf[2][2]*gw[1] + 2.0*gf[2]*ggw[1][2] + gf[1]*ggw[2][2] + f*gggw[1][2][2];
    ggg[1][1][2] = 2.0*ggf[1][2]*gw[1] + ggf[1][1]*gw[2] + 2.0*gf[1]*ggw[1][2] + gf[2]*ggw[1][1] + f*gggw[1][1][2];
    ggg[0][1][2] = ggf[0][1]*gw[2] + ggf[0][2]*gw[1] + ggf[1][2]*gw[0] + gf[0]*ggw[1][2] + gf[1]*ggw[0][2] + gf[2]*ggw[0][1] + f*gggw[0][1][2];
    
    return w*f;
  }
  
  double valueAndGradient1234(double g[3], double gg[3][3],
                              double ggg[3][3][3], double gggg[3][3][3][3],
                              float x, float y, float z){
    float vx = x - centerX;
    float vy = y - centerY;
    float vz = z - centerZ;
    
    double d = sqrt(vx*vx+vy*vy+vz*vz);
    if(d > support){
      g[0] = g[1] = g[2] = 
      gg[0][0] = gg[0][1] = gg[0][2] =
      gg[1][0] = gg[1][1] = gg[1][2] =
      gg[2][0] = gg[2][1] = gg[2][2] =
      ggg[0][0][0] = ggg[0][0][1] = ggg[0][0][2] =
      ggg[0][1][1] = ggg[0][1][2] = ggg[0][2][2] =
      ggg[1][1][1] = ggg[1][1][2] = ggg[1][2][2] =
      ggg[2][2][2] =
      gggg[0][0][0][0] = gggg[0][0][0][1] = gggg[0][0][0][2] =
      gggg[0][0][1][1] = gggg[0][0][1][2] = gggg[0][0][2][2] =
      gggg[0][1][1][1] = gggg[0][1][1][2] = gggg[0][1][2][2] =
      gggg[0][2][2][2] =
      gggg[1][1][1][1] = gggg[1][1][1][2] = gggg[1][1][2][2] =
      gggg[1][2][2][2] = gggg[2][2][2][2] = 0;
      
      return 0;
    }
    else if(d == 0){
      g[0] = g[1] = g[2] = 
      gg[0][0] = gg[0][1] = gg[0][2] =
      gg[1][0] = gg[1][1] = gg[1][2] =
      gg[2][0] = gg[2][1] = gg[2][2] =
      ggg[0][0][0] = ggg[0][0][1] = ggg[0][0][2] =
      ggg[0][1][1] = ggg[0][1][2] = ggg[0][2][2] =
      ggg[1][1][1] = ggg[1][1][2] = ggg[1][2][2] =
      ggg[2][2][2] =
      gggg[0][0][0][0] = gggg[0][0][0][1] = gggg[0][0][0][2] =
      gggg[0][0][1][1] = gggg[0][0][1][2] = gggg[0][0][2][2] =
      gggg[0][1][1][1] = gggg[0][1][1][2] = gggg[0][1][2][2] =
      gggg[0][2][2][2] =
      gggg[1][1][1][1] = gggg[1][1][1][2] = gggg[1][1][2][2] =
      gggg[1][2][2][2] = gggg[2][2][2][2] = 0;
      
      return weight(0, support)*poly(vx, vy, vz);
    }
    
    double w = weight(d, support);
    double gw[3], ggw[3][3], gggw[3][3][3], ggggw[3][3][3][3];
    weightG1234(gw, ggw, gggw, ggggw, d, support, vx, vy, vz);
    
    double f = poly(vx, vy, vz);
    double gf[3];
    polyG(gf, vx, vy, vz);
    double ggf[3][3];
    polyGG(ggf, vx, vy, vz);
    
    g[0] = gw[0]*f + w*gf[0];
    g[1] = gw[1]*f + w*gf[1];
    g[2] = gw[2]*f + w*gf[2];
    
    gg[0][0] = ggf[0][0]*w + 2.0*gf[0]*gw[0] + f*ggw[0][0];
    gg[0][1] = ggf[0][1]*w + gf[0]*gw[1] + gf[1]*gw[0] + f*ggw[0][1];
    gg[0][2] = ggf[0][2]*w + gf[0]*gw[2] + gf[2]*gw[0] + f*ggw[0][2];
    gg[1][1] = ggf[1][1]*w + 2.0*gf[1]*gw[1] + f*ggw[1][1];
    gg[1][2] = ggf[1][2]*w + gf[1]*gw[2] + gf[2]*gw[1] + f*ggw[1][2];
    gg[2][2] = ggf[2][2]*w + 2.0*gf[2]*gw[2] + f*ggw[2][2];
    
    ggg[0][0][0] = 3.0*ggf[0][0]*gw[0] + 3.0*gf[0]*ggw[0][0] + f*gggw[0][0][0];
    ggg[1][1][1] = 3.0*ggf[1][1]*gw[1] + 3.0*gf[1]*ggw[1][1] + f*gggw[1][1][1];
    ggg[2][2][2] = 3.0*ggf[2][2]*gw[2] + 3.0*gf[2]*ggw[2][2] + f*gggw[2][2][2];
    ggg[0][0][1] = 2.0*ggf[0][1]*gw[0] + ggf[0][0]*gw[1] + 2.0*gf[0]*ggw[0][1] + gf[1]*ggw[0][0] + f*gggw[0][0][1];
    ggg[0][0][2] = 2.0*ggf[0][2]*gw[0] + ggf[0][0]*gw[2] + 2.0*gf[0]*ggw[0][2] + gf[2]*ggw[0][0] + f*gggw[0][0][2];
    ggg[0][1][1] = 2.0*ggf[0][1]*gw[1] + ggf[1][1]*gw[0] + 2.0*gf[1]*ggw[0][1] + gf[0]*ggw[1][1] + f*gggw[0][1][1];
    ggg[0][2][2] = 2.0*ggf[0][2]*gw[2] + ggf[2][2]*gw[0] + 2.0*gf[2]*ggw[0][2] + gf[0]*ggw[2][2] + f*gggw[0][2][2];
    ggg[1][2][2] = 2.0*ggf[1][2]*gw[2] + ggf[2][2]*gw[1] + 2.0*gf[2]*ggw[1][2] + gf[1]*ggw[2][2] + f*gggw[1][2][2];
    ggg[1][1][2] = 2.0*ggf[1][2]*gw[1] + ggf[1][1]*gw[2] + 2.0*gf[1]*ggw[1][2] + gf[2]*ggw[1][1] + f*gggw[1][1][2];
    ggg[0][1][2] = ggf[0][1]*gw[2] + ggf[0][2]*gw[1] + ggf[1][2]*gw[0] + gf[0]*ggw[1][2] + gf[1]*ggw[0][2] + gf[2]*ggw[0][1] + f*gggw[0][1][2];
    
    gggg[0][0][0][0] = 6.0*ggf[0][0]*ggw[0][0] + 4.0*gf[0]*gggw[0][0][0] + f*ggggw[0][0][0][0];
    gggg[0][0][0][1] = 3.0*ggf[0][0]*ggw[0][1] + 3.0*ggf[0][1]*ggw[0][0] + 3.0*gf[0]*gggw[0][0][1] + gf[1]*gggw[0][0][0] + f*ggggw[0][0][0][1];
    gggg[0][0][0][2] = 3.0*ggf[0][0]*ggw[0][2] + 3.0*ggf[0][2]*ggw[0][0] + 3.0*gf[0]*gggw[0][0][2] + gf[2]*gggw[0][0][0] + f*ggggw[0][0][0][2];
    gggg[0][0][1][1] = ggf[0][0]*ggw[1][1] + ggf[1][1]*ggw[0][0] + 4.0*ggf[0][1]*ggw[0][1] + 2.0*gf[0]*gggw[0][1][1] + 2.0*gf[1]*gggw[0][0][1] + f*ggggw[0][0][1][1];
    gggg[0][0][1][2] = ggf[0][0]*ggw[1][2] + ggf[1][2]*ggw[0][0] + 2.0*ggf[0][1]*ggw[0][2] + 2.0*ggf[0][2]*ggw[0][1] + gf[1]*gggw[0][0][2] + gf[2]*gggw[0][0][1] + 2.0*gf[0]*gggw[0][1][2] + f*ggggw[0][0][1][2];
    gggg[0][0][2][2] = ggf[0][0]*ggw[2][2] + ggf[2][2]*ggw[0][0] + 4.0*ggf[0][2]*ggw[0][2] + 2.0*gf[0]*gggw[0][2][2] + 2.0*gf[2]*gggw[0][0][2] + f*ggggw[0][0][2][2];
    gggg[0][1][1][1] = 3.0*ggf[1][1]*ggw[0][1] + 3.0*ggf[0][1]*ggw[1][1] + 3.0*gf[1]*gggw[0][1][1] + gf[0]*gggw[1][1][1] + f*ggggw[0][1][1][1];
    gggg[0][1][1][2] = ggf[1][1]*ggw[0][2] + ggf[0][2]*ggw[1][1] + 2.0*ggf[0][1]*ggw[1][2] + 2.0*ggf[1][2]*ggw[0][1] + gf[0]*gggw[1][1][2] + gf[2]*gggw[0][1][1] + 2.0*gf[1]*gggw[0][1][2] + f*ggggw[0][1][1][2];
    gggg[0][1][2][2] = ggf[2][2]*ggw[0][1] + ggf[0][1]*ggw[2][2] + 2.0*ggf[0][2]*ggw[1][2] + 2.0*ggf[1][2]*ggw[0][2] + gf[0]*gggw[1][2][2] + gf[1]*gggw[0][2][2] + 2.0*gf[2]*gggw[0][1][2] + f*ggggw[0][1][2][2];
    gggg[0][2][2][2] = 3.0*ggf[2][2]*ggw[0][2] + 3.0*ggf[0][2]*ggw[2][2] + 3.0*gf[2]*gggw[0][2][2] + gf[0]*gggw[2][2][2] + f*ggggw[0][2][2][2];
    gggg[1][1][1][1] = 6.0*ggf[1][1]*ggw[1][1] + 4.0*gf[1]*gggw[1][1][1] + f*ggggw[1][1][1][1];
    gggg[1][1][1][2] = 3.0*ggf[1][1]*ggw[1][2] + 3.0*ggf[1][2]*ggw[1][1] + 3.0*gf[1]*gggw[1][1][2] + gf[2]*gggw[1][1][1] + f*ggggw[1][1][1][2];
    gggg[1][1][2][2] = ggf[1][1]*ggw[2][2] + ggf[2][2]*ggw[1][1] + 4.0*ggf[1][2]*ggw[1][2] + 2.0*gf[1]*gggw[1][2][2] + 2.0*gf[2]*gggw[1][1][2] + f*ggggw[1][1][2][2];
    gggg[1][2][2][2] = 3.0*ggf[2][2]*ggw[1][2] + 3.0*ggf[1][2]*ggw[2][2] + 3.0*gf[2]*gggw[1][2][2] + gf[1]*gggw[2][2][2] + f*ggggw[1][2][2][2];
    gggg[2][2][2][2] = 6.0*ggf[2][2]*ggw[2][2] + 4.0*gf[2]*gggw[2][2][2] + f*ggggw[2][2][2][2];
    
    return w*f;
  }
  
  inline double poly(float vx, float vy, float vz){
    return (cXX*vx+cXY*vy+cX)*vx +
           (cYY*vy+cYZ*vz+cY)*vy +
           (cZZ*vz+cZX*vx+cZ)*vz + c0;
  }
  
  inline void polyG(double g[3], float vx, float vy, float vz){
    g[0] = 2.0*cXX*vx + cXY*vy + cZX*vz + cX;
    g[1] = 2.0*cYY*vy + cYZ*vz + cXY*vx + cY;
    g[2] = 2.0*cZZ*vz + cZX*vx + cYZ*vy + cZ;
  }
  
  inline void polyGG(double gg[3][3], float vx, float vy, float vz){
    gg[0][0] = 2.0*cXX;
    gg[0][1] = cXY;
    gg[0][2] = cZX;
    
    gg[1][1] = 2.0*cYY;
    gg[1][2] = cYZ;
    
    gg[2][2] = 2.0*cZZ;
  }
  
  
  static inline double weight(double d, float R){
    //spline patched Gaussian
	  /*
    if(R < d)
      return 0;
    else{
      d /= R;
      if(d < 0.5f)
        return exp(-8*d*d);
      else{
        d = 1.0 - d;
        d = d*d;
        return EE16*d*d;
      }
    }*/

	  
    if(R < d){
      return 0.0;
    }
    //wendrand's functions
    else{
      d /= R;
      //double tmp = 1.0 - d;
      //tmp *= tmp;
      //tmp *= tmp;
      //return tmp*(4.0*d+1.0);
      //return (1.0 - d)*(1.0 - d);
      //return pow(1.0-d, 4.0)*(4.0*d+1.0);
      return pow(1.0-d, 6.0)*(35.0*d*d + 18.0*d + 3.0);
    }
  }
  
  static inline void weightG(double g[3],
                             double d, float R, float x, float y, float z){
    if(R < d || d == 0){
      g[0] = g[1] = g[2] = 0;
    }
    else{
		
      double d1 = d/R;
      
      //double dw = -20.0*d1*pow(1.0-d1, 3.0);
       
      double tmp = (d1-1.0)*(d1-1.0)*(d1-1.0)*(d1-1.0)*(d1-1.0);
      double dw = 56.0*d1*tmp*(5.0*d1+1.0);
      
      double Rd = R*d;
      double dx = x/Rd;
      double dy = y/Rd;
      double dz = z/Rd;
      
      g[0] = dw*dx;
      g[1] = dw*dy;
      g[2] = dw*dz;
      
      /*
      double t = -20.0*pow(1.0-d/R, 3.0)/(R*R);
      
      g[0] = t*x;
      g[1] = t*y;
      g[2] = t*z;*/
    }
  }
  
  static inline void weightGG(double gg[3][3],
                              double d, float R, float x, float y, float z){
    if(R < d || d == 0){
      gg[0][0] = gg[0][1] = gg[0][2] =
        gg[1][0] = gg[1][1] = gg[1][2] =
          gg[2][0] = gg[2][1] = gg[2][2] = 0;
    }
    else{
      double d1 = d/R;
      double dw = -20.0*d1*pow(1.0-d1, 3.0);
      double ddw = 20.0*(1.0-d1)*(1.0-d1)*(4.0*d1-1);
      
      double Rd = R*d;
      double dx = x/Rd;
      double dy = y/Rd;
      double dz = z/Rd;
      
      double Rd3 = R*d*d*d;
      double dxx = (y*y+z*z)/Rd3;
      double dxy = -x*y/Rd3;
      double dzx = -z*x/Rd3;
      double dyy = (z*z+x*x)/Rd3;
      double dyz = -y*z/Rd3;
      double dzz = (x*x+y*y)/Rd3;

      gg[0][0] = ddw*dx*dx + dw*dxx;
      gg[0][1] = gg[1][0] = ddw*dx*dy + dw*dxy;
      gg[0][2] = gg[2][0] = ddw*dx*dz + dw*dzx;
      
      gg[1][1] = ddw*dy*dy + dw*dyy;
      gg[1][2] = gg[2][1] = ddw*dy*dz + dw*dyz;
      
      gg[2][2] = ddw*dz*dz + dw*dzz;
    }
  }
  
  static inline void weightG12(double g[3], double gg[3][3],
                               double d, float R, float x, float y, float z){
    if(R < d || d == 0){
      g[0] = g[1] = g[2] =
      
      gg[0][0] = gg[0][1] = gg[0][2] =
      gg[1][1] = gg[1][2] = gg[2][2] = 0;
    }
    else{
      double d1 = d/R;
      /*
      double dw = -20.0*d1*pow(1.0-d1, 3.0);
      double ddw = 20.0*(1.0-d1)*(1.0-d1)*(4.0*d1-1);
      double dddw = 120.0*(2.0*d1*d1-3.0*d1+1.0);
      double ddddw = 120.0*(4.0*d1-3);
        */
      
      double tmp = (d1-1.0)*(d1-1.0)*(d1-1.0)*(d1-1.0);
      double ddw = 56.0*tmp*(35.0*d1*d1-4.0*d1-1.0);
      tmp *= (d1-1.0);
      double dw = 56.0*d1*tmp*(5.0*d1+1.0);
      
      double Rd = R*d;
      double dx = x/Rd;
      double dy = y/Rd;
      double dz = z/Rd;
      
      double Rd3 = Rd*d*d;
      double dxx = (y*y+z*z)/Rd3;
      double dxy = -x*y/Rd3;
      double dxz = -z*x/Rd3;
      double dyy = (z*z+x*x)/Rd3;
      double dyz = -y*z/Rd3;
      double dzz = (x*x+y*y)/Rd3;
      
      g[0] = dw*dx;
      g[1] = dw*dy;
      g[2] = dw*dz;
      
      gg[0][0] = ddw*dx*dx + dw*dxx;
      gg[0][1] = ddw*dx*dy + dw*dxy;
      gg[0][2] = ddw*dx*dz + dw*dxz;
      gg[1][1] = ddw*dy*dy + dw*dyy;
      gg[1][2] = ddw*dy*dz + dw*dyz;
      gg[2][2] = ddw*dz*dz + dw*dzz;
    }
  }
  
  static inline void weightG123(double g[3], double gg[3][3], double ggg[3][3][3],
                                double d, float R, float x, float y, float z){
    if(R < d || d == 0){
      g[0] = g[1] = g[2] = 
      
      gg[0][0] = gg[0][1] = gg[0][2] =
      gg[1][1] = gg[1][2] = gg[2][2] = 
      
      ggg[0][0][0] = ggg[0][0][1] = ggg[0][0][2] =
      ggg[0][1][1] = ggg[0][1][2] = ggg[0][2][2] =
      ggg[1][1][1] = ggg[1][1][2] = ggg[1][2][2] =
      ggg[2][2][2] = 0;
    }
    else{
      
      double d1 = d/R;
	  /*
      double dw = -20.0*d1*pow(1.0-d1, 3.0);
      double ddw = 20.0*(1.0-d1)*(1.0-d1)*(4.0*d1-1);
      double dddw = 120.0*(2.0*d1*d1-3.0*d1+1.0);
        */
	  
      double tmp = (d1-1.0)*(d1-1.0)*(d1-1.0);
      double dddw = 1680.0*d1*tmp*(7.0*d1-3.0);
      tmp *= (d1-1.0);
      double ddw = 56.0*tmp*(35.0*d1*d1-4.0*d1-1.0);
      tmp *= (d1-1.0);
      double dw = 56.0*d1*tmp*(5.0*d1+1.0);
      
      double Rd = R*d;
      double dx = x/Rd;
      double dy = y/Rd;
      double dz = z/Rd;
      
      double xx = x*x;
      double yy = y*y;
      double zz = z*z;
      
      double Rd3 = Rd*d*d;
      double dxx = (yy+zz)/Rd3;
      double dxy = -x*y/Rd3;
      double dxz = -z*x/Rd3;
      double dyy = (zz+xx)/Rd3;
      double dyz = -y*z/Rd3;
      double dzz = (xx+yy)/Rd3;
      
      double Rd5 = Rd3*d*d;
      double dxxx = -3.0*x*(yy+zz)/Rd5;
      double dxxy = -y*(-2.0*xx+yy+zz)/Rd5;
      double dxxz = -z*(-2.0*xx+yy+zz)/Rd5;
      double dxyy = -x*(xx-2.0*yy+zz)/Rd5;
      double dxyz = 3.0*x*y*z/Rd5;
      double dxzz = -x*(xx+yy-2.0*zz)/Rd5;
      double dyyy = -3.0*y*(xx+zz)/Rd5;
      double dyyz = -z*(xx-2.0*yy+zz)/Rd5;
      double dyzz = -y*(xx+yy-2.0*zz)/Rd5;
      double dzzz = -3.0*z*(xx+yy)/Rd5;
      
      g[0] = dw*dx;
      g[1] = dw*dy;
      g[2] = dw*dz;
      
      
      gg[0][0] = ddw*dx*dx + dw*dxx;
      gg[0][1] = ddw*dx*dy + dw*dxy;
      gg[0][2] = ddw*dx*dz + dw*dxz;
      gg[1][1] = ddw*dy*dy + dw*dyy;
      gg[1][2] = ddw*dy*dz + dw*dyz;
      gg[2][2] = ddw*dz*dz + dw*dzz;
      
      
      ggg[0][0][0] = dddw*dx*dx*dx + 3.0*ddw*dx*dxx + dw*dxxx;
      ggg[0][0][1] = dddw*dx*dx*dy + ddw*(2.0*dx*dxy+dy*dxx) + dw*dxxy;
      ggg[0][0][2] = dddw*dx*dx*dz + ddw*(2.0*dx*dxz+dz*dxx) + dw*dxxz;
      ggg[0][1][1] = dddw*dy*dy*dx + ddw*(2.0*dy*dxy+dx*dyy) + dw*dxyy;
      ggg[0][1][2] = dddw*dx*dy*dz + ddw*(dx*dyz+dy*dxz+dz*dxy) + dw*dxyz;
      ggg[0][2][2] = dddw*dz*dz*dx + ddw*(2.0*dz*dxz+dx*dzz) + dw*dxzz;
      ggg[1][1][1] = dddw*dy*dy*dy + 3.0*ddw*dy*dyy + dw*dyyy;
      ggg[1][1][2] = dddw*dy*dy*dz + ddw*(2.0*dy*dyz+dz*dyy) + dw*dyyz;
      ggg[1][2][2] = dddw*dz*dz*dy + ddw*(2.0*dz*dyz+dy*dzz) + dw*dyzz;
      ggg[2][2][2] = dddw*dz*dz*dz + 3.0*ddw*dz*dzz + dw*dzzz;
    }
  }
  
  static inline void weightG1234(double g[3], double gg[3][3],
                                 double ggg[3][3][3], double gggg[3][3][3][3],
                                 double d, float R, float x, float y, float z){
    if(R < d || d == 0){
      g[0] = g[1] = g[2] = 
      
      gg[0][0] = gg[0][1] = gg[0][2] =
      gg[1][1] = gg[1][2] = gg[2][2] = 
      
      ggg[0][0][0] = ggg[0][0][1] = ggg[0][0][2] =
      ggg[0][1][1] = ggg[0][1][2] = ggg[0][2][2] =
      ggg[1][1][1] = ggg[1][1][2] = ggg[1][2][2] =
      ggg[2][2][2] = 
      
      gggg[0][0][0][0] = gggg[0][0][0][1] = gggg[0][0][0][2] =
      gggg[0][0][1][1] = gggg[0][0][1][2] = gggg[0][0][2][2] =
      gggg[0][1][1][1] = gggg[0][1][1][2] = gggg[0][1][2][2] =
      gggg[0][2][2][2] =
      gggg[1][1][1][1] = gggg[1][1][1][2] = gggg[1][1][2][2] =
      gggg[1][2][2][2] = gggg[2][2][2][2] = 0;
    }
    else{
      double d1 = d/R;
      /*
      double dw = -20.0*d1*pow(1.0-d1, 3.0);
      double ddw = 20.0*(1.0-d1)*(1.0-d1)*(4.0*d1-1);
      double dddw = 120.0*(2.0*d1*d1-3.0*d1+1.0);
      double ddddw = 120.0*(4.0*d1-3);*/
       
      double tmp = (d1-1.0)*(d1-1.0);
      double ddddw = 1680.0*tmp*(35.0*d1*d1-26.0*d1+3.0);
      tmp *= (d1-1.0);
      double dddw = 1680.0*d1*tmp*(7.0*d1-3.0);
      tmp *= (d1-1.0);
      double ddw = 56.0*tmp*(35.0*d1*d1-4.0*d1-1.0);
      tmp *= (d1-1.0);
      double dw = 56.0*d1*tmp*(5.0*d1+1.0);
      
      double Rd = R*d;
      double dx = x/Rd;
      double dy = y/Rd;
      double dz = z/Rd;
      
      double xx = x*x;
      double yy = y*y;
      double zz = z*z;
      
      double Rd3 = Rd*d*d;
      double dxx = (yy+zz)/Rd3;
      double dxy = -x*y/Rd3;
      double dxz = -z*x/Rd3;
      double dyy = (zz+xx)/Rd3;
      double dyz = -y*z/Rd3;
      double dzz = (xx+yy)/Rd3;
      
      double Rd5 = Rd3*d*d;
      double dxxx = -3.0*x*(yy+zz)/Rd5;
      double dxxy = -y*(-2.0*xx+yy+zz)/Rd5;
      double dxxz = -z*(-2.0*xx+yy+zz)/Rd5;
      double dxyy = -x*(xx-2.0*yy+zz)/Rd5;
      double dxyz = 3.0*x*y*z/Rd5;
      double dxzz = -x*(xx+yy-2.0*zz)/Rd5;
      double dyyy = -3.0*y*(xx+zz)/Rd5;
      double dyyz = -z*(xx-2.0*yy+zz)/Rd5;
      double dyzz = -y*(xx+yy-2.0*zz)/Rd5;
      double dzzz = -3.0*z*(xx+yy)/Rd5;
      
      double Rd7 = Rd5*d*d;
      double dxxxx = 3.0*(4.0*xx-yy-zz)*(yy+zz)/Rd7;
      double dxxxy = x*y*(-6.0*xx+9.0*(yy+zz))/Rd7;
      double dxxxz = x*z*(-6.0*xx+9.0*(yy+zz))/Rd7;
      double dxxyy = (xx*(2.0*xx-11.0*yy+zz)+2.0*yy*yy+(yy-zz)*zz)/Rd7;
      double dxxyz = 3.0*y*z*(-4.0*xx+yy+zz)/Rd7;
      double dxxzz = (xx*(2.0*xx-11.0*zz+yy)+2.0*zz*zz+(zz-yy)*yy)/Rd7;
      double dxyyy = x*y*(-6.0*yy+9.0*(xx+zz))/Rd7;
      double dxyyz = 3.0*x*z*(xx-4.0*yy+zz)/Rd7;
      double dxyzz = 3.0*x*y*(xx-4.0*zz+yy)/Rd7;
      double dxzzz = x*z*(-6.0*zz+9.0*(xx+yy))/Rd7;
      double dyyyy = 3.0*(4.0*yy-xx-zz)*(xx+zz)/Rd7;
      double dyyyz = y*z*(-6.0*yy+9.0*(xx+zz))/Rd7;
      double dyyzz = (yy*(2.0*yy-11.0*zz+xx)+2.0*yy*yy+(zz-xx)*xx)/Rd7;
      double dyzzz = y*z*(-6.0*zz+9.0*(xx+yy))/Rd7;
      double dzzzz = 3.0*(4.0*zz-xx-yy)*(xx+yy)/Rd7;
      
      g[0] = dw*dx;
      g[1] = dw*dy;
      g[2] = dw*dz;
      
      
      gg[0][0] = ddw*dx*dx + dw*dxx;
      gg[0][1] = ddw*dx*dy + dw*dxy;
      gg[0][2] = ddw*dx*dz + dw*dxz;
      gg[1][1] = ddw*dy*dy + dw*dyy;
      gg[1][2] = ddw*dy*dz + dw*dyz;
      gg[2][2] = ddw*dz*dz + dw*dzz;
      
      
      ggg[0][0][0] = dddw*dx*dx*dx + 3.0*ddw*dx*dxx + dw*dxxx;
      ggg[0][0][1] = dddw*dx*dx*dy + ddw*(2.0*dx*dxy+dy*dxx) + dw*dxxy;
      ggg[0][0][2] = dddw*dx*dx*dz + ddw*(2.0*dx*dxz+dz*dxx) + dw*dxxz;
      ggg[0][1][1] = dddw*dy*dy*dx + ddw*(2.0*dy*dxy+dx*dyy) + dw*dxyy;
      ggg[0][1][2] = dddw*dx*dy*dz + ddw*(dx*dyz+dy*dxz+dz*dxy) + dw*dxyz;
      ggg[0][2][2] = dddw*dz*dz*dx + ddw*(2.0*dz*dxz+dx*dzz) + dw*dxzz;
      ggg[1][1][1] = dddw*dy*dy*dy + 3.0*ddw*dy*dyy + dw*dyyy;
      ggg[1][1][2] = dddw*dy*dy*dz + ddw*(2.0*dy*dyz+dz*dyy) + dw*dyyz;
      ggg[1][2][2] = dddw*dz*dz*dy + ddw*(2.0*dz*dyz+dy*dzz) + dw*dyzz;
      ggg[2][2][2] = dddw*dz*dz*dz + 3.0*ddw*dz*dzz + dw*dzzz;
      
      gggg[0][0][0][0] = dx*dx*(ddddw*dx*dx + 6.0*dddw*dxx) + ddw*(3.0*dxx*dxx + 4.0*dx*dxxx) + dw*dxxxx;
      gggg[0][0][0][1] = ddddw*dy*dx*dx*dx + 3.0*dddw*dx*(dx*dxy + dy*dxx) + ddw*(3.0*dxy*dxx + 3.0*dx*dxxy + dy*dxxx) + dw*dxxxy;
      gggg[0][0][0][2] = ddddw*dz*dx*dx*dx + 3.0*dddw*dx*(dx*dxz + dz*dxx) + ddw*(3.0*dxz*dxx + 3.0*dx*dxxz + dz*dxxx) + dw*dxxxz;
      gggg[0][0][1][1] = ddddw*dx*dx*dy*dy + dddw*(dx*dx*dyy + 4.0*dx*dy*dxy + dy*dy*dxx) + ddw*(2.0*dxy*dxy + 2.0*dx*dxyy + dxx*dyy + 2.0*dy*dxxy) + dw*dxxyy;
      gggg[0][0][1][2] = ddddw*dx*dx*dy*dz + dddw*(dx*(dx*dyz + 2.0*dy*dxz + 2.0*dz*dxy) + dy*dz*dxx) + ddw*(2.0*dxz*dxy + 2.0*dx*dxyz + dxx*dyz + dy*dxxz + dz*dxxy) + dw*dxxyz;
      gggg[0][0][2][2] = ddddw*dx*dx*dz*dz + dddw*(dx*dx*dzz + 4.0*dx*dz*dxz + dz*dz*dxx) + ddw*(2.0*dxz*dxz + 2.0*dx*dxzz + dxx*dzz + 2.0*dz*dxxz) + dw*dxxzz;
      gggg[0][1][1][1] = ddddw*dx*dy*dy*dy + 3.0*dddw*dy*(dy*dxy + dx*dyy) + ddw*(3.0*dxy*dyy + 3.0*dy*dxyy + dx*dyyy) + dw*dxyyy;
      gggg[0][1][1][2] = ddddw*dy*dy*dx*dz + dddw*(dy*(dy*dxz + 2.0*dx*dyz + 2.0*dz*dxy) + dx*dz*dyy) + ddw*(2.0*dxy*dyz + 2.0*dy*dxyz + dyy*dxz + dz*dxyy + dx*dyyz) + dw*dxyyz;
      gggg[0][1][2][2] = ddddw*dz*dz*dx*dy + dddw*(dz*(dz*dxy + 2.0*dy*dxz + 2.0*dx*dyz) + dx*dy*dzz) + ddw*(2.0*dxz*dyz + 2.0*dz*dxyz + dzz*dxy + dy*dxzz + dx*dyzz) + dw*dxyzz;
      gggg[0][2][2][2] = ddddw*dx*dz*dz*dz + 3.0*dddw*dz*(dz*dxz + dx*dzz) + ddw*(3.0*dxz*dzz + 3.0*dz*dxzz + dx*dzzz) + dw*dxzzz;
      gggg[1][1][1][1] = dy*dy*(ddddw*dy*dy + 6.0*dddw*dyy) + ddw*(3.0*dyy*dyy + 4.0*dy*dyyy) + dw*dyyyy;
      gggg[1][1][1][2] = ddddw*dz*dy*dy*dy + 3.0*dddw*dy*(dy*dyz + dz*dyy) + ddw*(3.0*dyz*dyy + 3.0*dy*dyyz + dz*dyyy) + dw*dyyyz;
      gggg[1][1][2][2] = ddddw*dy*dy*dz*dz + dddw*(dy*dy*dzz + 4.0*dy*dz*dyz + dz*dz*dyy) + ddw*(2.0*dyz*dyz + 2.0*dy*dyzz + dyy*dzz + 2.0*dz*dyyz) + dw*dyyzz;
      gggg[1][2][2][2] = ddddw*dy*dz*dz*dz + 3.0*dddw*dz*(dz*dyz + dy*dzz) + ddw*(3.0*dyz*dzz + 3.0*dz*dyzz + dy*dzzz) + dw*dyzzz;
      gggg[2][2][2][2] = dz*dz*(ddddw*dz*dz + 6.0*dddw*dzz) + ddw*(3.0*dzz*dzz + 4.0*dz*dzzz) + dw*dzzzz;
    }
  }
};

#endif