NB. J code to evaluate approximate polarization
NB. of the emergent radiation. 
NB. Usage: 'I Q P' =.  (a,b,Sig0,s) Ap_Pol mu 
mu=: int01 100  NB. 100 values of mu=cos(theta)
Ap_Pol=: dyad define
'a b Sig0 s'=. x
mu=. y
I=. a + b*mu
Q=. x QQ mu
P=. _100*Q%I  NB. negative of the percent polarization
>I;Q;P
)
QQ=: dyad define   NB. evaluate equation (31).
'a b Sig0 s'=. x
m2=. *: mu=. y
beta2=. *: beta=. mu%(1+s*mu)
llg=. (1-3*beta2)*lg=. ^. 1 + %beta
head=. 0.1875*(1-m2)*Sig0%(1+s*mu)
A=. (_1.5*beta)+(3*beta2)+beta*llg
B=. _0.25+(_1.5*beta2)+(3*beta*beta2)+beta2*llg
Q=. head*((a*A)+(b*B))
)
NB. Example: On a machine running J, Start J and enter:
NB. load'Anal_Pol.ijs'
NB. 'I Q P' =. (1,2,0.1,0) Ap_Pol mu
NB. load 'plot'
NB. plot mu;P
NB. 'I1 Q1 P1' =. (1,2,0.1,1) Ap_Pol mu
NB. 'pensize 3' plot mu;>P;P1