NB. SLAB GEOMETRY -- Lower boundry is open: I(T,-mu)=0.
NB. Usage:  idx MTs (tau grid) -> matrix operator
NB. idx MTs tau_grid -> spline matrix operator
NB. idx= 1,2,3,4,5 -> Lambda,Phi,M,Phi^4,N
load '/your_path/Ei.ijs'
load '/your_path/SpMD.ijs'

MTs=: 4 : 0
pm=. }:"(1) _1+ +: y<:/y
if. x=0 do. t=. y%mu else. t=. y end.
'I0 I1 I2 I3'=. x IIs 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)
NB. unit vectors of length rows, columns:
uc=. }. ur=. (#t)#1 
u1=. ur*/ b*(b2-h2)
u2=. -t*/((3*b2)-h2)
u3=. ((3*t2)*/b)-(t3*/uc)
U=. I0* u1+u2+u3
u1=. -ur*/ (3*b2)-h2
u2=. (6*t*/b)-(3*t2)*/uc
U=. U+ pm*I1*(u1+u2)
U=. U+ I2*3*(-t-/b)
U=. h6*"1 U- pm*I3
v1=. ur*/ a*(a2-h2)
v2=. -t*/((3*a2)-h2)
v3=. (3*t2*/a)-(t3*/uc)
V=. -I0* v1+v2+v3 
v1=. ur*/ (3*a2)-h2
v2=. ((3*t2)*/uc)-6*t*/a
V=. V+ pm*I1*(v1+v2)
V=. V+ I2*3*(t-/a)
V=. h6*"1 V+ pm*I3
C=. SpMD t
spl=. (U mm }:C)+(V mm }.C)
lin+spl
)

NB. make table of I-integrals for MTs 
IIs=: 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
I0;I1;I2;I3
)