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=: 4 : 0
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=: 3 : 0
'y y1 y2 y12 d1 d2'=. y
4 4$ bicu_wts mm y,(d1*y1),(d2*y2),(d1*d2*y12)
)

require'files'
mread=: (0&".)@('m'&fread)
bicu_wts=: mread 'Bicubic_wts.dat'
	
NB. Define matrix of coefficients for bcucof:
NB. (Explicit values if not read in from Bicubic_wts.dat)
NB. bicu_wts=. _1(43 111 133 157 189 190 238 239)}256#0
NB. bicu_wts=. 1(0 24 56 59 68 92 124 127 174 196 197 220)}bicu_wts
NB. bicu_wts=. 1(221 252 253 254 255)}bicu_wts
NB. bicu_wts=. _2(40 51 108 119 132 156 181 188 191 193 217)}bicu_wts
NB. bicu_wts=. _2(235 236 237 246 247 249 250)}bicu_wts
NB. bicu_wts=. 2(48 116 173 175 182 192 216 234 244 245 248 251)}bicu_wts
NB. bicu_wts=. _3(32 100 128 152 166 170 184 187 228 229)}bicu_wts
NB. bicu_wts=. 3(35 103 129 153 165 171 185 186 230 231)}bicu_wts
NB. bicu_wts=. _4(180 232 241 243)}bicu_wts
NB. bicu_wts=. 4(172 183 233 240 242)}bicu_wts
NB. bicu_wts=. _6(167 169 176 178 224 226)}bicu_wts
NB. bicu_wts=. 6(164 168 177 179 225 227)}bicu_wts
NB. bicu_wts=. _9(161 163)}bicu_wts
NB. bicu_wts=: 16 16 $ 9(160 162)}bicu_wts