NB. Iterate to solution for s(tau) & p(tau),
NB. starting with s(tau)=B(tau) and p=0 
NB. [Harrington, 1970, Ap. & Sp. Sci., v. 8,
NB.  p 227, equations (7)]

s_and_p0=: 4 : 0
'L0 M0 N0'=. x
'lam B'=. y
NB. the symbol ^:(n) means "iterate n times"
NB. if n is _ (infinity), iteration is continued
NB. until convergence (to machine accuracy)
NB. thus the statement below iterates sp12
(lam;B;L0;M0;N0)&sp12^:_ B;(0*B)
)

sp12=: 4 : '(x&sp1 y);(x&sp2 y)'

sp1=: 4 : 0
NB. Evaluate rhs of 1st integral equation
'lam B L0 M0 N0'=. x
's p'=. y
(lam*B)+(1-lam)*((L0 mm s)+3%~(M0 mm p))
)

sp2=: 4 : 0
NB. Evaluate rhs of 2nd integral equation
'lam B L0 M0 N0'=. x
's p'=. y
(3%8)*(1-lam)*((M0 mm s)+(N0 mm 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)  
)