function [vq] = TriInterp(x,y,v,xq,yq,ratio)
%Input arguments
%x: lon
%y: lat
%v: input values
%xq: lon of query points
%yq: lat of query points
%ratio: ration of interval between lon and lat ,default = 1
fprintf('######################################################### \n');
fprintf('########     Zhang-Feng Ma, Hohai University      ####### \n');
fprintf('########        Mi Jiang  , Hohai University      ####### \n');
fprintf('########         Triangulation&Interpolation      ####### \n');
fprintf('###################   9,July,2019    #################### \n');
[m,n] = size(xq);
X = x(:);Y = y(:)*ratio;Z = v(:);
xq = xq(:);yq = yq(:);
TRI = delaunayTriangulation(X,Y);
ID = pointLocation(TRI,xq,yq.*ratio);
Vx0 = X(TRI(ID,1));Vx1 = X(TRI(ID,2));Vx2 = X(TRI(ID,3));
Vy0 = Y(TRI(ID,1))./ratio;Vy1 = Y(TRI(ID,2))./ratio;Vy2 = Y(TRI(ID,3))./ratio;
Vz0 = Z(TRI(ID,1));Vz1 = Z(TRI(ID,2));Vz2 = Z(TRI(ID,3));
xkj = Vx1 - Vx0;
ykj = Vy1 - Vy0;
xlj = Vx2 - Vx0;
ylj = Vy2 - Vy0;
f = 1 ./ (xkj.*ylj - ykj.*xlj);
zj = Vz0;zk = Vz1;zl = Vz2;
zkj = zk - zj;
zlj = zl - zj;
a = -f.*(ykj.*zlj - zkj.*ylj);
b = -f.*(zkj.*xlj - xkj.*zlj);
c = -a.*Vx1 - b.*Vy1 + zk;
zout = a.* xq + b.* yq + c;
vq = reshape(zout,[m,n]);
% [ti,bc] = pointLocation(TRI,xq(:),yq(:));
% bc(isnan(ti),:)=[];ti(isnan(ti))=[];
% triVals = v(TRI(ti,:));
% vq = dot(bc',triVals')';
% vq = reshape(vq,size(xq));
end

