#       Given arrays X1A(length M) and X2A(length N) of independent variables
#       and an MxN array of function values YA, tabulated at the grid points
#       defined by X1A and X2A and given values X1 and X2 of the independent
#       variables, this routine returns an interpolated function value Y,
#       and an accuracy indication DY(based only on the interpolation in the
#       X1 direction, however).  From Numerical Recipes, p97, by Press et al.


procedure polin2 (x1a, x2a, ya, m, n, x1, x2, y, dy, status)

int m, n
real x1a[m], x2a[n], ya[m, n]
real x1, x2, y, dy
real yntmp[21], ymtmp[21]
int j, k
int status

begin

        status = OK

        do j = 1, m {
           do k = 1,n {
              yntmp[k] = ya[j, k]
           }
           call polint (x2a, yntmp, n, x2, ymtmp[j], dy, status)
        }

        call polint (x1a, ymtmp, m, x1, y, dy, status)

end