clear all; format compact

%%%% parameters and flags
n=20;
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

%1D differentiation matrices
scale=-2/Lx;
[x,DM] = chebdif(Nx,2); 
dx=DM(:,:,1)*scale;    
dxx=DM(:,:,2)*scale^2;    
x=(x-1)/scale; 

scale=-2/Ly;
[y,DM] = chebdif(Ny,2); 
dy=DM(:,:,1)*scale;    
dyy=DM(:,:,2)*scale^2;    
y=(y-1)/scale; 

scale=-2/Lz;
[z,DM] = chebdif(Nz,2); 
dz=DM(:,:,1)*scale;    
dzz=DM(:,:,2)*scale^2;    
z=(z-1)/scale; 

% 2D differentiation matrices
Dx=kron(dx,speye(Ny));
Dxx=kron(dxx,speye(Ny));
Dy=kron(speye(Nx),dy);
Dyy=kron(speye(Nx),dyy);

% 3D differentiation matrices,here we just add a direction z,and we are now in 3d
DDx=kron(speye(Nz),Dx);
DDxx=kron(speye(Nz),Dxx);
DDy=kron(speye(Nz),Dy);
DDyy=kron(speye(Nz),Dyy);
DDz=kron(dz,speye(Nx*Ny));
DDzz=kron(dzz,speye(Nx*Ny));

[X,Y,Z]=meshgrid(x,y,z);

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

% Numerical derivatives
fX=reshape(DDx*f(:),Ny,Nx,Nz);
fXX=reshape(DDxx*f(:),Ny,Nx,Nz);
fY=reshape(DDy*f(:),Ny,Nx,Nz);
fYY=reshape(DDyy*f(:),Ny,Nx,Nz);
fZ=reshape(DDz*f(:),Ny,Nx,Nz);
fZZ=reshape(DDzz*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
%This code is to print a picture in right size
print('-djpeg','-r80','3d.jpg');