NB. Direct solution for s(tau) & p(tau),
NB. given B(tau) and  lambda(tau)
NB. [Harrington, 1970, Ap. & Sp. Sci., v. 8,
NB.  p 227, equations (7)]

NB. (You need to choose an optical depth grid "t" first, and
NB.  and then run "1 MTsi t", etc. to create L0, M0 & N0.)
NB. Usage: 's p'=. (L0;M0;N0) s_and_p lam;B

s_and_p=: 4 : 0
'L0 M0 N0'=. x         NB. the transform arrays
'lam B'=. y            NB. Plank term B(tau) and lambda(tau) 
n=. #L0                NB. n is the number of depth points
idn=. =@i. n           NB. idn is n x n identity matrix
a=. idn - (1-lam)*L0
b=. _3%~(1-lam)*M0 
c=. (_3%8)*(1-lam)*M0
d=. idn - (3%8)*(1-lam)*N0
Q=. (a,"1 b),"2 c,"1 d NB. form grand coefficient matrix  Q=[a,b/c,d]
rhs=. (lam*(n$B)),(n#0)
sp=. rhs %. Q          NB. solve linear system: rhs %. (coeff matrix)
's p'=. (2,n)$ sp      NB. 1st n elements are "s"; 2nd n are "p"
)

NB. Create array of tau points. Usage:
NB. t = tau_grid tau1 npd T
NB. tau1 <-- first optical depth point
NB. npd <--  number of points/decade
NB. T <--    maximum depth of tau grid
tau_grid=: 3 : 0
'tau1 npd T'=: y
t=: T,~(i. +/t<T){ t=: 0,tau1*10^(npd %~ i. 500)  
)

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=: 3 : 0
'tau1 npd T'=: y
t=. T,~(i. +/t<T){ t=: 0,tau1*10^(npd %~ i. 500)  NB. tau points
nt=: +: #t       NB. number of tau points (not exceeding 999)
t, }. |. (2*T)-t
)