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=. :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=. :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 )