%
%  Compute the sparse grid quadrature values for a specified level ell and
%  using either Clenshaw-Curtis or trapezoidal quadrture points and weights. 
%

  clear all

  ell = 4;

  for j=1:ell
    if j==1
      n_ell = 1;
      xhat = 0;
      what = 2;
    else
      n_ell = 2^(j-1)+1;
      N = n_ell-1;
      [xhat,what] = clencurt(N);
    end
    x = (1/2)*(1-xhat);
    w = (1/2)*what;

    if j==1
      x1 = x;  y1 = x';
      w1 = w';
      d1 = w1;
      n1 = n_ell;
      wt11 = w1*w1';
      fvals11 = fun_exp(x1,y1,n1,n1);
      F11 =fvals11.*wt11;
      Q(j) = sum(sum(F11));
      wprev = [0; 1; 0];
    elseif j==2
      x2 = x;  y2 = x';
      w2 = w';
      d2 = w2-wprev;
      n2 = n_ell;
      wt21 = d2*d1';    wt12 = d1*d2';
      fvals12 = fun_exp(x1,y2,n1,n2);  fvals21 = fun_exp(x2,y1,n2,n1);
      F12 = fvals12.*wt12;   F21 = fvals21.*wt21;
      Q(j) = Q(j-1) + sum(sum(F12)) + sum(sum(F21));
      wprev = reshape([w; zeros(size(w))],[],1);
      wprev(2*n_ell) = [];
    elseif j==3
      x3 = x;  y3 = x';
      w3 = w';
      d3 = w3-wprev;
      n3 = n_ell;
      wt31 = d3*d1';   wt13 = d1*d3';  wt22 = d2*d2';
      fvals13 = fun_exp(x1,y3,n1,n3);  fvals31 = fun_exp(x3,y1,n3,n1);
      fvals22 = fun_exp(x2,y2,n2,n2);
      F13 = fvals13.*wt13;   F31 = fvals31.*wt31;  F22 = fvals22.*wt22;
      Q(j) = Q(j-1) + sum(sum(F13)) + sum(sum(F31)) + sum(sum(F22));
      wprev = reshape([w; zeros(size(w))],[],1);
      wprev(2*n_ell) = [];
    elseif j==4
      x4 = x;  y4 = x';
      w4 = w';
      d4 = w4-wprev;
      n4 = n_ell;
      wt41 = d4*d1';   wt14 = d1*d4';  wt23 = d2*d3';  wt32 = d3*d2';
      fvals14 = fun_exp(x1,y4,n1,n4);  fvals41 = fun_exp(x4,y1,n4,n1);
      fvals23 = fun_exp(x2,y3,n2,n3);  fvals32 = fun_exp(x3,y2,n3,n2);
      F14 = fvals14.*wt14;   F41 = fvals41.*wt41;  
      F23 = fvals23.*wt23;   F32 = fvals32.*wt32;
      Q(j) = Q(j-1) + sum(sum(F14)) + sum(sum(F41)) ...
                    + sum(sum(F23)) + sum(sum(F32));
      wprev = reshape([w; zeros(size(w))],[],1);
      wprev(2*n_ell) = [];
    end
  end

  %Int = (exp(1)-1)^2
  Int = exp(1) - 2
  Q   
  abs(Int*ones(size(Q)) - Q)

  Q1 = fun_exp(0.5,0.5,1,1)
  Q2 = Q1 + (1/6)*fun_exp(0.5,0,1,1) - (1/3)*fun_exp(0.5,0.5,1,1) + (1/6)*fun_exp(0.5,1,1,1) ...
          + (1/6)*fun_exp(0,0.5,1,1) - (1/3)*fun_exp(0.5,0.5,1,1) + (1/6)*fun_exp(1,0.5,1,1)