NB. Solve for Source function of two-level atom
NB. Usage: Two_Level tau;B(tau);epsilon(tau)
Two_Level=: monad define
'tau B eps'=. y
n=. # dtau=: DEL tau
B=. BIN B
eta=. BIN (1-eps)%eps
B%. (unit n)+ (eta%dtau)*MMs tau
)

unit=: =@i.  NB. unit n --> n x n unit matrix

NB. Construct logarithmic grid over the interval [0,2T], 
NB. symmetric about central palne at tau=T.
NB. tau1 <- first optical depth point, e.g. tau1=: 0.01
NB. npd <- number of points/decade, e.g. npd=: 5
NB. T <- Half-depth of slab, e.g. T=: 4
tau_grid2=: monad define
'tau1 npd T'=. y
t=. T,~(i. +/t<T){ t=: 0,tau1*10^(npd %~ i. 500)  NB. tau points
t, }. |. (2*T)-t
)

NB. Let's construct a grid dividing a slab of optical depth 10^5:
NB. tau=. tau_grid2 0.02 10 5e4
NB. $tau
NB. 131                  NB. 131 tau points
NB. B=. (#tau)#1         NB. B=1 throughout slab
NB. eps=. (#tau)#1e_4    NB. epsilon = 1e_4, eta=(1-eps)/eps
NB. ]S=. Two_Level tau;B;eps
NB. 0.0136976 0.0138147 0.0138643 0.0139246 0.0139978 0.0140866 0.0141941 0.0143242 ....
NB. the value at the edge should be 0.01 (sqrt(eps)) -- hmmm....
NB. (CEP method worst at surface, but good at depth.)
NB. load'plot'
NB. 'pensize 3' plot (BIN tau);S