NB. Matrix for emergent intensities from slab with
NB. source symeteric about the mid-plane.
load '/your_path/SpMDD.ijs'

IQMss=: 3 : 0
't mu'=. y
ts3=. ts* ts2=. *: ts=. +: {: t
'E0 Ef0 E1 Ef1 E2 Ef2 E3 Ef3'=. mu EEss t
t3=. t* t2=. t*t
a2=. *: a=. }: t
b2=. *: b=. }. t
h2=. *: h=. b-a
h6=. 6%~ hi=. % h
X=. hi*"1 (b*"1 E0)-E1
bs=. b-ts
X=. X+ hi*"1 (bs*"1 Ef0)+Ef1
Y=. -hi*"1 (a*"1 E0)-E1
as=. a-ts
Y=. Y- hi*"1 (as*"1 Ef0)+Ef1
lin=. (X,"(1) 0)+(0,"(1) Y)  NB. linear part
U=. (b*b2-h2)*"1 E0
U=. U- ((3*b2)-h2)*"1 E1
U=. h6*"1 U+ (3*b*"1 E2)-E3
u1=. (b*b2-h2)-ts*(3*b2)-h2
Uf=. (u1+(3*ts2*b)-ts3)*"1 Ef0
u1=. ((3*b2)-h2)+(_6*ts*b)+3*ts2
Uf=. Uf+ u1*"1 Ef1
U=. U+ h6*"1 Uf+ ((3*(b-ts))*"1 Ef2)+ Ef3
V=. -(a*a2-h2)*"1 E0 
V=. V+ ((3*a2)-h2)*"1 E1
V=. h6*"1 V- (3*a*"1 E2)-E3
v1=. (a*a2-h2)-ts*(3*a2)-h2
Vf=. -(v1+(3*ts2*a)-ts3)*"1 Ef0 
v1=. ((3*a2)-h2)+(_6*ts*a)+3*ts2
Vf=. Vf- v1*"1 Ef1
V=. V+ h6*"1 Vf+ (_3*(a-ts)*"1 Ef2)- Ef3
C=. >0{ ('f';0) SpMDD t
spl=. (U mm }:C)+(V mm }.C) NB. spline correction
lin+spl
)

NB. table of E-integrals (= mu^n int x^n exp(-x) dx)  
NB. emergent radiation at cos angles mu      
NB. Usage: mu EEsi tau --> E0;...
EEss=: 4 : 0
a=. }:"1 x%~/y    NB. table of tau_i/mu_j 
b=. }."1 x%~/y    NB. table of tau_(i+1)/mu_j
tf=. (+: {: y)-y  NB. far side: 2*tau(N)-tau(i)
d=. }:"1 x%~/tf   NB. 2*tau(N)-tau(k)
c=. }."1 x%~/tf   NB. 2*tau(N)-tau(k+1)
K0=. 3 : '-(^-y)'
K1=. 3 : '-(^-y)*(1+y)'
K2=. 3 : '-(^-y)*(2+(2*y)+y*y)'
K3=. 3 : '-(^-y)*(6+(6*y)+(3*y*y)+y^3)'
E0=. (K0 b) -K0 a    NB. near side contribution
Ef0=. (K0 d) -K0 c   NB. far side contribution
E1=. x* (K1 b) -K1 a
Ef1=. x* (K1 d) -K1 c
E2=. (*:x)* (K2 b) -K2 a
Ef2=. (*:x)* (K2 d) -K2 c
E3=. (x**:x)* (K3 b) -K3 a
Ef3=. (x**:x)* (K3 d) -K3 c
E0;Ef0;E1;Ef1;E2;Ef2;E3;Ef3
)