NB. Compute the polarized angular scattering cross-sections.
NB. Usage:  (diel_const,x) ang_scat phi;theta(s) [degrees]
NB. For this routine, diel_const should have *positive* imag. part
angle_scat=: dyad define
ci=. 0j1
cm=. %: 0{x
cnrm=. cm sphere yy=. 1{x
qext=. 0    NB. evaluate total extinction:
nci=. >: nc
  n=. 1 while. n <: nc do.
  nm=. <: n
  cim=. % ci^(>: n)
  rf=. (1+2*n)%(n*n+1)
  z=. cim*(nm{f + ci* nm{g)
  qext=. qext + rf* 1{ +. z
  n=. >: n 
end. 
qext=. qext% *: yy
t=. 0    NB. evaluate total scattering:
n=. 1 
while. n <: nc do.
  nm=. <: n
  z1=. *: | nm{f
  z2=. *: | nm{g
  t=. t + (nm{cnrm)*z1+z2
  n=. >: n 
end.
qsca=. t % *: yy
qabs=. qext - qsca
NB. Total cross-section: absorption, scattering, extinction:
Qtot=: qabs,qsca,qext      NB. Qext = Qabs + Qsca
rtd=. 180%1p1
phi=. (>0{y)%rtd
cosphi=. cos phi
sinphi=. sin phi
th=. >1{y
theta=. th % rtd
costh=. cos theta
pnmllg=. nci genlgp theta
NB. 1=parallel, 2=perpendicular, e.g., 
NB. index 12 <=> in parallel, out perpendicular
t11=. t12 =. t21 =. t22 =. s1 =. s2 =. 0
n=. 1 
while. n <: nc do.
   nm=. <: n
   cim=. % ci^(>: n)
   p1=. (n*costh*n{pnmllg)-((n+1)*(n-1){pnmllg)
   p2=. n{pnmllg
   t11=. t11 + cim*cosphi*((p2* nm{f)+(ci*p1* nm{g))* nm{cnrm
   t12=. t12 - cim*sinphi*((p1* nm{f)+(ci*p2* nm{g))* nm{cnrm
   t21=. t21 + cim*sinphi*((-p2* nm{(-f))+(ci*p1* nm{g))* nm{cnrm
   t22=. t22 + cim*cosphi*((-p1* nm{(-f))+(ci*p2* nm{g))* nm{cnrm
   s1=. s1 + cim*((p2* nm{g)+(ci*p1* nm{(-f)))* nm{cnrm
   s2=. s2 + cim*((p1* nm{g)-(ci*p2* nm{f))* nm{cnrm
   n=. >: n 
end.
snorm=. % o. *:yy
t11=. snorm* *: | t11 
t12=. snorm* *: | t12
t21=. snorm* *: | t21 
t22=. snorm* *: | t22 
tpa=. 0.25*t11+t12+t21+t22
p11=. 2*((*: |s1)+(*: |s2))%t
pl=. 2*((*: |s1)-(*: |s2))%(t*p11)
p33p11=. 4*Re(s1*(+s2))%(t*p11)
p34p11=. 4*Im(s2*(+s1))%(t*p11)
|: (10,(#th)) $ th,t11,t12,t21,t22,tpa,p11,pl,p33p11,p34p11 
)

sphere=: dyad define
ci=. 0j1
nc=: <. xc =. 2+y+4.05*y^%3
z=. x*y
nmx=. 15 + <.(xc >. |z) 
amat=. 0
n=. nmx 
while. n > nc do.
  amat=. (n%z)-(%amat+n%z)
  n=. <: n 
end.
n=. nc 
while. n > 1 do.
  amat=. amat,((n%z)-(%(_1{amat)+n%z))
  n=. <: n 
end. 
amat=. |. amat
hkl=. (>: nc) besh y
bj=. 0{ +. 0{ hkl
fl=. gl=. cnrml=. n=. 1 
while. n <: nc do.
  rf=. 2*n*n+1
  bjm=. bj
  bj=. 0{ +. n{hkl
  b=. (n%y)+ x*(n-1){amat 
  fl=. fl,((-ci^n)*rf*((-bjm)+b*bj)%((-(n-1){hkl)+b*n{hkl))
  b=. (n%y)+(((n-1){amat)%x)
  gl=. gl,((ci^n+1)*rf*((-bjm)+b*bj)%((-(n-1){hkl)+b*n{hkl))
  cnrml=. cnrml,(1+2*n)%(rf*n*n+1)
  n=. >: n 
end.
f=: }. fl
g=: }. gl
cnrm=: }. cnrml
)

NB. Generate the Hankel functions
besh=: dyad define
a=. (1 o. y)%y
by=. -(2 o. y)%y
by=. by,((by%y)-a)
n=. 1 
while. n < (x-1) do.
  z=. (1+2*n)*(_1{by)%y
  by=. by,(z - _2{by)
  n=. >: n 
end. 
nst=. x+ <. %:y +101
t3=. 0
t2=. 1e_35
i=. nst - 1 
while. i > (x- 2) do.
  t1=. ((1+2*i)*t2%y)-t3
  t3=. t2
  t2=. t1
  i=. <: i 
end.
bj=. t3,t2
i=. (x- 2) 
while. i > 0 do.
  z=. (1+2*i)*(_1{bj)%y
  bj=. bj,(z - _2{bj)
  i=. <: i 
end. 
bj=. |. bj
((a% 0{bj)*bj)+0j1*by
)

NB. Generate associated Legendre functions.
NB. Usage: (number in series) genlgp (angle in radians)
genlgp=: dyad define
m=. #y
costh=. 2 o. y 
pnmllg =. > (m#0) ; (m#1)
n=. 2 
while. n < x do.
  a=. (_1+2*n)*costh*_1{pnmllg
  pnmllg =. pnmllg,((%n-1)*a - n*_2{pnmllg) 
  n=. >: n 
end. 
pnmllg
)

NB. Example: refractive index 2.2+0.03i, x=(wavelength/2*pi*grain_radius)= 1.3
NB. Angle between polarization plane and scattering plane, phi = 60 degrees
NB. evaluate cross-section for scattering angles of 0,30,60,90,120,150 & 180 deg.

NB.   (2.2j0.03, 1.3) angle_scat 60; 0 30 60 90 120 150 180
NB.     0    0.027835 0.0835049   0.0835049   0.027835    0.05567  3.02459        0         1           0
NB.    30   0.0199242 0.0751056   0.0597727  0.0250352  0.0449594  2.44267  0.11368  0.993516 _0.00190776
NB.    60  0.00637654 0.0555536   0.0191296  0.0185179  0.0248944  1.35253 0.487713  0.872927  _0.0116204
NB.    90 0.000220796 0.0356158 0.000662388  0.0118719  0.0120927 0.657006 0.963483  0.265357  _0.0358646
NB.   120  0.00105114 0.0218078  0.00315341 0.00726925 0.00832039 0.452052 0.747335 _0.663088  _0.0424848
NB.   150  0.00329705 0.0147074  0.00989114 0.00490246 0.00819951 0.445484 0.195794  _0.98053  _0.0150345
NB.   180  0.00421011 0.0126303   0.0126303 0.00421011 0.00842023 0.457476        0        _1           0

NB. the total (integrated over all angles) cross-section is stored in Qtot:
NB.    Qtot
NB. 0.0427298 0.462589 0.505318