NB. Isotropic scattering in a slab with a mid-plane source. NB. Usage: (T; hnmu) Uslab k; nscat NB. Here, T=slab thickness, hnmu= number of (positive) mu bins NB. k= number of Monte Carlo scattering photons, scatter "nscat" times NB. since this is a lambda-iteration, we must set nscat >> T^2 Mid_slab=: 4 : 0 'T hnmu'=. x T2=. -: T NB. optical depth to midplane nmu=. +: hnmu NB. total (positive and negative) mu bins 'k nscat'=. y NB. number of initial photons and scatterings Int=. k# 1 NB. Set Int=1 for the "k" initial photons z=. k#0 NB. initially emitted photons have z=0 fret=: _1+ 2*int01 nmu NB. frets to divide -1Su;Sd NB. plot emergent intensities NB. they should be equal -- difference indicates statistical error NB. Let's compare to the "exact" solution at the mu's of bin centers: NB. load '/your_path/Mid_plane.ijs' NB. get matrix transform routine NB. using a routine assuming symmetric slab: tau grid only goes to mid-plane: NB. tau=. tau_grid2 0.001 10 0.5 NB. define a 32-point tau grid over [0,1]. NB. 'Is Ins Itot Fs Fns'=. Mid_plane tau;mubin NB. run code NB. plot mubin;>(1p1*Itot);((1-remain)%~(-:Su+Sd)) NB. plot comparison NB. (1-remain)%~ corrects for the photons still trapped in the slab NB. Also, we multiply the matrix transform result by pi: physical vs. NB. "astrophysical flux" (i.e., per steradian) used in matrix equations.