load'/your_path/tau_grid.ijs'
load'/your_path/MTsi.ijs'
load'/your_path/IQMsi.ijs'
load'/your_path/IQs.ijs'
load'plot'

NB. Solve problem of scattering semi-infinite grey
NB. atmosphere with polarization, constant lambda.   
NB. Emergent flux is set to unity.
NB. Usage: (tau_grid parameters) si_pol lambda 
Grey_si=: 4 : 0
lambda=. y
tau=: tau_grid x
FM0=. 0{ FM=. 2 MTsi tau  
NB. (the first row of FM gives the surface flux)
FM40=. 0{ FM4=. 4 MTsi tau
L0=. 1 MTsi tau
M0=. 3 MTsi tau
N0=. 5 MTsi tau
NB. solve for s & p:
's p'=. (L0;M0;N0;FM0;FM40) Grey_s_p lambda 
Flx=: (FM mm s)+(FM4 mm p)  NB. flux array
mu=: int01 20
IQM=. IQMsi tau;mu
'I Q pol'=. (IQM;mu) IQs s;p
NB. plot I(mu) and -10 times the polarization as function of mu
plot mu;>I;_10*pol
s;p;I;Q;pol
)

NB. Direct solution for s(tau) & p(tau),
NB. Usage: (Transform matrices) S_and_P lambda 
Grey_s_p=: 4 : 0
'L0 M0 N0 FM0 FM40'=. x NB. the transform arrays
scat=. (3%8)*1-y        NB. (3/8)*(1-lambda)   
Idn=. =@i. n=. #L0      NB. Idn is n x n identity matrix
a=. FM0+"(1) (L0- Idn)  NB. add surface flux eqn to each row 
b=. FM40+"(1) M0%3 
c=. scat*M0
d=. - Idn - scat*N0
Q=. (a,"1 b),"2 c,"1 d  NB. form grand coefficient matrix Q=[a,b/c,d]
rhs=. (n#1),(n#0)       NB. right hand side, assuming F = 1
sp=. rhs %. Q           NB. solve system: rhs %. (coeff matrix)
's p'=. (2,n)$ sp       NB. 1st n elements are "s"; 2nd n are "p"
)

NB. For example, case of lambda=0.2 
NB. Run "Grey_si":
NB. z=: (0.01 8 20) Grey_si 0.2

NB. after running this script, the variables are stored in "z"
NB. unpack them with:  's p I Q pol'=. z
NB. then list them with:  |: >tau;s;p
NB. and with:  |: >mu;I;Q;pol
NB. check flux constancy with
NB. plot tau;Flx