NB. Monte Carlo polarized scattering of beam by a slab, for particles
NB. defined by a read-in Mie scattering polarization matrix. 
NB. This version uses the rejection method to pick phi angles and an 
NB. inverse interpolation to pick the theta(alpha) angles.
NB. (Should be OK even for strongly peaked phase functions.)
NB. Usage:(half_tau;mu0;hnmu;nang) PB_mixed (# photon),(scat. kept),(# scatterings)
PB_mixed=: dyad define
file=. '/your_path/P-Matrix.dat'           NB. full path to scattering matrix 
'angdeg S11 S12 S33 S34'=: |: mread file   NB. read in scattering matrix table
load'/your_path/splines.ijs'               NB. load spline routines
angrad=: (1p1%180)*angdeg                  NB. convert angles to radians
S11b=: -: angrad sinteg (S11* sin angrad)  NB. should be normalized to unity
'S11 S12 S33 S34'=: S11b%~ S11,S12,S33,:S34 
S12b=: -: angrad sinteg (S12* sin angrad)  NB. integral of S12(mu) 
S11D=: angrad spline S11   NB. get 2nd derivatives for spline interpolation
S12D=: angrad spline S12
S33D=: angrad spline S33
S34D=: angrad spline S34
F11=: -:angrad spintg S11* sin angrad      NB. Cumulative S11 element
S12b=: {: F12=: -:angrad spintg S12* sin angrad NB. Cumulative S12
aa=:  _1++:int01 50           NB. Grid of alpha values (_1< alpha <1)
eqn51=: dyad : '(1+y*S12b)%~ ((0{x)+y*(1{x))'   NB. equation 51 solved for r 
TAB51=: (>F11;F12)&eqn51"0 aa NB. Table of r as function of alpha & theta

'T0 mu0 hnmu nang'=: x        NB. T0= midplane optical depth; mu0= beam mu   
nmu=. +: hnmu                 NB. total number of mu bins (+ and -)
fr1=: _1+ 2*int01 nmu         NB. frets in mu over [-1,1]
fr2=: 1p1*int01 nang          NB. frets in angle over [0,pi]
bfr=: (}:}.fr1);(}:}.fr2)     NB. boxed frets; mu frets trimmed
Sorted=: 0$~(>(1+#)&.>bfr),4  NB. place holder for sorted results
SS=: (1,$Sorted)$ ,Sorted     NB. SortSum framework, new axis to acumulate SS
'n nkeep nscat'=: y           NB. number of initial photons and scatterings
nsc=: >:nscat                 NB. Oth Round not scattered, just passes through
S=. |: (n,4)$1 0 0 0          NB. Initial Stokes parameters (unpolarized)
z=. n#(T0-1e_8)               NB. photons enter at top of slab          
mux=. %: 1- *:mu0             NB. direction cosine along x-axis
C=. |: (n,3)$ - mux,0,mu0     NB. direction cosigns of entering beam (-) downward!
remn=: 1                      NB. intensity not yet escaped after each scatter
'S z C'=. Round^:(nsc) S;z;C  NB. do "nscat"+1 escape/scatter cycles 
SS=: n%~ }. SS                NB. trim off empty 1st item & normalize
remn=: }. remn                NB. trim off initial placeholder
Sorted=: +/SS                 NB. sum the scatterings
upind=: hnmu+ i. hnmu         NB. mu indicies of upward radiation
dnind=: |. i. hnmu            NB. reversed mu indicies of downward radiation
Sup=: |:"2 |: upind{ Sorted   NB. seperate upward and downward streams
Sdn=: |:"2 |: dnind{ Sorted   NB. reflect in x-axis
  NB. To compare downward (mu<0) with upward (mu>0) radiation, we must
  NB. reverse the direction of the z-axis, i.e. 45 deg --> 135 deg. This
  NB. is done by reversing the sign of the Stokes U parameter:
Sdn=: (1 1 _1 1)* Sdn         NB. doen the sign of V remain unchanged ?
SNu=: |:"2 |:"_1 upind{"3 SS  NB. Stokes parameters of upward radiation
SNd=: (1 1 _1 1)*"1 dnind{"3 SS NB. Stokes of downward, U fliped
SNd=: |:"2 |:"_1 SNd          NB. reorder axes
SNu=: (>:nkeep){. SNu         NB. keep details of first "nkeep" scatterings
SNd=: (>:nkeep){. SNd         NB. (others are included in sum but not saved)
mubin=: upind{(%nmu)+ }: fr1  NB. the value of mu at the bin centers
abins=: BIN fr2               NB. central angles of the phi bins
mubin;abins;Sup;Sdn;SNu;SNd
)

Round =: monad define              NB. do one cycle of escape and scattering 
'S z C'=. y
dst=. (T0%|mu)-z%mu=: 2{C          NB. dst= optical depth to edge of slab
b=. 1p1<~ phi=. atan2 |: (1 0){C   NB. phi is azimuthal angle of escape trajectory
phi=. ((-.b)*2p1- phi)+ b*phi      NB. map phi into [0,1p1] (b=0 if phi>1p1) 
Sout=. S*"1 extn=. ^-dst           NB. photons (S) exiting at mu
sgny=. _1+ +:b                     NB. sgny=-1 if phi>pi otherwise =+1
Sout=. (sgny* 2{Sout) 2}Sout       NB. reverse U if phi>pi   
Sp=. bfr Sort mu;phi;Sout          NB. sort into bins by mu and sum  
SS=: SS, Sp                        NB. appends nth scattering to leading axis
S=. S *"1 (1-extn)                 NB. fractional intensities if scattered again
remn=: remn, (n%~ +/0{S)           NB. integrated intensity still in slab (global) 
z=. z+ mu* Dt (1-extn)*rand n      NB. z = coordinates of new scat point
'S C'=. C SCAT S                   NB. scatter photons into new direction C
S;z;C
)

NB. Get angles for *random* scattering of S'=(I Q U V) array with
NB. dir. cosines dc', scatter through th and ph, and transform S'->S
NB. and dc'-->dc   Usage: dc' SCAT S' --> S;dc
SCAT=: dyad define
'THETA phi'=. Thphi=. MonCar y  NB. sample distribution of scattering angles
dc=. x Swerve Thphi NB. find corresponding direction cosines of outgoing beam
'th ph'=. outang=. DCinv dc     NB. and the angles corresponding to the dc
phi=. map phi       NB. change range of phi from 0<phi<2*pi to -pi<phi<pi
S0=.  phi Rot y     NB. rotate S' from meridian plane to scattering plane
S0=. THETA RayN S0  NB. change S0 due to Mie scattering through THETA
cTH=. cos THETA     NB. THETA's are scattering angles in scattering planes
sTH=. sin THETA     
cthp=. 2{x          NB. 3rd dir. cos. is just cosine of incoming theta'
cth=. 2{dc          NB. likewise for cosine of outgoing theta
sth=. sin th
NB. i2 is angle of scattering plane with outgoing meridian plane:
i2=. (*phi)* acos abslim1 (cthp-cTH*cth)%(sTH*sth)
NB. need (*phi) since if phi<0, we rotate in the opposite direction
S0=. i2 Rot S0      NB. rotate S0 to outgoing meridian plane
S0; dc
)

NB. MonCar S --> Theta & phi for random scattering
MonCar=: monad define
phi=. PHI y                     NB. get phi from probability dist.
q=. (1{y)%(0{y)  NB. q = Q/I
u=. (2{y)%(0{y)  NB. u = U/I
alpha=. (q*cos +: phi)+(u*sin +: phi)   NB. compute alpha(phi)'s
THETA=. intrp_TH alpha             NB. intrp_TH gets THETA's from prob. dist.
>THETA;phi             
)

NB. PHI S --> phi angles drawn from probability distribution
PHI=: monad define
p=. POL y                       NB. fractional polarization for S's
chi=. CHI y                     NB. and their polarization angles
2p1| chi+ rj_PHI p              NB. rj_PHI gets phi's from prob. dist.
)

NB. Find random 0<phi'<2*pi angles for various polarizations
NB. "p" for the case of Mie scattering, using the rejection
NB. method. Usage: rj-PHI p --> phi'
rj_PHI=: monad define
inr=. I. bi=. (n=. #y)#1   NB. inr are indices still to do
phip=. n#0                 NB. placeholder for phi' results
NB. p_count=: 0                NB. initialize counter
while. n>0 do.
  r=. 2p1* rand n          NB. trial value of phi'
  b2=: (inr{y) p_test r    NB. call the test function
  inn=. I. b2              NB. subarray index of those which pass
  idx=. inn{inr            NB. corresponding index in whole array
  phip=. (inn{r) idx}phip  NB. insert successful phi's
  inr=. I. bi=. 0 idx} bi  NB. remove those accepted from "inr"
  n=. +/bi                 NB. number left to do
NB. p_count=: >: p_count      NB. see how many cycles needed
end. phip
)
NB. Test: Is r0*(1+|p*S12b) <= 1 + p*S12b*cos(2*r) ?
NB. Usage: "p p_test r" where r is over [0,2*pi]
NB. (Returns binary result.)
p_test=: dyad : '((1+ |x*S12b)*(rand #y))<: 1+ S12b*x*cos +:y'

NB. Get random sample of theta angles, for various polarizations
NB. "alpha", for Mie scattering, by *reverse interpolation*. 
NB. Usage: interp_TH alpha --> theta
intrp_TH=: monad define
rd=. rand #y     NB. random draws from the theta cumulative function
toobig=. (TAB51;angrad) lintrpX rd  NB. thetas for 51 alphas at each rd
(aa;toobig) lintrpY y               NB. & interpolate for the correct alpha
)

NB. Fractional polarization from (I,Q,U,V) vector or (4,n) array
NB. added 1<. to POL=: {. %"1~ [: L }.  so result can't exceed unity.
POL=: 1<. {. %"1~ [: L }.

NB. Angle of plane of polarization wrt reference plane of
NB. (I,Q,U,V) vector or (4,n) array
CHI=: [: -: [: atan2r 2 1 { ]

NB. This computes new direction cosines from scattering angles phi and THeta.
NB. Usage: dc' Swerve Theta,:phi -->dc  (where dc' is incoming and dc is out)
NB. dc' is [3 n] array and Theta,:phi a [2 n] array of angles
NB. assumes Theta=0 ==> no deflection, Theta=pi ==> backscatter
NB. phi is angle counterclockwise from z-axis plane along direction of propagation.
Swerve=: dyad define           NB. x = DC's of incoming photons
v1=. norm (1 0 2){ (_1 1 0)*x  NB. [vectors perp to dc',z planes: dc' cross (0 0 1)]
q1=. v1 Qrot 0{y    NB. quaternion to rotate by Theta in (dc',z_axis) plane, and
q2=. (-x) Qrot 1{y  NB. (-x) to rotate *counter-clock* by phi about +dc' axis
q2=. q2 QM q1       NB. multiply quaternions for sequence of both rotations
norm q2 Pr x        NB. do rotation: norm fixes case of dc'=(0 0 1)
)

NB. Verb to transform (I,Q,U,V) under counter-clockwise rotation 
NB. of the basis vectors (looking along propagation direction)
NB. Rotate (4,n) set of Stokes parameters by angles th(n)
NB. Usage: th(n) Rot S(4,n) --> S'(4,n)
Rot=: dyad define
c=. cos +:x
s=. sin +:x
q=. (c* 1{y)+ s* 2{y
u=. (c* 2{y)- s* 1{y
>(0{y);q;u;(3{y)
)

NB. Operation which transforms (I,Q,U,V) under Mie scattering
NB. but keeps total intensity I constant: 
NB. Usage:  theta RayN S ---> S'
RayN=: dyad define  
s11=. (angrad;S11;S11D) spltrp x
s12=. (angrad;S12;S12D) spltrp x
s33=. (angrad;S33;S33D) spltrp x
s34=. (angrad;S34;S34D) spltrp x
Ix=. 0{y
scale=. Ix% (s11* Ix)+(s12* 1{y)
Qsp=. (s12* Ix)+(s11* 1{y)
Usp=. (s33* 2{y)+(s34* 3{y)
Vsp=. (-s34* 2{y)+(s33* 3{y)
>Ix;(scale*Qsp);(scale*Usp);(scale*Vsp)
)

NB. Get array of direction cosines from array of (theta,phi) angles
NB. Angles in radians, 0<=theta<pi, 0=<phi<2*pi or  -pi<phi<pi
DC=: monad define
st=. sin 0{y       NB. 0{y => theta(s)
dc=. cos 0{y
xc=. st*cos 1{y    NB. 1{y => phi(s)
yc=. st*sin 1{y
xc,yc,:dc
)

NB. Convert direction cosines to (theta,phi) pairs.
NB. results in radians, 0<=theta<pi, 0=<phi<2*pi
DCinv=: monad define
theta=. acos 2{y
theta,: atan2r (1 0){y
)

NB. Given a disordered array A of results depending on 2 variables, sort
NB. it into bins where the bin edges are provided in the arrays fr1 & fr2,
NB. while summing the contents of each bin. The vectors a1 and a2 are the
NB. variables attached to the values of A by which we sort them.
NB. The output has dimensions (1+#fr1,1+#fr2).
NB. Usage:  Sp=.  (fr1;fr2) Sort a1;a2;A
Sort=: dyad define
b=. |: > x I. &.> }:y         NB. get cell indicies of each a1,a2 pair
idx=. <"1 ~. b                NB. nub of indicies in order of appearance
w=. b +//. |: >{:y            NB. use key to sum values of A in each cell
Sh=. 0$~(>(1+#)&.>x),(#|:w)   NB. make empty frame for output
w idx} Sh                     NB. fill frame with with sums
)

L=: +/ &. (*:"_) NB. Gives the lengths of vector arrays
dot=: [: +/ *    NB. dot products of arrays of vectors: [3 n] dot [3 n]
cross=: ((1: |.[)*(_1: |. ]))-((_1: |.[)*(1:|.]))  NB. cross product

NB. Multiply quaternions: q1 QM q2 --> q3
QM=: dyad define
s1=. 0{ x
v1=. }. x
s2=. 0{ y
v2=. }. y
s=. (s1*s2) - v1 dot v2
s,(s1*"1 v2)+(s2*"1 v1)+(v1 cross v2)
)

NB. Quaternion for rotation of points by theta about unit vector u
NB. (Rotates *points*; use -theta to rotate coordinate axes)
NB. In right handed coordinates, looking in the +z direction, a
NB. positive rotation (from x to y) is *clockwise* -- we must use
NB. -theta to rotate counter clockwise (or make the axis -z). 
NB. Usage: u Qrot theta --> q  then  q Pr P --> P'
NB. u may be [3 n] array and theta [n] array.
Qrot=: dyad define
s=. cos -: y
f=. sin -: y
s, (f*"1 x) 
)

NB. The conjuate of the quaternions: Qb q --> q_bar
Qb=: monad : '(1 _1 _1 _1)* y'

NB. Rotate point using quaternion q
NB. Usage: q Pr P --> P' where q is [4 n] array of quaternions,
NB. P is [3 n] array of point coordinates, and P' also [3 n]
NB. Pr=: dyad : '}. x QM (0,y) QM Qb x'
Pr=: [: }. [ QM (0"_ , ]) QM [: Qb [

norm=: ] %"1 L      NB. normalize to unit length

Dt=: [: -[:^. 1-]   NB. "Dt rand n" maps [0,1] into travel on 0<tau<inf.
ang=: acos@(1&+@(_2&*@]))  NB. ang [0,1] --> 0<ang<1p1 (including sine(theta))
NB. abslim=: dyad : '(*y)* x<. |y'  NB. limits range of y to (-x) <= y <= (+x)
abslim=: ([: * ])* [<. [: |]
abslim1=: * * 1<. |
NB. atan2r=: monad : '2p1|((0>1{y)*1p1)+ 2p1+ atan %/y=. Re y'
atan2r=: 13 : '2p1|((0>1{y)*1p1)+ 2p1+ atan %/y=. Re y'
map=: monad : '(1- +:1p1<y)*1p1-|(y- 1p1)'

NB. I, Fractional polarization & angle from (I,Q,U,V) vector (or array)
FracPol=: monad define
pol=. POL y  
chi=. CHI y
cir=. ({: % {.) y
(0{y), pol, (180*chi%1p1),: cir 
)

NB. Run PB_mixed ncycle times to build up statistics.
PBDmix_CYC=: monad define
'T0 mu0 hnmu nang nphoton nscat nkeep ncycle'=. y
zres=. (T0;mu0;hnmu;nang) PB_mixed nphoton;nkeep;nscat
'mubin abin SSu SSd SSNu SSNd'=. zres
file=. 'PBDmix_',(":T0),'_',(":mu0),'_',(":hnmu),'_',(":nang),'_',(":nphoton)
file=. file,'_',(":nscat),'_',(":nkeep),'_',(":ncycle),'.dat'
((":1),LF) (1!:2) < Cyc_Log          
(zres=. 1;zres;remn) writeB PB_dir,file
remain=. remn
count=. 2
while. count<: ncycle do.
  'mubin abin Su Sd SNu SNd'=. (T0;mu0;hnmu;nang) PB_mixed nphoton;nkeep;nscat
  SSu=. SSu+ Su           
  SSd=. SSd+ Sd           
  SSNu=. SSNu+ SNu           
  SSNd=. SSNd+ SNd  
  remain=. remain + remn
  zres=. count;mubin;abin;count%~ &.> SSu;SSd;SSNu;SSNd;remain
  zres writeB PB_dir,file
  ((": count),LF) (1!:2) < Cyc_Log          
  count=. >: count
end.
zres
)

NB. linear interpolation: (x;y) lintrpX x0 --> y(x0)
NB. x may be an array, x0 vector of different length
lintrpX=: dyad define
'xx yy'=. x
g0=. (}: yy)% d=. DEL |:xx
g1=. (}. yy)% d
u=. (g0* }.|:xx) - g1* }:|:xx
n=. y IdX~"(1 1) xx
(n{"(1 1)|:u)+ y*"1 n{"(1 1)|:g1-g0
)
IdX=: ] I.~ [: }: [: }. [

NB. linear interpolation: (x;y) lintrpX x0 --> y(x0)
NB. x0 is a vector and y a ((#x),(#x0)) array
NB. output (#x0) values are each based on different y.
lintrpY=: dyad define
'xx yy'=. x
g0=. (}: yy)% d=. DEL xx
g1=. (}. yy)% d
u=. (g0* }.xx) - g1* }:xx
n=. y IdX~"(0 1) xx
(n{"(0 1)|:u)+ y*n{"(0 1) |:g1-g0
)

PB_dir=: '/your_path/'
Cyc_Log=: PB_dir,'PBDmix_CYC.log'

NB. Example: Slab of half-thickness 1.0, beam at mu=0.4. 
NB. Sort scattered radiation into 20 mu bins and 24 phi angles.
NB. scatter 3e5 photons 35 times keep the 1st 10; repeat 100 times:
NB.   Z=. PBDmix_CYC 1 0.4 20 24 3e5 35 10 100
NB.   'count mubin abin SSu SSd SSNu SSNd remain'=. Z
NB. get polarization fraction & angle from Stokes parameters:
NB.   Wu=. FracPol SSu
NB.   Wd=. FracPol SSd
NB. Results symmetric about the x-z plane; fill in by reflection:
NB.  abin2=. abin,(|. 2p1 -abin) NB. extend to [0, 2*pi]
NB.  W2u=. Wu,"1 |."1 Wu
NB.  W2d=. Wd,"1 |."1 Wd
NB. plot fraction polarization vs. mu:
NB.  'pensize 3' plot mubin;|:(1{Wu)
NB. plot fraction polarization vs. phi (in degrees):
NB.  'pensize 3' plot (180*abin2%1p1); 1{W2u 
NB. get polarization fraction and angle for escape after each scattering
NB.  WNu=. FracPol"3 SSNu
NB.  WNd=. FracPol"3 SSNd
NB.  WN2u=. WNu,"1 |."1 WNu
NB.  WN2d=. WNd,"1 |."1 WNd
NB. plot polarization of once scattered radiation vs. angle:
NB.  'pensize 2' plot (180*abin2%1p1); 1{ 1{WN2u
NB. polarization of radiation scattered six times vs. angle:
NB.  'pensize 2' plot (180*abin2%1p1); 1{ 6{WN2u
NB. Show polarization as a surface:
NB. 'surface' plot mubin;(180*abin2%1p1); 1{W2u
NB. and radiation emerging from bottom of slab:
NB. 'surface' plot mubin;(180*abin2%1p1); 1{W2d
NB. Emergent intensity as function of phi & mu:
NB. 'surface' plot mubin;(180*abin2%1p1);(mubin%~ 0{W2u)