@doc raw"""
    myinterpolate3(xgrd1, xgrd2, xgrd3, ygrd, xeval1, xeval2, xeval3)

Trilineary project `ygrd` on (`xgrd1`,`xgrd2`,`xgrd3`) and use it to
interpolate value at (`xeval1`,`xeval2`,`xeval3`).

# Example
```jldoctest
julia> xgrd = [1.0,6.0];
julia> f((x,y,z)) = x+y+z;
julia> ygrd = f.(collect(Iterators.product(xgrid,xgrid,xgrid));
julia> xeval = [3.0,5.0];
julia> mylinearinterpolate3(xgrd,xgrd,xgrd,ygrd,xeval,xeval,xeval)
2x2x2 Array{Float64,3}:
[:,:,1] =
 9.0 11.0
11.0 13.0
[:,:,2] =
11.0 13.0
13.0 15.0
```
"""
myinterpolate3(xgrd1::AbstractVector, xgrd2::AbstractVector,
    xgrd3::AbstractVector, ygrd::AbstractArray,
    xeval1::AbstractVector, xeval2::AbstractVector, xeval3::AbstractVector) = 
    mylinearinterpolate3(xgrd1, xgrd2, xgrd3, ygrd, xeval1, xeval2, xeval3)

# use Akima interpolation 
# myinterpolate3(xgrd1::AbstractVector, xgrd2::AbstractVector,
#     xgrd3::AbstractVector, ygrd::AbstractArray,
#     xeval1::AbstractVector, xeval2::AbstractVector, xeval3::AbstractVector) = 
#                     eval_Akima3(xgrd1, xgrd2, xgrd3,
#                                 poly_coeff_Akima3(xgrd1, xgrd2, xgrd3, F) ,
#                                 xeval1, xeval2, xeval3)
#-------------------------------------------------------------------------------
## Linear interpolations ##
#-------------------------------------------------------------------------------
# of one function
function mylinearinterpolate(xgrd::AbstractVector,
    ygrd::AbstractVector,
    xeval::AbstractVector)
    yeval= Array{eltype(ygrd),1}(undef,length(xeval))
    mylinearinterpolate!(yeval,xgrd,ygrd,xeval)

    return yeval
end

function mylinearinterpolate!(yeval::AbstractVector,xgrd::AbstractVector,
    ygrd::AbstractVector,
    xeval::AbstractVector)

    n_xgrd = length(xgrd)
    #yeval= Array{eltype(ygrd),1}(undef,length(xeval))

    n_xgrd = length(xgrd)
    @views for i in eachindex(xeval)
        xi = xeval[i]
        if xi .> xgrd[n_xgrd-1]
            iL = n_xgrd - 1
        elseif xi .< xgrd[2]
            iL = 1
        else
            iL      = locate(xi, xgrd)
        end
        iR = iL+1
        xL = xgrd[iL]
        wR = (xi .- xL)./ (xgrd[iR] .- xL)
        wL = 1.0-wR
        yeval[i] = wL .* ygrd[iL] .+ wR .* ygrd[iR]
    end

end

# of two functions on the same grids
function mylinearinterpolate_mult2(xgrd::AbstractVector,
    ygrd1::AbstractVector,
    ygrd2::AbstractVector,
    xeval::AbstractVector)
    yeval1= Array{eltype(ygrd1),1}(undef,length(xeval))
    yeval2= Array{eltype(ygrd1),1}(undef,length(xeval))
    mylinearinterpolate_mult2!(yeval1, yeval2, xgrd, ygrd1, ygrd2, xeval)
    return yeval1, yeval2
end

function mylinearinterpolate_mult2!(yeval1::AbstractVector,
    yeval2::AbstractVector,
    xgrd::AbstractVector,
    ygrd1::AbstractVector,
    ygrd2::AbstractVector,
    xeval::AbstractVector)


    n_xgrd = length(xgrd)
    @views @inbounds begin
        for i in eachindex(xeval)
        xi = xeval[i]
        if xi .> xgrd[n_xgrd-1]
            iL = n_xgrd - 1
        elseif xi .< xgrd[2]
            iL = 1
        else
            iL      = locate(xi, xgrd)
        end
        iR = iL+1
        xL = xgrd[iL]
        wR = (xeval[i] .- xL)./ (xgrd[iR] .- xL)
        #wL = 1.0-wR
        y1L = ygrd1[iL]
        y2L = ygrd2[iL]
        yeval1[i] = y1L .+ wR .* (ygrd1[iR] - y1L)
        yeval2[i] = y2L .+ wR .* (ygrd2[iR] - y2L)
    end
end
    
end

# of three functions on the same grids
function mylinearinterpolate_mult3(xgrd::AbstractVector,
    ygrd1::AbstractVector,
    ygrd2::AbstractVector,
    ygrd3::AbstractVector,
    xeval::AbstractVector)

    yeval1= Array{eltype(ygrd1),1}(undef,length(xeval))
    yeval2= Array{eltype(ygrd2),1}(undef,length(xeval))
    yeval3= Array{eltype(ygrd3),1}(undef,length(xeval))
    mylinearinterpolate_mult3!(yeval1, yeval2, yeval3, xgrd, ygrd1, ygrd2, ygrd3, xeval)
    return yeval1, yeval2, yeval3
end
function mylinearinterpolate_mult3!(yeval1::AbstractVector,
    yeval2::AbstractVector,
    yeval3::AbstractVector,
    xgrd::AbstractVector,
    ygrd1::AbstractVector,
    ygrd2::AbstractVector,
    ygrd3::AbstractVector,
    xeval::AbstractVector)
    n_xgrd = length(xgrd)
    for i in eachindex(xeval)
        xi = xeval[i]
        if xi .> xgrd[n_xgrd-1]
            iL = n_xgrd - 1
        elseif xi .< xgrd[2]
            iL = 1
        else
            iL      = locate(xi, xgrd)
        end
        iR = iL+1
        xL = xgrd[iL]
        wR = (xi .- xL)./ (xgrd[iR] .- xL)
        y1L = ygrd1[iL]
        y2L = ygrd2[iL]
        y3L = ygrd3[iL]
        yeval1[i] = y1L .+ wR .* (ygrd1[iR] - y1L)
        yeval2[i] = y2L .+ wR .* (ygrd2[iR] - y2L)
        yeval3[i] = y3L .+ wR .* (ygrd3[iR] - y3L)
    end
end

#---------------- Billinear interpolation -------------------------------------
function mylinearinterpolate2(xgrd1::AbstractVector, xgrd2::AbstractVector, ygrd::AbstractArray,
    xeval1::AbstractVector, xeval2::AbstractVector)

    yeval   = zeros(eltype(xeval1),length(xeval1),length(xeval2))
    n_xgrd1 = length(xgrd1)
    n_xgrd2 = length(xgrd2)
    weight1 = Array{eltype(xeval1),1}(undef,length(xeval1))
    weight2 = Array{eltype(xeval2),1}(undef,length(xeval2))
    ind1    = zeros(Int,length(xeval1))
    ind2    = zeros(Int,length(xeval2))
    @views for i in eachindex(xeval1)
        xi = xeval1[i]
        if xi .> xgrd1[n_xgrd1-1]
            iL = n_xgrd1 - 1
        elseif xi .< xgrd1[2]
            iL = 1
        else
            iL      = locate(xi, xgrd1)
        end
        ind1[i]      = copy(iL)
        weight1[i]   = copy((xi .- xgrd1[iL])./ (xgrd1[iL.+1] .- xgrd1[iL]))
    end

    @views for i in eachindex(xeval2)
        xi = xeval2[i]
        if xi .> xgrd2[n_xgrd2-1]
            iL = n_xgrd2 - 1
        elseif xi .< xgrd2[2]
            iL = 1
        else
            iL      = locate(xi, xgrd2)
        end
        ind2[i]      = copy(iL)
        weight2[i]   = copy((xi .- xgrd2[iL])./ (xgrd2[iL.+1] .- xgrd2[iL]))
    end

    for j in eachindex(xeval2)
        w2R = weight2[j]
        w2L = 1.0-w2R
        for i in eachindex(xeval1)
            w1R = weight1[i]
            w1L = 1.0-w1R
            aux = w2L*(w1L*ygrd[ind1[i],ind2[j]] + w1R*ygrd[ind1[i]+1,ind2[j]]) +
                  w2R*(w1L*ygrd[ind1[i],ind2[j]+1] + w1R*ygrd[ind1[i]+1,ind2[j]+1])
            yeval[i,j] = aux[1]
        end
    end

    return yeval
end

function mylinearinterpolate3(xgrd1::AbstractVector, xgrd2::AbstractVector,
    xgrd3::AbstractVector, ygrd::AbstractArray,
    xeval1::AbstractVector, xeval2::AbstractVector, xeval3::AbstractVector)
    ## Define variables ##
    yeval   = Array{eltype(xeval1),3}(undef,(length(xeval1),length(xeval2),length(xeval3)))
    n_xgrd1 = length(xgrd1)
    n_xgrd2 = length(xgrd2)
    n_xgrd3 = length(xgrd3)
    weight1r = Array{eltype(xeval1),1}(undef,length(xeval1))
    weight2r = Array{eltype(xeval2),1}(undef,length(xeval2))
    weight3r = Array{eltype(xeval3),1}(undef,length(xeval3))
    ind1    = Array{Int,1}(undef,length(xeval1))
    ind2    = Array{Int,1}(undef,length(xeval2))
    ind3    = Array{Int,1}(undef,length(xeval3))
    ## find left element of xeval1 on xgrd1
    for i in eachindex(xeval1)
        xi = xeval1[i]
        if xi .>= xgrd1[n_xgrd1-1]
            iL = n_xgrd1 - 1
        elseif xi .< xgrd1[2]
            iL = 1
        else
            iL      = locate(xi, xgrd1)
        end
        ind1[i]      = copy(iL)
        xL          = xgrd1[iL]
        weight1r[i]  = copy((xi .- xL)./ (xgrd1[iL.+1] .- xL))
    end

    ## find left element of xeval2 on xgrd2
    for i in eachindex(xeval2)
        xi = xeval2[i]
        if xi .>= xgrd2[n_xgrd2-1]
            iL = n_xgrd2 - 1
        elseif xi .< xgrd2[2]
            iL = 1
        else
            iL      = locate(xi, xgrd2)
        end
        ind2[i]      = copy(iL)
        xL           = xgrd2[iL]
        weight2r[i]  = copy((xi .- xL)./ (xgrd2[iL.+1] .- xL))
    end

    ## find left element of xeval3 on xgrd3
    for i in eachindex(xeval3)
        xi = xeval3[i]
        if xi .>= xgrd3[n_xgrd3-1]
            iL = n_xgrd3 - 1
        elseif xi .< xgrd3[2]
            iL = 1
        else
            iL      = locate(xi, xgrd3)
        end
        ind3[i]      = copy(iL)
        xL           = xgrd3[iL]
        weight3r[i]  = copy((xi .- xL)./ (xgrd3[iL.+1] .- xL))
    end
@inbounds  begin
    for l in eachindex(xeval3)
        i3  = ind3[l]
        w3r = weight3r[l]
        w3l = 1.0 - w3r
        for j in eachindex(xeval2)
            i2 =ind2[j]
            w2r = weight2r[j]
            w2l = 1.0 - w2r
            for k in eachindex(xeval1)
                i1  = ind1[k]
                w1r = weight1r[k]
                w1l = 1.0 - w1r
               yeval[k,j,l] =  w1l.*(w2l.*(w3l.*ygrd[i1,i2,i3]     + w3r.*ygrd[i1,i2,i3+1])    +
                                      w2r.*(w3l.*ygrd[i1,i2+1,i3]   + w3r.*ygrd[i1,i2+1,i3+1])) +
                                w1r.*(w2l.*(w3l.*ygrd[i1+1,i2,i3]   + w3r.*ygrd[i1+1,i2,i3+1])  +
                                      w2r.*(w3l.*ygrd[i1+1,i2+1,i3] + w3r.*ygrd[i1+1,i2+1,i3+1]))
            end
        end
    end
end
    return yeval
end




#-------------------------------------------------------------------------
# 1-D Akima interpolation
#-------------------------------------------------------------------------

function poly_coeff_Akima1(xgrd1::Vector, F::AbstractVector)
    # determine polynomial coefficients for local polynomials using Akima's fomrula
    T       = eltype(xgrd1[1]+F[1]) # Makes sure that this works with forwardDiff, because type of ForwardDiffDual is broken (not AbstractFloat type)
            
    n = length(xgrd1) 
        coeffs = Array{T}(undef,4,n-1)
        t = Array{T}(undef,n,1)
        m = Array{T}(undef,n+3,1)
        
        dx  = diff(xgrd1)
        dx2 = dx.^2
        dx3 = dx.*dx2
        # First order divided differences x-dim
        m[3:n+1] = diff(F)./dx
        # Extrapolation
        m[n+2]   = 2.0 .* m[n+1] - m[n]
        m[n+3]   = 2.0 .* m[n+2] - m[n+1]
        m[2]     = 2.0 .* m[3] - m[4]
        m[1]     = 2.0 .* m[2] - m[3]
        w = abs.(diff(m;dims=1)).+10*eps()
        t = (w[3:n+2].*m[2:n+1] + w[1:n].*m[3:n+2])./(w[3:n+2] + w[1:n])
    
        D = Matrix{T}(I,4,4)
        @inbounds for i = 1:n-1
            D[3:4,:] .=[
                1.0 dx[i] dx2[i] dx3[i];   
                0.0 1.0 2*dx[i] 3*dx2[i]
                ]
            coeffs[:,i] = D\ [F[i]; t[i]; F[i+1]; t[i+1]] 
        end
        return coeffs
    end
    
    function eval_Akima1(xgrd1::AbstractVector,  coeffs::AbstractArray, xeval::AbstractVector)
        # Evaluate local polynomials on a grid of evaluation points
        T       = eltype(xgrd1[1]+xeval[1]+coeffs[1])
        n_xgrd = length(xgrd1)
        yeval = Array{T,1}(undef,length(xeval))
        @views for i in eachindex(xeval)
            xi = xeval[i]
            if xi .> xgrd1[n_xgrd-1]
                iL   = n_xgrd - 1
            elseif xi .< xgrd1[2]
                iL   = 1
            else
                iL   = locate(xi, xgrd1)
            end
            xd        = xi .- xgrd1[iL]
            yeval[i] = dot(coeffs[:,iL],[1 xd xd.^2 xd.^3])
        end
        return yeval
    end
    #-------------------------------------------------------------------------
    # 2-D Akima interpolation
    #-------------------------------------------------------------------------
    function poly_coeff_Akima2(xgrd1::Vector, xgrd2::Vector, F::AbstractMatrix)
        T       = eltype(xgrd1[1]+xgrd2[1]+F[1]) # Makes sure that this works with forwardDiff, because type of ForwardDiffDual is broken (not AbstractFloat type)
        
        # Determine polynomial coefficients
        nx      = length(xgrd1)
        ny      = length(xgrd2)
        coeffs  = Array{T}(undef,16, (nx-1),(ny-1))
    
        mx      = Array{T}(undef,nx+3,ny)
        my      = Array{T}(undef,nx,ny+3)
        mxy     = zeros(T,nx+3,ny+3)
    
        tx      = Array{T}(undef,nx,ny)
        ty      = Array{T}(undef,nx,ny)
        txy     = Array{T}(undef,nx,ny)
    
        dx      = diff(xgrd1)
        dx2     = dx.^2
        dx3     = dx.*dx2
    
        dy      = diff(xgrd2)
        dy2     = dy.^2
        dy3     = dy.*dy2
    
        # First order divided differences x-dim
        mx[3:nx+1,:] = diff(F;dims=1)./dx
        # Extrapolation
        mx[nx+2,:]  = 2.0 .* mx[nx+1,:] .- mx[nx,:]
        mx[nx+3,:]  = 2.0 .* mx[nx+2,:] .- mx[nx+1,:]
        mx[2,:]     = 2.0 .* mx[3,:]    .- mx[4,:]
        mx[1,:]     = 2.0 .* mx[2,:]    .- mx[3,:]
    
        # First order divided differences y-dim
        my[:,3:ny+1] = diff(F;dims=2)./dy'
        # Extrapolation
        my[:,ny+2]  = 2.0 .* my[:,ny+1] .- my[:,ny]
        my[:,ny+3]  = 2.0 .* my[:,ny+2] .- my[:,ny+1]
        my[:,2]     = 2.0 .* my[:,3]    .- my[:,4]
        my[:,1]     = 2.0 .* my[:,2]    .- my[:,3]
    
        # Second order divided differences (croos derivative)
        # Assume cross derivative is zero outside grid
        mxy[3:nx+1,:] = diff(my;dims=1)./dx
        mxy[:,3:ny+1] = diff(mx;dims=2)./dy'
        
    
        #calculate weights 
        wx = abs.(diff(mx;dims=1)).+10*eps()
        wy = abs.(diff(my;dims=2)).+10*eps()
        # calculate tangential values and cross derivatives
        tx = (wx[3:nx+2,:].*mx[2:nx+1,:] + wx[1:nx,:].*mx[3:nx+2,:])./(wx[3:nx+2,:] + wx[1:nx,:])
        ty = (wy[:,3:ny+2].*my[:,2:ny+1] + wy[:,1:ny].*my[:,3:ny+2])./(wy[:,3:ny+2] + wy[:,1:ny])
    
        txy = (wx[3:nx+2,:].*(wy[:,3:ny+2].*mxy[2:nx+1,2:ny+1] .+ wy[:,1:ny].*mxy[2:nx+1,3:ny+2]).+
               wx[1:nx,:]  .*(wy[:,3:ny+2].*mxy[3:nx+2,2:ny+1] .+ wy[:,1:ny].*mxy[3:nx+2,3:ny+2])) ./
             ((wy[:,3:ny+2] .+ wy[:,1:ny]).*(wx[3:nx+2,:] .+ wx[1:nx,:]))
    
        Dx = Matrix{T}(I,4,4)
        Dy = Matrix{T}(I,4,4)
        for i = 1:nx-1
            Dx[3,:] = [1.0 dx[i] dx2[i] dx3[i]]   
            Dx[4,:] = [0.0 1.0 2*dx[i] 3*dx2[i]]
            for j=1:ny-1
                Dy[3,:] = [1.0 dy[j] dy2[j]  dy3[j]]   
                Dy[4,:] = [0.0 1.0 2*dy[j] 3*dy2[j]]
                coeffs[:,i,j] = Dy\ [F[i,j]     tx[i,j]     F[i+1,j]      tx[i+1,j]; 
                                    ty[i,j]     txy[i,j]    ty[i+1,j]     txy[i+1,j];
                                    F[i,j+1]    tx[i,j+1]   F[i+1,j+1]    tx[i+1,j+1];    
                                    ty[i,j+1]   txy[i,j+1]  ty[i+1,j+1]   txy[i+1,j+1]]/(Dx')
            end
        end
        return coeffs
    end
    
    function eval_Akima2(xgrd1::AbstractVector, xgrd2::AbstractVector, coeffs::AbstractArray,
        xeval1::AbstractVector, xeval2::AbstractVector)
        T       = eltype(xgrd1[1]+xgrd2[1]+xeval1[1]+xeval2[1]+coeffs[1])
        # Evaluate local polynomials on a mesh of evaluation points
        yeval   = zeros(T,length(xeval1),length(xeval2))
        n_xgrd1 = length(xgrd1)
        n_xgrd2 = length(xgrd2)
        xd1     = Array{T,1}(undef,length(xeval1))
        xd2     = Array{T,1}(undef,length(xeval2))
        ind1    = zeros(Int,length(xeval1))
        ind2    = zeros(Int,length(xeval2))
       
        @views for i in eachindex(xeval1)
            xi = xeval1[i]
            if xi .> xgrd1[n_xgrd1-1]
                iL  = n_xgrd1 - 1
            elseif xi .< xgrd1[2]
                iL  = 1
            else
                iL  = locate(xi, xgrd1)
            end
            xd1[i]  = xi .- xgrd1[iL]
            ind1[i] = copy(iL)
        end
        @views for i in eachindex(xeval2)
            xi = xeval2[i]
            if xi .> xgrd2[n_xgrd2-1]
                iL  = n_xgrd2 - 1
            elseif xi .< xgrd2[2]
                iL  = 1
            else
                iL  = locate(xi, xgrd2)
            end
            ind2[i] = copy(iL)
            xd2[i]  = xi .- xgrd2[iL]
        end
        # evaluate polynomial
        for j in eachindex(xeval2)
            for i in eachindex(xeval1)
                c     = coeffs[:,ind1[i],ind2[j]][:]
                count = 0
                for h = 0:3
                    aux_x = xd1[i].^h
                    for k = 0:3
                        count+=1
                        yeval[i,j] += c[count] .* aux_x .* xd2[j].^k
                    end
                end 
            end
        end
        return yeval
    end
    #-------------------------------------------------------------------------
    # 3-D Akima interpolation, following Dubey and Upadhyay (1989, Computer Physics Communications)
    #-------------------------------------------------------------------------
    
    myAkimaInterp3(xgrd1::AbstractVector, xgrd2::AbstractVector,
    xgrd3::AbstractVector, F::AbstractArray,
    xeval1::AbstractVector, xeval2::AbstractVector, xeval3::AbstractVector) = 
                        eval_Akima3(xgrd1, xgrd2, xgrd3,
                                    poly_coeff_Akima3(xgrd1, xgrd2, xgrd3, F) ,
                                    xeval1, xeval2, xeval3)
    
    function poly_coeff_Akima3(xgrd1::Vector, xgrd2::Vector, xgrd3::Vector, F::AbstractArray)
        T       = eltype(xgrd1[1]+xgrd2[1]+xgrd3[1]+F[1]) # Makes sure that this works with forwardDiff, because type of ForwardDiffDual is broken (not AbstractFloat type)
        nx      = length(xgrd1)
        ny      = length(xgrd2)
        nz      = length(xgrd3)
        coeffs  = Array{T}(undef,64,(nx-1),(ny-1),(nz-1))

        mx      = Array{T}(undef,nx+3,ny,nz)
        my      = Array{T}(undef,nx,ny+3,nz)
        mz      = Array{T}(undef,nx,ny,nz+3)
    
        mxy     = zeros(T,nx+3,ny+3,nz)
        mxz     = zeros(T,nx+3,ny,nz+3)
        myz     = zeros(T,nx,ny+3,nz+3)
    
        mxyz     = zeros(T,nx+3,ny+3,nz+3)
    
        tx      = Array{T}(undef,nx,ny)
        ty      = Array{T}(undef,nx,ny)
        tz      = Array{T}(undef,nx,ny)
        txy     = Array{T}(undef,nx,ny)
        txz     = Array{T}(undef,nx,ny)
        tyz     = Array{T}(undef,nx,ny)
        txyz    = Array{T}(undef,nx,ny)
    
        dx      = diff(xgrd1)
        dx2     = dx.^2
        dx3     = dx.*dx2
    
        dy      = diff(xgrd2)
        dy2     = dy.^2
        dy3     = dy.*dy2
    
        dz      = diff(xgrd3)
        dz2     = dz.^2
        dz3     = dz.*dz2
    
        # First order divided differences x-dim
        mx[3:nx+1,:,:] = diff(F;dims=1)./dx
        # Extrapolation
        mx[nx+2,:,:]  = 2.0 .* mx[nx+1,:,:] .- mx[nx,:,:]
        mx[nx+3,:,:]  = 2.0 .* mx[nx+2,:,:] .- mx[nx+1,:,:]
        mx[2,:,:]     = 2.0 .* mx[3,:,:]    .- mx[4,:,:]
        mx[1,:,:]     = 2.0 .* mx[2,:,:]    .- mx[3,:,:]
    
        # First order divided differences y-dim
        my[:,3:ny+1,:] = diff(F;dims=2)./dy'
        # Extrapolation
        my[:,ny+2,:]  = 2.0 .* my[:,ny+1,:] .- my[:,ny,:]
        my[:,ny+3,:]  = 2.0 .* my[:,ny+2,:] .- my[:,ny+1,:]
        my[:,2,:]     = 2.0 .* my[:,3,:]    .- my[:,4,:]
        my[:,1,:]     = 2.0 .* my[:,2,:]    .- my[:,3,:]
    
        # First order divided differences z-dim
        mz[:,:,3:nz+1] = diff(F;dims=3)./reshape(dz,(1,1,nz-1))
        # Extrapolation
        mz[:,:,nz+2]  = 2.0 .* mz[:,:,nz+1] .- mz[:,:,nz]
        mz[:,:,nz+3]  = 2.0 .* mz[:,:,nz+2] .- mz[:,:,nz+1]
        mz[:,:,2]     = 2.0 .* mz[:,:,3]    .- mz[:,:,4]
        mz[:,:,1]     = 2.0 .* mz[:,:,2]    .- mz[:,:,3]
    
        # Second order divided differences (first cross derivative)
        # Assume cross derivative is zero outside grid 
        # (innefficient many calculations, when filling the gaps)
        mxy[3:nx+1,:,:] = diff(my;dims=1)./dx
        mxy[:,3:ny+1,:] = diff(mx;dims=2)./dy'
    
        mxz[3:nx+1,:,:] = diff(mz;dims=1)./dx
        mxz[:,:,3:nz+1] = diff(mx;dims=3)./reshape(dz,(1,1,nz-1))
        
        myz[:,3:ny+1,:] = diff(mz;dims=2)./dy'
        myz[:,:,3:nz+1] = diff(my;dims=3)./reshape(dz,(1,1,nz-1))
    
        # Third order divided differences (first cross derivative)
        mxyz[3:nx+1,:,:] = diff(myz;dims=1)./dx
        mxyz[:,3:ny+1,:] = diff(mxz;dims=2)./dy'
        mxyz[:,:,3:nz+1] = diff(mxy;dims=3)./reshape(dz,(1,1,nz-1))
    
        #calculate weights 
        wx = abs.(diff(mx;dims=1)).+10*eps()
        wy = abs.(diff(my;dims=2)).+10*eps()
        wz = abs.(diff(mz;dims=3)).+10*eps()
        # calculate tangential values 
        tx = (wx[3:nx+2,:,:].*mx[2:nx+1,:,:] + wx[1:nx,:,:].*mx[3:nx+2,:,:])./(wx[3:nx+2,:,:] + wx[1:nx,:,:])
        ty = (wy[:,3:ny+2,:].*my[:,2:ny+1,:] + wy[:,1:ny,:].*my[:,3:ny+2,:])./(wy[:,3:ny+2,:] + wy[:,1:ny,:])
        tz = (wz[:,:,3:nz+2].*mz[:,:,2:nz+1] + wz[:,:,1:nz].*mz[:,:,3:nz+2])./(wz[:,:,3:nz+2] + wz[:,:,1:nz])
    
        # and second order (i.e., cross) derivatives
        txy = (wx[3:nx+2,:,:].*(wy[:,3:ny+2,:].*mxy[2:nx+1,2:ny+1,:] .+ wy[:,1:ny,:].*mxy[2:nx+1,3:ny+2,:]).+
               wx[1:nx,:,:]  .*(wy[:,3:ny+2,:].*mxy[3:nx+2,2:ny+1,:] .+ wy[:,1:ny,:].*mxy[3:nx+2,3:ny+2,:])) ./
             ((wy[:,3:ny+2,:] .+ wy[:,1:ny,:]).*(wx[3:nx+2,:,:] .+ wx[1:nx,:,:]))
    
        txz = (wx[3:nx+2,:,:].*(wz[:,:,3:nz+2].*mxz[2:nx+1,:,2:nz+1] .+ wz[:,:,1:nz].*mxz[2:nx+1,:,3:nz+2]).+
               wx[1:nx,:,:]  .*(wz[:,:,3:nz+2].*mxz[3:nx+2,:,2:nz+1] .+ wz[:,:,1:nz].*mxz[3:nx+2,:,3:nz+2])) ./
               ((wz[:,:,3:nz+2] .+ wz[:,:,1:nz]).*(wx[3:nx+2,:,:] .+ wx[1:nx,:,:]))
        
        tyz = (wz[:,:,3:nz+2].*(wy[:,3:ny+2,:].*myz[:,2:ny+1,2:nz+1] .+ wy[:,1:ny,:].*myz[:,3:ny+2,2:nz+1]).+
               wz[:,:,1:nz]  .*(wy[:,3:ny+2,:].*myz[:,2:ny+1,3:nz+2] .+ wy[:,1:ny,:].*myz[:,3:ny+2,3:nz+2])) ./
             ((wy[:,3:ny+2,:] .+ wy[:,1:ny,:]).*(wz[:,:,3:nz+2] .+ wz[:,:,1:nz]))
        # and third order (i.e., tripple cross) derivatives
        
        txyz = (wx[3:nx+2,:,:].*(wy[:,3:ny+2,:].*(wz[:,:,3:nz+2].*mxyz[2:nx+1,2:ny+1,2:nz+1] .+
                                                  wz[:,:,1:nz]  .*mxyz[2:nx+1,2:ny+1,3:nz+2] ) .+ 
                                 wy[:,1:ny,:]  .*(wz[:,:,3:nz+2].*mxyz[2:nx+1,3:ny+2,2:nz+1] .+
                                                  wz[:,:,1:nz]  .*mxyz[2:nx+1,3:ny+2,3:nz+2] )).+
                wx[1:nx,:,:]  .*(wy[:,3:ny+2,:].*(wz[:,:,3:nz+2].*mxyz[3:nx+2,2:ny+1,2:nz+1] .+
                                                  wz[:,:,1:nz]  .*mxyz[3:nx+2,2:ny+1,3:nz+2] ) .+ 
                                 wy[:,1:ny,:]  .*(wz[:,:,3:nz+2].*mxyz[3:nx+2,3:ny+2,2:nz+1] .+
                                                  wz[:,:,1:nz]  .*mxyz[3:nx+2,3:ny+2,3:nz+2] ))) ./
                ((wy[:,3:ny+2,:] .+ wy[:,1:ny,:]).*(wx[3:nx+2,:,:] .+ wx[1:nx,:,:]).*(wz[:,:,3:nz+2] .+ wz[:,:,1:nz]))
        # this it what costs time can we do this better? (i.e. vectorize instead of )
        Dx = Matrix{T}(I,4,4)
        Dy = Matrix{T}(I,4,4)
        Dz = Matrix{T}(I,4,4)
        DxInv=Vector{Matrix{T}}(undef,nx-1)
        DyInv=Vector{Matrix{T}}(undef,ny-1)
        DzInv=Vector{Matrix{T}}(undef,nz-1)
        @inbounds for k = 1:nz-1
            Dz[3,:] = [1.0 dz[k] dz2[k]  dz3[k]]   
            Dz[4,:] = [0.0 1.0 2*dz[k] 3*dz2[k]]
            DzInv[k]= Dz \ I
        end
        @inbounds for j=1:ny-1
            Dy[3,:] = [1.0 dy[j] dy2[j]   dy3[j]]   
            Dy[4,:] = [0.0 1.0  2*dy[j] 3*dy2[j]]
            DyInv[j] = Dy \ I
        end
        @inbounds for i = 1:nx-1
            Dx[3,:] = [1.0 dx[i] dx2[i] dx3[i]]   
            Dx[4,:] = [0.0 1.0 2*dx[i] 3*dx2[i]]    
            DxInv[i] = I / Dx'
        end
        A = Array{T,5}(undef,(8,2,nx,ny,nz-1))
        A[1,1,:,:,:] .= F[1:nx,  1:ny,  1:nz-1]
        A[2,1,:,:,:] .= tz[1:nx, 1:ny,  1:nz-1]    
        A[3,1,:,:,:] .= F[1:nx,  1:ny,  2:nz]  
        A[4,1,:,:,:] .= tz[1:nx, 1:ny,  2:nz] 
        A[5,1,:,:,:] .= ty[1:nx, 1:ny,  1:nz-1]  
        A[6,1,:,:,:] .= tyz[1:nx,1:ny,  1:nz-1]  
        A[7,1,:,:,:] .= ty[1:nx, 1:ny,  2:nz]  
        A[8,1,:,:,:] .= tyz[1:nx,1:ny,  2:nz]

        A[1,2,:,:,:] .= tx[1:nx,  1:ny,  1:nz-1]
        A[2,2,:,:,:] .= txz[1:nx, 1:ny,  1:nz-1]
        A[3,2,:,:,:] .= tx[1:nx,  1:ny,  2:nz]     
        A[4,2,:,:,:] .= txz[1:nx, 1:ny,  2:nz]     
        A[5,2,:,:,:] .= txy[1:nx, 1:ny,  1:nz-1]  
        A[6,2,:,:,:] .= txyz[1:nx,1:ny,  1:nz-1] 
        A[7,2,:,:,:] .= txy[1:nx, 1:ny,  2:nz]    
        A[8,2,:,:,:] .= txyz[1:nx,1:ny,  2:nz]   
    
         for k = 1:nz-1
             for j=1:ny-1
                DzDyInv =(kron(DyInv[j],DzInv[k]))
                B = [ A[:,:,:,j,k];
                      A[:,:,:,j+1,k] ]
                B = reshape(DzDyInv*reshape(B,16,2*nx),(16,2,nx))
                 @views for i = 1:nx-1
                    coeffs[:,i,j,k] = ([B[:,:,i] B[:,:,i+1]] * DxInv[i])[:]
                    # AX = [
                    # F[i,  j,  k]      tx[i,j,k]       F[i+1,j,k]       tx[i+1,j,k];
                    # tz[i, j,  k]      txz[i,j,k]      tz[i+1,j,k]      txz[i+1,j,k];
                    # F[i,  j,  k+1]    tx[i,j,k+1]     F[i+1,j,k+1]     tx[i+1,j,k+1];
                    # tz[i, j,  k+1]    txz[i,j,k+1]    tz[i+1,j,k+1]    txz[i+1,j,k+1]; 
                    # ty[i, j,  k]      txy[i,j,k]      ty[i+1,j,k]      txy[i+1,j,k];
                    # tyz[i,j,  k]      txyz[i,j,k]     tyz[i+1,j,k]     txyz[i+1,j,k];
                    # ty[i, j,  k+1]    txy[i,j,k+1]    ty[i+1,j,k+1]    txy[i+1,j,k+1];
                    # tyz[i,j,  k+1]    txyz[i,j,k+1]   tyz[i+1,j,k+1]   txyz[i+1,j,k+1];
                    # F[i,  j+1,k]      tx[i,j+1,k]     F[i+1,j+1,k]     tx[i+1,j+1,k];
                    # tz[i, j+1,k]      txz[i,j+1,k]    tz[i+1,j+1,k]    txz[i+1,j+1,k];
                    # F[i,  j+1,k+1]    tx[i,j+1,k+1]   F[i+1,j+1,k+1]   tx[i+1,j+1,k+1];
                    # tz[i, j+1,k+1]    txz[i,j+1,k+1]  tz[i+1,j+1,k+1]  txz[i+1,j+1,k+1];     
                    # ty[i, j+1,k]      txy[i,j+1,k]    ty[i+1,j+1,k]    txy[i+1,j+1,k];
                    # tyz[i,j+1,k]      txyz[i,j,k]     tyz[i+1,j,k]     txyz[i+1,j,k];
                    # ty[i, j+1,k+1]    txy[i,j+1,k+1]  ty[i+1,j+1,k+1]  txy[i+1,j+1,k+1];
                    # tyz[i,j+1,k+1]    txyz[i,j+1,k+1] tyz[i+1,j+1,k+1] txyz[i+1,j+1,k+1]
                    # ]
                end
            end
        end
  
        return coeffs
    end
    
    function eval_Akima3(xgrd1::AbstractVector, xgrd2::AbstractVector, xgrd3::AbstractVector, coeffs::AbstractArray,
        xeval1::AbstractVector, xeval2::AbstractVector, xeval3::AbstractVector)
        # Evaluate local polynomials on a mesh of evaluation points
        T       = eltype(xgrd1[1]+xgrd2[1]+xgrd3[1]+xeval1[1]+xeval2[1]+xeval3[1]+coeffs[1])
        ## Define variables ##
        yeval   = zeros(T,(length(xeval1),length(xeval2),length(xeval3)))
        n_xgrd1 = length(xgrd1)
        n_xgrd2 = length(xgrd2)
        n_xgrd3 = length(xgrd3)
        pp1     = Array{T,2}(undef,(length(xeval1),4))
        pp2     = Array{T,2}(undef,(length(xeval2),4))
        pp3     = Array{T,2}(undef,(length(xeval3),4))
        ind1    = Array{Int,1}(undef,length(xeval1))
        ind2    = Array{Int,1}(undef,length(xeval2))
        ind3    = Array{Int,1}(undef,length(xeval3))
       
        @views for i in eachindex(xeval1)
            xi = xeval1[i]
            if xi .> xgrd1[n_xgrd1-1]
                iL  = n_xgrd1 - 1
            elseif xi .< xgrd1[2]
                iL  = 1
            else
                iL  = locate(xi, xgrd1)
            end
            xd1  = xi .- xgrd1[iL]
            ind1[i] = copy(iL)
            pp1[i,:]= [1 xd1 xd1.^2 xd1.^3]
        end
        @views for i in eachindex(xeval2)
            xi = xeval2[i]
            if xi .> xgrd2[n_xgrd2-1]
                iL  = n_xgrd2 - 1
            elseif xi .< xgrd2[2]
                iL  = 1
            else
                iL  = locate(xi, xgrd2)
            end
            ind2[i] = copy(iL)
            xd2  = xi .- xgrd2[iL]
            pp2[i,:]= [1 xd2 xd2.^2 xd2.^3]
        end
        @views for i in eachindex(xeval3)
            xi = xeval3[i]
            if xi .> xgrd3[n_xgrd3-1]
                iL  = n_xgrd3 - 1
            elseif xi .< xgrd3[2]
                iL  = 1
            else
                iL  = locate(xi, xgrd3)
            end
            ind3[i] = copy(iL)
            xd3  = xi .- xgrd3[iL]
            pp3[i,:]= [1 xd3 xd3.^2 xd3.^3]
        end
        # evaluate polynomial
        
        @inbounds @views for k in eachindex(xeval3) 
            slice1 = reshape(pp3[k,:]'*reshape(coeffs[:,:,:,ind3[k]],(4,16*(n_xgrd1-1)*(n_xgrd2-1))),16,(n_xgrd1-1),(n_xgrd2-1))
            @inbounds @views for j in eachindex(xeval2)
                slice2= reshape(pp2[j,:]'*reshape(slice1[:,:,ind2[j]],(4,4*(n_xgrd1-1))),(4,n_xgrd1-1))
                @inbounds @views for i in eachindex(xeval1)
                    yeval[i,j,k] = pp1[i,:]'*slice2[:,ind1[i]]
                end
            end
        end
        return yeval
    end
    