"""
    interpolation2Dsym(r::Vector{Float64}, z::Vector{Float64}, v::Vector{Float64}, Rq::Vector{Float64}, Zq::Vector{Float64})

Interpolates the values given in `v` at the points `r, z` and interpolates
them at the query points 'Rq, Zq' and returns `Vq`.

`v`, `r`, and `z` are vectors of size N×1.
`Vq`, `Rq`, and `Rz` are vectors of size M×1.

"""
function interpolation2Dsym(r::Vector{Float64}, z::Vector{Float64}, v::Vector{Float64}, Rq::Vector{Float64}, Zq::Vector{Float64})
    Vq = zeros(Float64, length(Rq)) # init result vector

    # loop over every query point
    for i in 1:length(Rq)
        rq = Rq[i]; zq = Zq[i]; # extract current point
        corners = _findcorners(r, z, rq, zq); # get corners of interpolation quad

        # check if all corners were found
        if sum(corners .== 0) == 0
            X = _X(r, z, corners); # get matrix representation of corners
            V = [v[corners[1]]; v[corners[2]]; v[corners[3]]; v[corners[4]]]; # get vector of corner values
            Vq[i] = _interppoint(X, V, rq, zq); # add interpolated value to vector
        else
            # if no interpolation quad could be found
            Vq[i] = NaN;
        end
    end

    return Vq
end
interpolation2Dsym(r::Vector{Float64}, z::Vector{Float64}, v::Vector{Float64}, Rq::LinRange{Float64}, Zq::LinRange{Float64}) = interpolation2Dsym(r,   z,   v, collect(Rq), collect(Zq))
interpolation2Dsym(r::Vector{Float64}, z::Vector{Float64}, v::Vector{Float64}, Rq::LinRange{Float64}, Zq::Vector{Float64})   = interpolation2Dsym(r,   z,   v, collect(Rq), Zq)
interpolation2Dsym(r::Vector{Float64}, z::Vector{Float64}, v::Vector{Float64}, Rq::Vector{Float64},   Zq::LinRange{Float64}) = interpolation2Dsym(r,   z,   v, Rq,          collect(Zq))
interpolation2Dsym(p::CPosition2Dsym,                      v::Vector{Float64}, Rq::Vector{Float64},   Zq::Vector{Float64})   = interpolation2Dsym(p.r, p.z, v, Rq,          Zq)
interpolation2Dsym(p::CPosition2Dsym,                      v::Vector{Float64}, Rq::LinRange{Float64}, Zq::LinRange{Float64}) = interpolation2Dsym(p.r, p.z, v, collect(Rq), collect(Zq))
interpolation2Dsym(p::CPosition2Dsym,                      v::Vector{Float64}, Rq::LinRange{Float64}, Zq::Vector{Float64})   = interpolation2Dsym(p.r, p.z, v, collect(Rq), Zq)
interpolation2Dsym(p::CPosition2Dsym,                      v::Vector{Float64}, Rq::Vector{Float64},   Zq::LinRange{Float64}) = interpolation2Dsym(p.r, p.z, v, Rq,          collect(Zq))


"""
Finds the closest corner points of the enclosing quad element.
"""
function _findcorners(r::Vector{Float64}, z::Vector{Float64}, rq::Float64, zq::Float64)
    d = sqrt.( (rq .- r).^2 + (zq .- z).^2 );
    corners = zeros(Int64, 4);
    for i in 1:length(d)
        r[i]<rq && z[i]<zq ? (corners[1]==0 || d[i]<d[corners[1]] ? corners[1] = i : continue) :
        r[i]>rq && z[i]<zq ? (corners[2]==0 || d[i]<d[corners[2]] ? corners[2] = i : continue) :
        r[i]>rq && z[i]>zq ? (corners[3]==0 || d[i]<d[corners[3]] ? corners[3] = i : continue) :
        r[i]<rq && z[i]>zq ? (corners[4]==0 || d[i]<d[corners[4]] ? corners[4] = i : continue) : continue
    end
    return corners
end


"""
Interpolates at the query point `rq, zq` given the matrix representation of the
corner coordinates `X` (16×2) and the vector representation of the corner values
`V` (4×1).
"""
function _interppoint(X::Matrix{Float64}, V::Vector{Float64}, rq::Float64, zq::Float64; λ=0.1, rtol=1e-4, N=250)
    ξi = [0.0; 0.0] # initial guess in natural coordinates

    # loop to find natural coordinates corresponding to [rq; zq]
    counter = 1
    while true
        counter += 1; if counter > N; break; end

        Ni = _N(ξi) # shape functions
        xi = Ni * X[1:2:end, 1] # coordinates of guess ξi
        dNdξi = _dNdξ(ξi) # derivative of shape functions wrt. natural coordinates
        dxdξi = dNdξi * X # derivative of coordinates wrt. natural coordinates
        dξdxi = inv(dxdξi) # derivative of natural coordinates wrt. coordinates

        # update natural coordinate guess
        ξi += dξdxi * ([rq; zq] - xi) * λ

        # break condition
        sqrt((rq-xi[1])^2 + (zq-xi[2])^2) < rtol ? break : continue
    end

    # calculate interpolated value with shape functions
    N = transpose(_N(ξi)[1, 1:2:end])
    return N * V
end


"""
Generates the matrix representation of the corner coordinates `X` (16×2).
"""
function _X(r::Vector{Float64}, z::Vector{Float64}, corners::Vector{Int64})
    return [
        r[corners[1]] 0;
        0 r[corners[1]];
        z[corners[1]] 0;
        0 z[corners[1]];

        r[corners[2]] 0;
        0 r[corners[2]];
        z[corners[2]] 0;
        0 z[corners[2]];

        r[corners[3]] 0;
        0 r[corners[3]];
        z[corners[3]] 0;
        0 z[corners[3]];

        r[corners[4]] 0;
        0 r[corners[4]];
        z[corners[4]] 0;
        0 z[corners[4]];
    ]
end


"""
Generates the shape functions for the mapping of local coordinates (r, z) to
natural coordinates (ξ₁, ξ₂) as a matrix `N` (2×8).
"""
function _N(ξ::Vector{Float64})
    N1 = 1/4 * (1 - ξ[1]) * (1 - ξ[2]);
    N2 = 1/4 * (1 + ξ[1]) * (1 - ξ[2]);
    N3 = 1/4 * (1 + ξ[1]) * (1 + ξ[2]);
    N4 = 1/4 * (1 - ξ[1]) * (1 + ξ[2]);
    return [N1 0 N2 0 N3 0 N4 0; 0 N1 0 N2 0 N3 0 N4]
end


"""
Generates the matrix representation of the gradients of the shape functions wrt.
the natural coordinates (ξ₁, ξ₂) as `dNdξ` (2×16).
"""
function _dNdξ(ξ::Vector{Float64})
    dN1dξ1 = - 1/4 * (1 - ξ[2]); dN1dξ2 = - 1/4 * (1 - ξ[1]);
    dN2dξ1 = + 1/4 * (1 - ξ[2]); dN2dξ2 = - 1/4 * (1 - ξ[1]);
    dN3dξ1 = + 1/4 * (1 - ξ[2]); dN3dξ2 = + 1/4 * (1 - ξ[1]);
    dN4dξ1 = - 1/4 * (1 - ξ[2]); dN4dξ2 = + 1/4 * (1 - ξ[1]);

    return [
        dN1dξ1 dN1dξ2 0 0 dN2dξ1 dN2dξ2 0 0 dN3dξ1 dN3dξ2 0 0 dN4dξ1 dN4dξ2 0 0;
        0 0 dN1dξ1 dN1dξ2 0 0 dN2dξ1 dN2dξ2 0 0 dN3dξ1 dN3dξ2 0 0 dN4dξ1 dN4dξ2
    ]
end
