function [X,Z,dx,dz,Xxi,Xeta,Zxi,Zeta,CHI] = grid_Kurt(Nx,Nz,x0,x1,g,z1)

    z0 = g(1);
    
    %Space step
    dx = (x1-x0)/(Nx-1);
    dz = (z1-z0)/(Nz-1);
    
    %First derivative 6th order
    a1 = 0.70127127127127;  
    
    D1x = (1/60*diag(ones(Nx-3,1),3)-9/60*diag(ones(Nx-2,1),2)+45/60*diag(ones(Nx-1,1),1)-45/60*diag(ones(Nx-1,1),-1)+9/60*diag(ones(Nx-2,1),-2)-1/60*diag(ones(Nx-3,1),-3));
    D1x(1:6,1:9) = [-21600/13649, 43200/13649*a1-7624/40947, -172800/13649*a1+ ...
	     715489/81894, 259200/13649*a1-187917/13649, -172800/13649* ...
	     a1+735635/81894, 43200/13649*a1-89387/40947, 0, 0, 0; ...
	     -8640/12013*a1+7624/180195, 0, 86400/12013*a1-57139/12013, ...
	     -172800/12013*a1+745733/72078, 129600/12013*a1-91715/12013, ...
	     -34560/12013*a1+240569/120130, 0, 0, 0; 17280/2711*a1-715489/162660, -43200/2711*a1+57139/5422, 0, 86400/2711*a1-176839/8133, -86400/2711*a1+242111/10844, 25920/2711*a1-182261/27110, 0, 0, 0; -25920/5359*a1+187917/53590, 86400/5359*a1-745733/64308, -86400/5359*a1+176839/16077, 0, 43200/5359*a1-165041/32154, -17280/5359*a1+710473/321540, 72/5359, 0, 0; 34560/7877*a1-147127/47262, -129600/7877*a1+91715/7877, 172800/7877*a1-242111/15754, -86400/7877*a1+165041/23631, 0, 8640/7877*a1, -1296/7877, 144/7877, 0; -43200/43801*a1+89387/131403, 172800/43801*a1-240569/87602, -259200/43801*a1+182261/43801, 172800/43801*a1-710473/262806, -43200/43801*a1, 0, 32400/43801, -6480/43801, 720/43801];
    D1x(Nx-5:Nx,Nx-8:Nx) = rot90( -D1x(1:6,1:9),2);
    D1x=D1x/dx;
    
    D1z = (1/60*diag(ones(Nz-3,1),3)-9/60*diag(ones(Nz-2,1),2)+45/60*diag(ones(Nz-1,1),1)-45/60*diag(ones(Nz-1,1),-1)+9/60*diag(ones(Nz-2,1),-2)-1/60*diag(ones(Nz-3,1),-3));
    D1z(1:6,1:9) = [-21600/13649, 43200/13649*a1-7624/40947, -172800/13649*a1+ ...
	     715489/81894, 259200/13649*a1-187917/13649, -172800/13649* ...
	     a1+735635/81894, 43200/13649*a1-89387/40947, 0, 0, 0; ...
	     -8640/12013*a1+7624/180195, 0, 86400/12013*a1-57139/12013, ...
	     -172800/12013*a1+745733/72078, 129600/12013*a1-91715/12013, ...
	     -34560/12013*a1+240569/120130, 0, 0, 0; 17280/2711*a1-715489/162660, -43200/2711*a1+57139/5422, 0, 86400/2711*a1-176839/8133, -86400/2711*a1+242111/10844, 25920/2711*a1-182261/27110, 0, 0, 0; -25920/5359*a1+187917/53590, 86400/5359*a1-745733/64308, -86400/5359*a1+176839/16077, 0, 43200/5359*a1-165041/32154, -17280/5359*a1+710473/321540, 72/5359, 0, 0; 34560/7877*a1-147127/47262, -129600/7877*a1+91715/7877, 172800/7877*a1-242111/15754, -86400/7877*a1+165041/23631, 0, 8640/7877*a1, -1296/7877, 144/7877, 0; -43200/43801*a1+89387/131403, 172800/43801*a1-240569/87602, -259200/43801*a1+182261/43801, 172800/43801*a1-710473/262806, -43200/43801*a1, 0, 32400/43801, -6480/43801, 720/43801];
    D1z(Nz-5:Nz,Nz-8:Nz) = rot90( -D1z(1:6,1:9),2);
    D1z=D1z/dz;
    
    %Construct derivatives in xi and eta direction
    Ix = speye(Nx,Nx);
    Iz = speye(Nz,Nz);
    D1x = sparse(D1x);
    D1z = sparse(D1z);
    
    Dxi = kron(D1x,Iz);
    Deta = kron(Ix,D1z);
    
    %Define the logical domain
    xi = x0:dx:x1;
    eta = g(1):dz:z1;
    [XI ETA] = meshgrid(xi,eta);
    
    % Compute derivative, curvature
    gprime = D1x*g;
    gbis = D1x*gprime;
    g = g';
    gprime = gprime';
    gbis = gbis';
    
    CHI = -gbis./(1+gprime.^2).^(3/2);
    
%     plot(xi,CHI)
%     pause
%     

    % Compute grid in physical domain.
    X = zeros(size(XI));
    Z = zeros(size(ETA));
    
    i = Nz;
    s = 0;
        X(i,:) = xi+(g-z0).*gprime./(1+sqrt(1+gprime.^2));
        plot(xi,X(i,:));grid
        pause
    
    
    for i = 1:Nz
        s = (eta(i)-z1)/(z0-z1);
        X(i,:) = xi+(g-z0).*(1-s.^2).*gprime./(1-s.^2+(1+s.^2).*sqrt(1+gprime.^2));
        Z(i,:) = z1+(g-z1).*2.*s.*sqrt(1+gprime.^2)./(1-s.^2+(1+s.^2).*sqrt(1+gprime.^2));
    end
    
    %%%%%%%%%%%%%%%%%%DEBUGING%%%%%%%%%%%%%%%%%%%%
    close all
    Y = zeros(size(X));
    surf(X,Z,Y);
    view(0,90);
    axis equal
    pause
    
    %Compute metric derivatives
    X = reshape(X,Nz*Nx,1);
    Z = reshape(Z,Nz*Nx,1);
    
    Xxi = Dxi*X;
    Xeta = Deta*X;
    Zxi = Dxi*Z;
    Zeta = Deta*Z;
    
    
    
    % Scale factors
%     ETA1 = sqrt(Xxi.^2 + Zxi.^2);
%     ETA2 = sqrt(Xeta.^2 + Zeta.^2);
%     
%     J = ETA1.*ETA2;
%     J2 = Xxi.*Zeta - Zxi.*Xeta;
%     error = abs(J-J2);
%     X = reshape(X,Nz,Nx);
%     Z = reshape(Z,Nz,Nx);
%     error = reshape(error,Nz,Nx);
%     figure
%     surf(X,Z,error)
%     view(0,90);
%     pause
%     norm(J2);
%     orthogonality_error = norm(J-J2)
%     max_orthogonality_error = max(max(abs(J-J2)))
%     rel_error = norm(J-J2)/norm(J2)
%     
%     error = abs(Zxi.*Zeta + Xxi.*Xeta);
%     error = reshape(error,Nz,Nx);
%     figure
%     surf(X,Z,error);
    
    X = reshape(X,Nz,Nx);
    Z = reshape(Z,Nz,Nx);
    Xxi = reshape(Xxi,Nz,Nx);
    Zxi = reshape(Zxi,Nz,Nx);
    Xeta = reshape(Xeta,Nz,Nx);
    Zeta = reshape(Zeta,Nz,Nx);
    

end