NB. 2-D bicubic interpolation in table.
NB. no search -- computes table indices, assuming that
NB. the points are logarithmic and evenly spaced.
NB. also evaluates derivatives at interpolation point
NB. lhs is, e.g.,   z=. readb'O3_loss.dat'

log10=: 10&^.
pcder=: [: }. 0 1 2 3"_ * ]  NB. Coefficients of derivatives of polynomials

bcutabs=: dyad define
'ne t loss dlosst dlossN'=. x    NB. look-up table (loss includes derivatives)
'ne0 te'=. y                     NB. interpolation point
nt=. #t
nn=. #ne
idelt=. %-/log10 (1 0){t         NB. inverse of the log step
tcon=. - idelt * log10 0{t       NB. offset of 1st point
dt=. 1|xt=. tcon+idelt*log10 te  NB. fractional index of te
nte=. 0>.(<. xt)<. (nt-2)        NB. integral index just below te
t1=. nte{t                       NB. t below te
t2=. (nte+1){t                   NB. t above te
ideln=. %-/log10 (1 0){ne
ncon=. - ideln * log10 0{ne
dn=. 1|xn=. ncon+ideln*log10 ne0
nne=. 0>.(<. xn)<. (nn-2)        NB. integral index just below ne
n1=. nne{ne                      NB. ne below ne0
n2=. (nne+1){ne                  NB. ne above ne0
v=. (te-t1)% d1=. t2 - t1
u=. (ne0-n1)% d2=. n2 - n1
eps=. 1,d1,d2,(d1*d2)
iLL=. <nne,nte            NB. index of lower left point
iLR=. <nne,(>:nte)
iUR=. <(>:nne),(>:nte)
iUL=. <(>:nne),nte
c=. (iLL{loss),(iLR{loss),(iUR{loss),iUL{loss
c=. (4,4)$ ,bicu_wts&mm ,eps* |: (4,4)$c
z=. (c p. u) p. v
c1=. pcder c
z1=. d1%~ (c1 p. u) p. v
c2=. pcder |: c
z2=. d2%~ (c2 p. v) p. u
NB. interpolate in dlosst & dlossN grids (no derivatives):
c=. (iLL{dlosst),(iLR{dlosst),(iUR{dlosst),iUL{dlosst
c=. (4,4)$ ,bicu_wts&mm ,eps* |: (4,4)$c
zdt=. (c p. u) p. v
c=. (iLL{dlossN),(iLR{dlossN),(iUR{dlossN),iUL{dlossN
c=. (4,4)$ ,bicu_wts&mm ,eps* |: (4,4)$c
zdN=. (c p. u) p. v
z,z1,z2,zdt,:zdN
)

NB. Bicubic interpolation in loss grid without evaluation of derivatives.
bcu_simple=: dyad define
'ne t loss'=. x                  NB. "loss" grid includes derivatives
'ne0 te'=. y                     NB. interpolation point
nt=. #t
nn=. #ne
idelt=. %-/log10 (1 0){t         NB. inverse of the log step
tcon=. - idelt * log10 0{t       NB. offset of 1st point
dt=. 1|xt=. tcon+idelt*log10 te  NB. fractional index of te
nte=. 0>.(<. xt)<. (nt-2)        NB. integral index just below te
t1=. nte{t                       NB. t below te
t2=. (nte+1){t                   NB. t above te
ideln=. %-/log10 (1 0){ne
ncon=. - ideln * log10 0{ne
dn=. 1|xn=. ncon+ideln*log10 ne0
nne=. 0>.(<. xn)<. (nn-2)        NB. integral index just below ne
n1=. nne{ne                      NB. ne below ne0
n2=. (nne+1){ne                  NB. ne above ne0
v=. (te-t1)% d1=. t2 - t1
u=. (ne0-n1)% d2=. n2 - n1
eps=. 1,d1,d2,(d1*d2)
iLL=. <nne,nte                   NB. index of lower left point
iLR=. <nne,(>:nte)
iUR=. <(>:nne),(>:nte)
iUL=. <(>:nne),nte
c=. (iLL{loss),(iLR{loss),(iUR{loss),iUL{loss
c=. (4,4)$ ,bicu_wts&mm ,eps* |: (4,4)$c
(c p. u) p. v
)

NB. define bicubic coefficient matrix:
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 
)