load'/your_path/EIn.ijs'
load'/your_path/MTmp.ijs'
load'/your_path/IQMmp.ijs'
load'/your_path/IQs.ijs'
NB. we need points near the mid-plane:
load'/your_path/tau_grid2.ijs'  
load'plot'
NB. Usage: Mid_plane (tau array;mu out array)  
Mid_plane=: 3 : 0
'tau mu'=. y
T2=. {:tau   NB. tau runs from 0 to midplane at tau=T2
NB. source term is 1st scattering from mid-plane source               
NB. which turns out to be 1st exponential integral               
tt=. (}:tau),(T2- -:1{tau)  NB. to avoid Ei(0)=infty
star=: 8p1%~ 1&EIn T2-tt    NB. S*= E_1(T2-tau)/8*pi
F=. 2 MTmp tau
LL=. 1 MTmp tau
Id=. =i. #tau           NB. identity matrix
S=: star %.(Id - LL)    NB. source function (solution)
E=. IQMmp y             NB. matrix for emergent intensity
Itop=. E mm S           NB. Scattered intensity out top
Ins=. (^-T2%mu)%4p1*mu  NB. Intensity never scattered
Flx=: F mm S            NB. flux array 
Fs=. +: 1p1* {.Flx      NB. flux scattered out top + bottom 
Fns=. 2&EIn T2          NB. flux escaping slab without being scattered
err=: 1-(Fs+Fns)        NB. error = 1-(Ftot+Fns)
Itop;Ins;(Itop+Ins);Fs;Fns    
)

NB. Example:
NB. Get tau's for slab of total depth = 1.6:
NB. tau=. tau_grid2 0.0005 12 0.4 NB. T2=0.8
NB. mu=. int01 40
NB. 'Is Ins Itot Fs Fns'=. Mid_plane tau;mu    
NB.  NB. look at flux conservation:
NB.  Fs,Fns,err --> 0.798829 0.200852 0.000319381                 
NB. 'pensize 2' plot tau;>S;Flx     NB. plot source fn & flux
NB. 'pensize 2' plot mu;>Is;Ins;Itot NB. plot intensities vs mu 
NB. Compare with Monte Carlo calculation:
NB. 'mubin Su Sd'=. (1.6;10) Mid_slab 4e6; 35
NB. 'Is Ins Itot Fs Fns'=. Mid_plane tau;mubin    
NB. Sud=. -:Su+Sd
NB. 'pensize 2' plot mubin;>(1p1*Is);(1p1*Ins);(1p1*Itot);Sud