function [lambda,lambdagrad,lambdacurl] = ...
  Plygn_Wachspress_ValGradCurl(vrtx,pt) 
%% Polygon: Wachspress coordinates: value, gradient, curl 
% James Liu, Zhuoran Wang, ColoState; 2012/07--2017/08 

%% Polygon info 
m = size(vrtx,1);

%% Edge lengths and outward unit normal vectors 
LenEg = zeros(m,1);
NmlEg = zeros(m,2);
for j=1:m
  j1 = mod(j,m) + 1;
  tmpx = vrtx(j1,1) - vrtx(j,1);
  tmpy = vrtx(j1,2) - vrtx(j,2);
  LenEg(j) = sqrt(tmpx.^2+tmpy.^2);
  NmlEg(j,1) =  tmpy./LenEg(j);
  NmlEg(j,2) = -tmpx./LenEg(j);
end

%% Elementwise Wachspress coordinates: value, gradient, curl  
% Distances to the edges and scaled normal vectors 
p = zeros(m,2);
for j=1:m
  dist = dot(vrtx(j,:)-pt, NmlEg(j,:));
  p(j,1) = NmlEg(j,1)/dist;
  p(j,2) = NmlEg(j,2)/dist;
end
% More 
w = zeros(m,1);
R = zeros(m,2);
for j=1:m
  k = mod(j-2,m) + 1;
  w(j) = p(k,1).*p(j,2) - p(k,2).*p(j,1);  % Determinant 
  R(j,:) = p(k,:) + p(j,:); 
end
wsum = sum(w);
% Wachspress coordinates: value 
lambda = w/wsum;
% Wachspress coordinates: gradient 
lambdagrad = zeros(m,2);
for j=1:m
  lambdagrad(j,1) = lambda(j) * ( R(j,1) - dot(lambda,R(:,1)) );
  lambdagrad(j,2) = lambda(j) * ( R(j,2) - dot(lambda,R(:,2)) );
end
% Wachspress coordinates: curl 
lambdacurl = zeros(m,2);
for j=1:m
  lambdacurl(j,1) = -lambdagrad(j,2);  % Here to put the neg.sign
  lambdacurl(j,2) =  lambdagrad(j,1);  
end

return;