NB. Find Matrix to give emergent intensity from source s(x)
NB. from a semi-infinite plane-parallel atmosphere,i.e., 
NB. the integral \int_0^\infty s(x) exp(-x/mu) dx/mu.
NB. Let there be n_x points in the optical depth array x
NB. Also define an array of angle cosines "mu". Let s(x_i) be
NB. the source function. Execute  E=. IQMsi x;mu
NB. Then E mm s gives the array of intensities I(mu_i). 
load '/your_path/SpMD.ijs'
 
IQMsi=: 3 : 0
't mu'=. y
'E0 E1 E2 E3 EN0 ENs'=. mu EEsi 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
X=. X - ({:hi)*ENs
Y=. -hi*"1 (a*"1 E0)-E1
Y=. Y + EN0+ ({:hi)*ENs
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)
end=: 6%~({:h)*ENs
spl=. spl+ end*/ {:}:C  NB. spline correction
lin+spl
)

NB. table of E-integrals (= mu^n int x^n exp(-x) dx)  
NB. for emergent radiation at cos angles mu      
NB. Usage: mu EEsi tau --> E0;...
EEsi=: 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
EN0=. K0 xN=. ({:y)% x      NB. tau(N)-infty
ENs=. (x* K1 xN)-({:y)*EN0
E0;E1;E2;E3;EN0;ENs
)