#include "libscl.h"
using namespace scl;
using namespace std;

int main(int argc, char** argp, char** envp)
{
 
  bool unit_ball = false;

  
  REAL x0 = 0.0; 
  REAL r = 0.25;
  REAL inc = 0.02;
  REAL left = x0 - r;
  REAL rite = x0 + r;
 

  if (unit_ball) {
    x0 = 0.0; 
    r = 1.0;
    inc = 0.02;
    left = x0 - r;
    rite = x0 + r;
  }

  string filename = "bump_fuction.dat";
  ofstream fout;
  fout.open(filename.c_str());
  if (!fout) error("Error, bump_function, cannot open " + filename);

  cout << starbox("/x, phi, d1phi, d2phi//") << '\n';

  for (REAL x=left; x<=rite; x+=inc) {
    cout << fmt('f',10,5,x)
         << fmt('f',10,5,phi(x,x0,r))
         << fmt('f',10,5,d1phi(x,x0,r))
         << fmt('f',10,5,d2phi(x,x0,r))
         << '\n';
    fout << fmt('f',10,5,x)
         << fmt('f',10,5,phi(x,x0,r))
         << fmt('f',10,5,d1phi(x,x0,r))
         << fmt('f',10,5,d2phi(x,x0,r))
         << '\n';
  }

  fout.clear(); fout.close();

  cout << starbox("/x, first derivative check//") << '\n';

  for (REAL x=left; x<=rite; x+=inc) {
    REAL h = 1.0e-08;
    cout << fmt('f',10,5,x)
         << fmt('f',10,5,phi(x,x0,r))
         << fmt('f',10,5,d1phi(x,x0,r))
         << fmt('f',10,5,(phi(x+h,x0,r) - phi(x,x0,r))/h)
         << '\n';
  }

  cout << starbox("/x, second derivative check//") << '\n';

  for (REAL x=left; x<=rite; x+=inc) {
    REAL h = 1.0e-08;
    cout << fmt('f',10,5,x)
         << fmt('f',10,5,phi(x,x0,r))
         << fmt('f',10,5,d2phi(x,x0,r))
         << fmt('f',10,5,(d1phi(x+h,x0,r) - d1phi(x,x0,r))/h)
         << '\n';
  }

  cout << starbox("/multivariate x, phi, d1phi, d2phi//") << '\n';

  INTEGER d = 2;
  const realmat X0(d,1,0.0);
  const realmat R(d,1,0.25);
  const realmat INC(d,1,0.02);
  const realmat Left = X0 - R;
  const realmat Rite = X0 + R;

  realmat X = X0;
  INTEGER i = 1;
  INTEGER j = 2;

  for (REAL x=Left[i]; x<=Rite[i]; x+=INC[i]) {
    X[i] = x;
    X[j] = x;
    cout << fmt('f',10,5,X[i])
         << fmt('f',10,5,phi(X,X0,R))
         << fmt('f',10,5,d1phi(X,X0,R)[i])
         << fmt('f',10,5,d2phi(X,X0,R)(i,i))
         << fmt('f',10,5,d1phi(X,X0,R)[j])
         << fmt('f',10,5,d2phi(X,X0,R)(i,j))
         << fmt('f',10,5,d2phi(X,X0,R)(j,j))
         << '\n';
  }

  cout << starbox("/multivariate x, first derivative check//") << '\n';

  for (REAL x=Left[i]; x<=Rite[i]; x+=INC[i]) {
    X[i] = x;
    X[j] = x;
    REAL h = 1.0e-08;
    realmat H(d,1,0.0);
    H[i] = h;
    cout << fmt('f',10,5,X[i])
         << fmt('f',10,5,phi(X,X0,R))
         << fmt('f',10,5,d1phi(X,X0,R)[i])
         << fmt('f',10,5,(phi(X+H,X0,R) - phi(X,X0,R))/h)
         << '\n';
  }
  cout << '\n';
  for (REAL x=Left[i]; x<=Rite[i]; x+=INC[i]) {
    X[i] = x;
    X[j] = x;
    REAL h = 1.0e-08;
    realmat H(d,1,0.0);
    H[j] = h;
    cout << fmt('f',10,5,X[i])
         << fmt('f',10,5,phi(X,X0,R))
         << fmt('f',10,5,d1phi(X,X0,R)[j])
         << fmt('f',10,5,(phi(X+H,X0,R) - phi(X,X0,R))/h)
         << '\n';
  }


  cout << starbox("/multivariate x, second derivative check//") << '\n';

  for (REAL x=Left[i]; x<=Rite[i]; x+=INC[i]) {
    X[i] = x;
    X[j] = x;
    REAL h = 1.0e-08;
    realmat H(d,1,0.0);
    H[i] = h;
    cout << fmt('f',10,5,X[i])
         << fmt('f',10,5,phi(X,X0,R))
         << fmt('f',10,5,d2phi(X,X0,R)(i,i))
         << fmt('f',10,5,(d1phi(X+H,X0,R)[i] - d1phi(X,X0,R)[i])/h)
         << '\n';
  }
  cout << '\n';
  for (REAL x=Left[i]; x<=Rite[i]; x+=INC[i]) {
    X[i] = x;
    X[j] = x;
    REAL h = 1.0e-08;
    realmat H(d,1,0.0);
    H[j] = h;
    cout << fmt('f',10,5,X[i])
         << fmt('f',10,5,phi(X,X0,R))
         << fmt('f',10,5,d2phi(X,X0,R)(j,j))
         << fmt('f',10,5,(d1phi(X+H,X0,R)[j] - d1phi(X,X0,R)[j])/h)
         << '\n';
  }
  cout << '\n';
  for (REAL x=Left[i]; x<=Rite[i]; x+=INC[i]) {
    X[i] = x;
    X[j] = x;
    REAL h = 1.0e-08;
    realmat H(d,1,0.0);
    H[j] = h;
    cout << fmt('f',10,5,X[i])
         << fmt('f',10,5,phi(X,X0,R))
         << fmt('f',10,5,d2phi(X,X0,R)(i,j))
         << fmt('f',10,5,(d1phi(X+H,X0,R)[i] - d1phi(X,X0,R)[i])/h)
         << '\n';
  }


  cout << starbox("/univariate quad rule//");

  realmat a;
  realmat w;
  INTEGER n = 51;

  cout << '\n';
  cout << "x0 = " << x0 << '\n';
  cout << "r = " << r << '\n';
  cout << "n = " << n << '\n';
  cout << '\n';

  bfquad(x0,r,n,a,w);

  for (INTEGER i=1; i<=n; ++i) {
    cout << fmt('d',5,i) << fmt('f',10,5,a[i]) << fmt('f',10,5,w[i]) << '\n';
  }

  cout << '\n';
  REAL sum0 = 0.0;
  for (INTEGER i=1; i<=n; ++i) sum0 += w[i]*phi(a[i],x0,r);

  REAL sum1 = 0.0;
  for (REAL x=left; x<=rite; x+=inc) sum1 += phi(x,x0,r)*inc;

  cout << fmt('f',15,8,sum0) << fmt('f',15,8,sum1) << '\n';

  cout << '\n';
  sum0 = 0.0;
  for (INTEGER i=1; i<=n; ++i) sum0 += a[i]*w[i]*d2phi(a[i],x0,r);

  sum1 = 0.0;
  for (REAL x=left; x<=rite; x+=inc) sum1 += x*d2phi(x,x0,r)*inc;

  cout << fmt('f',15,8,sum0) << fmt('f',15,8,sum1) << '\n';

  REAL sumw = 0.0;
  for (INTEGER i=1; i<=n; ++i) sumw += w[i];

  cout << '\n';

  cout << "Sum of the weights" << fmt('f',15,8,sumw) << '\n';

  cout << starbox("/multivariate quad rule//");

  d = 2;
  n = 10;

  realmat XX0(d,1,0.0);
  XX0[1] = 2.0;

  realmat RR(d,1,0.5);

  vector<realmat> A;
  bfquad(XX0,RR,n,A,w);

  INTEGER ngrid = w.size();

  cout << '\n';
  cout << "d = " << d << '\n';
  cout << "n = " << n << '\n';
  cout << "ngrid = " << ngrid << '\n';
  cout << "XX0 = " << XX0 << '\n';
  cout << "RR = " << RR << '\n';
  cout << "A = " << '\n';
  for (INTEGER j=0; j<ngrid; ++j) {
    for (INTEGER i=1; i<=d; ++i) cout << fmt('f',7,3,A[j][i]);
    cout << '\n';
  }
  cout << '\n';
  cout << "w = " << w << '\n';

  sumw = 0.0;
  for (INTEGER i=1; i<=ngrid; ++i) sumw += w[i];

  cout << '\n';

  cout << "Sum of the weights" << fmt('f',15,8,sumw) << '\n';
  cout << '\n';

  {
  //heat equation: (d/dt) u(x,t) = alpha * (d^2/x^2) u(x,t), alpha = 1
  //range        : 0 < x < pi, 0 < t < 3
  //boundary cond: u(0,t) = u(pi,t) = T, u(x,0) = sin(x), T = 50
  //solution     : u(x,t) = T + exp(-alpha*t)*sin(x)
  realmat x0(2,1,1.0);
  realmat radius(2,1,0.75);
  INTEGER nquad = 150;
  vector<realmat> a;
  realmat w;
  bfquad(x0,radius,nquad,a,w);
  INTEGER npoints = w.size();
  REAL lhs = 0.0; REAL rhs = 0.0;
  for (INTEGER j=1; j<=npoints; ++j) {
    realmat x = a[j-1];
    REAL u = 50.0 + exp(-x[2])*sin(x[1]);
    realmat H = d2phi(x,x0,radius);
    rhs += u*H(1,1)*w[j];
    realmat D = d1phi(x,x0,radius);
    lhs -= u*D[2]*w[j];
  }
  cout << "lhs = "<< lhs << "\nrhs = " << rhs << '\n';
  }

  return 0;
}
