NB. Bicubic interpolation
NB. Equivalent to bcuint.for of Numerical Recipes, p. 120

NB. XA defines the cell in grid coordinates x1 & x2: the edges 
NB. of the bounding rectangle are x1L(ower), x1U(pper), x2L, 
NB. x2U. y contains the four function values in the following
NB. order: y(x1L,x2L),y(x1U,x2L),y(x1U,x2U),y(x1L,x2U).
NB. y1 contains values of dy/dx1 at these same four points.
NB. y2 contains dy/dx2 values, y12 contains d^2y/dx1 dx2.
NB. Usage: (y;y1;y2;y12;x1L x1U x2L x2U) bcuint x1 x2

bcuint=: dyad define
NB. rhs may be (n,2) array of x1,x2 pairs
'x1 x2'=. |: y
'y y1 y2 y12 XA'=. x
'x1L x1U x2L x2U'=. XA
t=. (x1-x1L)% d1=. x1U-x1L
u=. (x2-x2L)% d2=. x2U-x2L
c=. bcucof y;y1;y2;y12;d1;d2
z=. (|: c p."1 u) p."2 t
f=. [: }. 0 1 2 3"_ * ]
c1=. |: f c
z1=. d1%~ (|: c1 p."1 t) p."2 u
c2=. f &. |: c
z2=. d2%~ (|: c2 p."1 u) p."2 t
NB. result is (n,3) array of y, dy/dx1, dy/dx2 triplets: 
|: z,z1,:z2
)

NB. Coefficients for bicubic interpolation.
NB. Equivalent to bcucof.for of Numerical Recipes, p. 119
bcucof=: monad define
'y y1 y2 y12 d1 d2'=. y
4 4$ bicu_wts mm y,(d1*y1),(d2*y2),(d1*d2*y12)
)

NB. Define matrix of coefficients for bcucof:
bicu_wts=: ". ;. _2 ]  0 : 0
 1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0
_3  0  0  3  0  0  0  0 _2  0  0 _1  0  0  0  0
 2  0  0 _2  0  0  0  0  1  0  0  1  0  0  0  0
 0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0
 0  0  0  0 _3  0  0  3  0  0  0  0 _2  0  0 _1
 0  0  0  0  2  0  0 _2  0  0  0  0  1  0  0  1
_3  3  0  0 _2 _1  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0 _3  3  0  0 _2 _1  0  0
 9 _9  9 _9  6  3 _3 _6  6 _6 _3  3  4  2  1  2
_6  6 _6  6 _4 _2  2  4 _3  3  3 _3 _2 _1 _1 _2
 2 _2  0  0  1  1  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  2 _2  0  0  1  1  0  0
_6  6 _6  6 _3 _3  3  3 _4  4  2 _2 _2 _2 _1 _1
 4 _4  4 _4  2  2 _2 _2  2 _2 _2  2  1  1  1  1 
)