function w_map = Warpfactor(N,coordinates)

% compute scaled warp function at order N based on coordinates interpolation nodes

r_GL = GaussLobattoNodes(N + 1);
r_eq = linspace(-1,1,N + 1);

LagPols = LagrangePolynomialBasis(r_eq);

N_p = length(coordinates);
LagPolsValues = zeros(N + 1,N_p);
for i = 1:(N + 1)
    LagPolsValues(i,:) = polyval(LagPols{i},coordinates); % in row vector of values of Lag. polynomial evalueated in coordinates
end

w_map = (r_GL - r_eq) * LagPolsValues;

sing = find(abs(coordinates) == 1); % singularities in blend functions for vertices
reg = setdiff(1:N_p,sing); % points where we can scale w(r)
w_map(reg) = w_map(reg) ./ (1 - (coordinates(reg)) .^ 2);
w_map(sing) = 0;

end