```/*
DIGIT2.I
2D versions of digitize and interp functions (based on mesh_loc)

\$Id\$
*/
/*    Copyright (c) 1996.  The Regents of the University of California.

func digit2 (y0,x0, y,x,reg, pt=)
/* DOCUMENT digit2(y0,x0, y,x)
-or- digit2(y0,x0, y,x,reg)

return the index of the zone of the point or points (X0,Y0)
in the quadrilateral mesh (X,Y) with the optional region
array REG.  The result has the same dimensions as the input
X0 and Y0 arrays.  The result is <=0 at points outside the mesh.

By default, the zone index is an index into an (M-1)-by-(N-1)
array, if X and Y are M-by-N.  However, if the keyword pt= is
non-nil and non-zero, the return value is the index into an
M-by-N array in which the first row and column are non-existent
(like the optional REG array).

*/
{
zero= array(0, dimsof(x0,y0));
ndx= mesh_loc(y0+zero,x0+zero, y,x,reg);
ndx+= zero;  /* works around bug- mesh_loc can return scalar */
if (pt) return ndx;
m= dimsof(y)(2);
return ndx - m - (ndx-1)/m;
}

func interp2 (y0,x0, z,y,x,reg, outside=)
/* DOCUMENT z0= interp2(y0,x0, z,y,x)
-or- z0= interp2(y0,x0, z,y,x,reg)

return the bilinear interpolate of the function Z(X,Y) at the
points (X0,Y0).  The X, Y, and optional REG arrays specify a
quadrilateral mesh as for the plm function.  The Z values are
specified at the vertices of this mesh, so Z must have the
same dimensions as X and Y.

Points outside the mesh get the value 0.0, unless the outside
keyword is non-nil, in which case they get that value.

*/
{
scalar= !dimsof(x0,y0)(1);
if (scalar) { x0= [x0];  y0= [y0]; }
ndx= digit2(y0,x0, y,x,reg, pt=1);

/* first handle points inside mesh */
if (numberof(list)) {
ndx= ndx(list);
zero= array(0., dimsof(x0,y0));
x0= (x0+zero)(list);
y0= (y0+zero)(list);
/* here are corners of the interpolation tets */
m= dimsof(y)(2);
x00= [x(ndx-m-1), y(ndx-m-1), z(ndx-m-1)];
x10= [x(ndx-m), y(ndx-m), z(ndx-m)];
x11= [x(ndx), y(ndx), z(ndx)];
x01= [x(ndx-1), y(ndx-1), z(ndx-1)];
/* form the centroid, make centroid the origin */
xc= 0.25*(x00+x10+x01+x11);
x0= xc(,1)-x0;
y0= xc(,2)-y0;
/* form the median vectors */
xm= 0.5*[x11+x10-x01-x00, x11+x01-x10-x00, x11+x00-x10-x01];
xu= xm(,1,1);  yu= xm(,2,1);  zu= xm(,3,1);
xv= xm(,1,2);  yv= xm(,2,2);  zv= xm(,3,2);
xw= xm(,1,3);  yw= xm(,2,3);  zw= xm(,3,3);
x00= x10= x11= x01= xm= [];
/* form various cross products */
cwu= xw*yu - yw*xu;
cwv= xw*yv - yw*xv;
cuv= xu*yv - yu*xv;
cup= xu*y0 - yu*x0;
cvp= xv*y0 - yv*x0;
cwp= xw*y0 - yw*x0;
/* compute the discriminant */
cuv*= 0.5;
cwu*= 2.0;
cwv*= 2.0;
bu= cwp - cuv;
bv= cwp + cuv;
tmpa= bu*bu;
tmpb= bv*bv;
use= double(tmpa<=tmpb);
dis= sqrt(max(use*(tmpa-cwu*cvp) + (1.-use)*(tmpb-cwv*cup), 0.0));
/* only one solution is the correct one, but there is no way to
know which it is without computing both (want to allow bowtied
zones and other pathologies) */
if (numberof(list)) {
tmpa= -bu(list)-dis(list);
tmpb= cwu(list);
epsa= double(!tmpa)*1.e-99;
epsb= double(!tmpb)*1.e-99;
u11= tmpa/(tmpb+epsb);
u21= cvp(list)/(tmpa+epsa);
}
if (numberof(list)) {
tmpa= -bu(list)+dis(list);
tmpb= cwu(list);
epsa= double(!tmpa)*1.e-99;
epsb= double(!tmpb)*1.e-99;
u22= tmpa/(tmpb+epsb);
u12= cvp(list)/(tmpa+epsa);
}
u11= u21= u12= u22= [];
if (numberof(list)) {
tmpa= -bv(list)-dis(list);
tmpb= cwv(list);
epsa= double(!tmpa)*1.e-99;
epsb= double(!tmpb)*1.e-99;
v21= tmpa/(tmpb+epsb);
v11= cup(list)/(tmpa+epsa);
}
if (numberof(list)) {
tmpa= -bv(list)+dis(list);
tmpb= cwv(list);
epsa= double(!tmpa)*1.e-99;
epsb= double(!tmpb)*1.e-99;
v12= tmpa/(tmpb+epsb);
v22= cup(list)/(tmpa+epsa);
}
v11= v21= v12= v22= [];
/* compute the two z values and select the proper one */
xc= xc(,3);
z1= zw*2.*u1*v1 + zv*v1 + zu*u1 + xc;
z2= zw*2.*u2*v2 + zv*v2 + zu*u2 + xc;
}

/* just punt on points outside mesh */