load'/your_path/MTs.ijs'
load'/your_path/IQMs.ijs'
load'/your_path/IQs.ijs'
load'/your_path/tau_grid2.ijs'
NB. Solve problem of scattering slab illuminated by beam.
NB. Incident intensity makes (cos) angle imu to normal.
NB. Usage: Beam (tau array; mu_in of beam; mu out array)  
Beam=: 3 : 0
'tau mu_i mu'=. y
NB. source term = point of 1st scattering of incident beam
NB. normalized by (1/mu_in) as slanted beam is spread out
NB. source (per steradian) = (1/4*pi)(1/mu_i)e^(-tau/mu_i)
star=: 4p1%~ (^-tau%mu_i)%mu_i  
F=. 2 MTs tau            NB. get Phi transform
LL=. 1 MTs tau           NB. get Lambda transform
Id=. =i. #tau            NB. identity matrix
S=: star %.(Id - LL)     NB. source function (solution)
E=. IQMs tau;mu          NB. get matrix for emergent intensity
Itop=. E mm S            NB. Scattered intensity out top
Ibot=. E mm |.S          NB. and out bottom
Flx=: F mm S             NB. flux array [(+) =up, (-) =down]
Ftot=. 1p1*({. - {:)Flx  NB. flux scattered out top + bottom (bot<0)
Fck=. 1 - ^(-{:tau)%mu_i NB. should equal fraction removed by slab
err=. (Ftot-Fck)%Fck     NB. error in emergent flux
(Ftot,Fck,err);Itop;Ibot
)

NB. Example:
NB. Get tau's for slab of total depth = 3:
NB. tau=. tau_grid2 0.001 10 1.5
NB. mu=. int01 40
NB. run routine for mu(incident)=0.864665 (30 deg)
NB. 'zz Itop Ibot'=. Beam tau;0.864665;mu
NB.  zz   NB. look at flux conservation:
NB. 0.96888 0.968869 1.11415e_5
NB. load'plot'
NB. 'pensize 2' plot tau;>S;Flx     NB. plot source fn & flux
NB. 'pensize 2' plot mu;>Itop;Ibot  NB. plot emergent intensities
NB. 1p1* ({. Flx),({: Flx)          NB. (physical) flux from top, bottom
NB. 0.651784 _0.317095 
NB. compare to Monte Carlo run:
NB. load '/your_path/MC_beam.ijs'
NB. 'mubin Su Sd'=. (3;0.864665;10) MC_beam 5e6;40
NB. 'zz Itop Ibot'=. Beam tau;0.864665;mubin  NB. run "Beam" with mu=mubin
NB. 'pensize 2' plot mubin;>Su;Sd;(1p1*Itop);(1p1*Ibot) NB. compare outputs