clear all; format compact

% parameters and flags
n=20;   % to change easily the size of the problem
Nx=n+1; % gridpoints in x 
Ny=n+2; % gridpoints in y  
Nz=n+3; % gridpoints in z

Lx=2*pi % domain size in x
Ly=2*pi % domain size in y
Lz=2*pi; % % domain size in z

N=Nx*Ny*Nz % number of degrees of freedom

% 1D and 2D differentiation matrices
[d.x,d.xx,d.wx,x]=dif1D('cheb',0,Lx,Nx,3);
[d.y,d.yy,d.wy,y]=dif1D('cheb',0,Ly,Ny,3);
[d.z,d.zz,d.wz,z]=dif1D('cheb',0,Lz,Nz,3);
[D,l,X,Y,Z,I,NN]=dif2D(d,x,y);

% 3D differentiation matrices
DD.x=kron(speye(Nz),D.x);
DD.xx=kron(speye(Nz),D.xx);
DD.y=kron(speye(Nz),D.y);
DD.yy=kron(speye(Nz),D.yy);
DD.z=kron(d.z,speye(Nx*Ny));
DD.zz=kron(d.zz,speye(Nx*Ny));

% Build the grid itself with meshgrid
[X,Y,Z]=meshgrid(x,y,z);

% Analytical derivatives
f=cos(X).*sin(Y).*cos(Z);
fx=-sin(X).*sin(Y).*cos(Z);
fxx=-cos(X).*sin(Y).*cos(Z);
fy=cos(X).*cos(Y).*cos(Z);
fyy=-cos(X).*sin(Y).*cos(Z);
fz=-cos(X).*sin(Y).*sin(Z);
fzz=-cos(X).*sin(Y).*cos(Z);

% Numerical derivatives
fX=reshape(DD.x*f(:),Ny,Nx,Nz);
fXX=reshape(DD.xx*f(:),Ny,Nx,Nz);
fY=reshape(DD.y*f(:),Ny,Nx,Nz);
fYY=reshape(DD.yy*f(:),Ny,Nx,Nz);
fZ=reshape(DD.z*f(:),Ny,Nx,Nz);
fZZ=reshape(DD.zz*f(:),Ny,Nx,Nz);

% Validation
err_x=max(max(max(abs(fx-fX))))
err_xx=max(max(max(abs(fxx-fXX))))
err_y=max(max(max(abs(fy-fY))))
err_yy=max(max(max(abs(fyy-fYY))))
err_z=max(max(max(abs(fz-fZ))))
err_zz=max(max(max(abs(fzz-fZZ))))

% Loop to show the function f
for ind=1:Nz
    mesh(X(:,:,ind),Y(:,:,ind),f(:,:,ind));
    axis([0 Lx 0 Ly -1 1]);
    xlabel('x');ylabel('y');zlabel('f');
    title(['z=' num2str(z(ind))]);
    drawnow;pause(0.1)
end