load'/your_path/EIn.ijs'
load'/your_path/MTss.ijs'
load'/your_path/IQMss.ijs'
load'/your_path/tau_grid.ijs'
load'plot'
NB. Radiation from a uniform slab.
NB. For slab of thickness T, tau should span [0,T/2]
NB. (We are using the equations for a symmetric slab)
NB. Usage: Uniform (tau array;mu out array)  
Uniform=: 3 : 0
'tau mu'=. y
T=. +: {: tau           NB. total slab thickness
NB. normalize emission to give pi*(F_top+F_bot) = unity .
norm=. % T*4p1          NB. emission /unit vol/steradian
NB. The source is uniformly distributed:         
star=: (#tau)#norm 
F=. 2 MTss tau          NB. Phi-transform matrix
LL=. 1 MTss tau         NB. Lambda-transform matrix
Id=. =i. #tau           NB. identity matrix
S=: star %.(Id - LL)    NB. solve linear system --> source function 
E=. IQMss y             NB. matrix to get emergent intensity
Itop=. E mm S           NB. Scattered intensity out top
Flx=: 1p1* F mm S       NB. (physical) flux array [(+)=up, (-)=down]
Ftot=. ({. - {:) Flx    NB. flux scattered out top + bottom (bot<0)
Ftot;Itop               
)

NB. Example:
NB. Get tau's for slab of total depth = 3:
NB.     tau=. tau_grid 0.001 10 1.5
NB.     mu=. int01 40
NB.     'Ftot Itop '=. Uniform tau;mu 
NB.     Ftot   NB. look at total flux:
NB.  1.00001 
NB.     'pensize 2' plot tau;>S;Flx     NB. plot source fn & flux
NB.     'pensize 2' plot mu;Itop        NB. plot emergent intensity
NB.     ({. Flx),({: Flx)               NB. flux from top, bottom
NB. 0.500004 _0.500004