NB. Matrix for emergent intensities for slab.
load '/your_path/SpMD.ijs'

IQMs=: 3 : 0
't mu'=. y
'E0 E1 E2 E3'=. mu EE 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
Y=. -hi*"1 (a*"1 E0)-E1
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
V=. -(a*a2-h2)*"1 E0 
V=. V+ ((3*a2)-h2)*"1 E1
V=. h6*"1 V- (3*a*"1 E2)-E3
C=. SpMD 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 various mu (= cos angle)       
NB. Usage: mu EE tau --> E0;...
EE=: 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
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
E1=. x* (K1 b) -K1 a
E2=. (*:x)* (K2 b) -K2 a
E3=. (x**:x)* (K3 b) -K3 a
E0;E1;E2;E3
)