NB. Matrix Transforms for semi-infinite atmosphere.
NB. linear extrapolation from last point to infinity.
NB. idx=1,2,3,4,5 -> Lambda,Phi,M,Phi^4,N
NB. Usage:  idx MTsi tau_grid -> matrix operator
load '/your_path/Ei.ijs'
load '/your_path/SpMD.ijs'

MTsi=: 4 : 0
pm=. }:"(1) _1+ +: t<:/t=. y
'I0 I1 I2 I3 IN0 INs'=. x IIsi t
t3=. t* t2=. t*t
a2=. *: a=. }: t
b2=. *: b=. }. t
h2=. *: h=. b-a
h6=. 6%~ hi=. % h
X=. hi*"1 ((-t-/b)*I0)-pm*I1
Y=. -hi*"1 ((-t-/a)*I0)-pm*I1
lin=. (X,"(1) 0)+(0,"(1) Y)
e2=. -({:hi)*INs
e1=. IN0+ ({:hi)*INs
eN=. |: e1(N-1)}e2(N-2)}(N,(N=.#t))$0
lin=. lin+ eN       NB. linear approximation
NB. unit vectors of length rows, columns:
uc=. }. ur=. N#1 
u1=. ur*/ b*(b2-h2)
u1=. u1- t*/((3*b2)-h2)
u1=. u1+ ((3*t2)*/b)-(t3*/uc)
U=. I0*u1
u1=. -ur*/ (3*b2)-h2
u1=. u1+ (6*t*/b)-(3*t2)*/uc
U=. U+ pm*I1*u1
U=. U+ I2*3*(-t-/b)
U=. h6*"1 U- pm*I3
v1=. ur*/ a*(a2-h2)
v1=. v1- t*/((3*a2)-h2)
v1=. v1+ (3*t2*/a)-(t3*/uc)
V=. -I0*v1 
v1=. ur*/ (3*a2)-h2
v1=. v1+ ((3*t2)*/uc)-6*t*/a
V=. V+ pm*I1*v1
V=. V+ I2*3*(t-/a)
V=. h6*"1 V+ pm*I3
C=. SpMD t
spl=. (U mm }:C)+(V mm }.C)
e3=. 6%~({:h)* INs
spl=. spl+ e3*/(N-2){C  NB. spline correction
lin+spl
)

NB. make table of I-integrals for MTsi           
IIsi=: 4 : 0
idx=. x
a=. | y -/ (}:y) NB. |tau(k)-tau(i)|
b=. | y -/ (}.y) NB. |tau(k+1)-tau(i)|
pm=. }:"(1) _1+ +: y<:/y
if. -.(idx=2)+:(idx=4) do. pm=. |pm end.
if. idx=1 do.
NB. the kernels for Lambda transform
K0=. 3 : '_2%~ 2&Ei"(0)y'
K1=. 3 : '2%~ (3&Ei"(0)y)-(^-y)'
K2=. 3 : '-(4&Ei"(0)y)+ 2%~y*(^-y)'
K3=. 3 : '(3*5&Ei"(0)y)- 2%~(^-y)*(3+y+y*y)'
end.
if. idx=2 do.
NB. the kernels for Phi-transform
K0=. 3 : '_2*3&Ei"(0)y'
K1=. 3 : '+:(2*4&Ei"(0)y)- ^-y'
K2=. 3 : '(_12*5&Ei"(0)y)+2*(^-y)*(1-y)'
K3=. 3 : '(48*6&Ei"(0)y)-2*(^-y)*(6+y*y)'
end.
if. idx=3 do.
NB. the kernels for M-transform
K0=. 3 : '(_2%~ 2&Ei"(0)y)+(3%2)*4&Ei"(0)y'
K1=. 3 : '-:(3&Ei"(0)y)+(_9*5&Ei"(0)y)+2*(^-y)'
K2=. 3 : '(-(4&Ei"(0)y))+(18*6&Ei"(0)y)-(^-y)*(3-y)'
K3=. 3 : '(3*5&Ei"(0)y)+(_90*7&Ei"(0)y)+(^-y)*(15+(_2*y)+y*y)'
end.
if. idx=4 do.
NB. the kernels for Phi^(4)
K0=. 3 : '((_2%3)*3&Ei"(0)y)+(2*5&Ei"(0)y)'
K1=. 3 : '((4%3)*4&Ei"(0)y)+(_8*6&Ei"(0)y)+(4%3)*(^-y)'
K2=. 3 : '(_4*5&Ei"(0)y)+(40*7&Ei"(0)y)-(4%3)*(^-y)*(4-y)'
K3=. 3 : '(16*6&Ei"(0)y)+(_240*8&Ei"(0)y)+(4%3)*(^-y)*(24+(_3*y)+y*y)'
end.
if. idx=5 do.
NB. the kernels for N_tau:
K0=. 3 : '((_5%3)*2&Ei"(0)y)+(4*4&Ei"(0)y)-(3*6&Ei"(0)y)'
K1=. 3 : '((5%3)*3&Ei"(0)y)+(_12*5&Ei"(0)y)+(15*7&Ei"(0)y)-(2%3)*(^-y)'
K2=. 3 : '((_10%3)*4&Ei"(0)y)+(48*6&Ei"(0)y)+(_90*8&Ei"(0)y)+(2%3)*(^-y)*(6-y)'
K3=. 3 : '(10*5&Ei"(0)y)+(_240*7&Ei"(0)y)+(630*9&Ei"(0)y)-(2%3)*(^-y)*(63+(_5*y)+y*y)'
end.
I0=. pm*(K0 b) -K0 a
I1=. pm*(K1 b) -K1 a
I2=. pm*(K2 b) -K2 a
I3=. pm*(K3 b) -K3 a
NB. integrals from tau(N) to infinity:
IN0=. -K0 tiN=. ({:y)-y  NB. tiN = tau(N)-tau(i)
INs=. (-K1 tiN)- tiN*IN0 
I0;I1;I2;I3;IN0;INs
)