NB. Rayleigh scattering from a point source at x=y=0, z=z_s of slab.
NB. Usage: (T; z_s; hnmu; nxy; dxy) PPaSrj k; nkeep; nscat
NB. Here, T=slab thickness, z_s= location of source along z-axis
NB. hnmu= number of (positive) mu bins, nxy= number of positive x and y 
NB. bins, total array will be (2*nxy x nxy), dxy= size of x & y bins
NB. k= number of Monte Carlo scattering photons, scatter "nscat" times
NB. we keep seperate arrays for the first nkeep scattering as well as
NB. the sum of all nscat. We need to set nkeep < nscat.
NB. The albedo is unity - other albedos can recovered from seperate scatterings.
PPaSrj=: 4 : 0
'T z_s hnmu nxy dxy'=: x         NB. T is total optical depth through slab
T2=: -: T                        NB. optical depth to midplane
nmu=: +: hnmu                    NB. total (positive and negative) mu bins
fr1=. _1+ 2*int01 nmu            NB. frets to divide up the [-1,1] range of mu
fr2=. nxy*dxy*(_1+2*int01 +:nxy) NB. divide up x-coord of emergent point
fr3=. nxy*dxy*int01 nxy          NB. y-coordinate bins: by symmetry, (y>=0) only 
bfr=: (}:}.fr1);fr2;fr3          NB. boxed frets; mu frets trimmed
Sorted=: 0$~(>: #&>bfr),3        NB. place holder for sorted results 
SS=: (1,$Sorted)$ ,Sorted        NB. SortSum framework, new axis to acumulate SS
'k nkeep nscat'=: y              NB. number of initial photons and scatterings
S=. |: (k,3)$1 0 0               NB. Initial Stokes parameters (unpolarized)
z=. k#z_s                        NB. initially emitted photons at z=z_s.           
P=. >0;0;z                       NB. starting x,y,z coordinates of all photons 
C=. DCs k                        NB. get direction cosines of random directions    
remn=: 1                         NB. intensity not yet escaped after each scatter
'S z C P n0'=. Step^:nkeep S;z;C;P;0  NB. do "nkeep" escape/scatter cycles     
SS=: k%~ }. SS                   NB. trim off empty 1st item & normalize 
Step^:(nscat-nkeep) S;z;C;P;1    NB. do remaining (nscat-nkeep) escape/scatter cycles     
Sorted=: Sorted% k               NB. normalize to total (+/+/+/ 0{Sorted)= 1.
Sorted=: Sorted+ +/ }.SS         NB. sum the (nkeep-1) scatterings & add to rest
upind=. hnmu+ i. hnmu            NB. mu indicies of upward radiation
dnind=. |. i. hnmu               NB. reversed mu indicies of downward radiation
Sup=: trim upind{ Sorted         NB. seperate upward and downward streams  
Sdn=: trim dnind{ Sorted         NB. trim empty borders & 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)* Sdn       
Su=. trim2 upind{"4 SS           NB. Stokes parameters of upward radiation
Sd=. (1 1 _1)*"1 dnind{"4 SS     NB. Stokes of downward, U fliped
Sd=. trim2 Sd                    NB. trim & reorder axes  
mubin=. upind{(%nmu)+ }: fr1     NB. the value of mu at the bin centers
xbins=. ({.fr2),(BIN fr2),{:fr2  NB. x-coordinates of bin centers + ends
ybins=. (BIN fr3),{:fr3          NB. y-coordinates of bin centers + far end
mubin;xbins;ybins;Sup;Sdn;Su;Sd
)

NB. One cycle of escape and re-scatter.
Step=: 3 : 0      
'S z C P key'=. y
dst=. (T2%|mu)-z%mu=: 2{C       NB. dst= optical depth to edge; extn= extinction
Pout=. P+ dst *"1 C             NB. coordinates of exiting photon's escape point
Rout=. L }: Pout                NB. radial distance from z-axis of escape point
phi_out=. atan2 |: }: Pout      NB. angular coordinate of escape point
phi_ray=. atan2 |: }: C         NB. angular coordinate of escape trajectory
rphi=. 2p1| phi_out - phi_ray   NB. angle between position & trajectory
xout=. Rout* cos rphi           NB. x-coordinate of escape point
yout=. Rout* sin rphi           NB. y-coordinate of escape point
sgny=. 1- +: 0>yout             NB. sgny=-1 if yout<0 otherwise =+1
Sout=. S*"1 extn=. ^-dst        NB. flux of escaping photons         
Sout=. (sgny* 2{Sout) 2}Sout    NB. reverse U if yout < 0
Sp=. bfr Sort mu;xout;(|yout);Sout  NB. sum into bins; positive y only 
if. key=0 do.                   NB. for scattering <= nkeep
SS=: SS,Sp end.                 NB. append to earlier esapes 
if. key=1 do.                   NB. for scattering > nkeep
Sorted=: Sorted+ Sp end.        NB. add to the earlier escapes 
S=. S*"1 (1-extn)               NB. intensities of scattered fraction
remn=: remn, (k%~ +/0{S)        NB. integrated intensity still in slab (global)
Dtau=. Dt (1-extn)*rand k       NB. travel distance of scattered (fractional) photons 
z=. 2{ P=. P+ Dtau *"1 C        NB. new position of next scattering
'S C'=. C SCAT S
S;z;C;P;key
)

NB. Get angles for *random* scattering of S'=(I Q U) 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=: 4 : 0
'THETA phi'=. Thphi=. MonCar y  NB. get (random) 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 Rayleigh 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=: 3 : 0
phi=. PHI y                     NB. get random phi
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)
>(rj_TH alpha); phi
)

NB. PHI S --> random phi angles
PHI=: 3 : 0
p=. POL y                       NB. fractional polarization
chi=. CHI y                     NB. angle of polarization 
2p1| chi+ rj_PHI p              NB. solve for random phi
)

NB. Find random 0<phi'<2*pi angles for various polarizations
NB. "p" for the case of Rayleigh scattering, using the
NB. rejection method. Usage: rj-PHI p --> phi'
rj_PHI=: 3 : 0
inr=. I. bi=. (n=. #y)#1   NB. inr are indices still to do
phip=. n#0                 NB. placeholder for phi' results
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 done from inr
  n=. +/bi                 NB. number left to do
end. phip
)
NB. Is r0*(1+p/2) <= 1 - (p/2)cos(2*r) ?
p_test=: 4 : '((1+ -:x)*(rand #y))<: 1- -:x*cos +:y'

NB. Fractional polarization from (I,Q,U) vector or (3,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) vector or (3,n) array
CHI=: [: -: [: atan2r 2 1 { ]

NB. Get random _1<mu<1 & theta angles, for various polarizations
NB. "alpha", for Rayleigh scattering, using the rejection method. 
NB. Usage: rj-TH alpha --> theta
rj_TH=: 3 : 0
inr=. I. bi=. (n=. #y)#1   NB. inr are indices still to do
mus=. n#0                  NB. placeholder for mu results
while. n>0 do. 
  r=. _1+ +: rand n
  b2=: (inr{y) th_test r   NB. call th_test to test r values
  inn=. I. b2              NB. subarray index of those which pass    
  idx=. inn{inr            NB. corresponding index in whole array
  mus=. (inn{r) idx}mus    NB. insert passed mu's into mus array 
  inr=. I. bi=. 0 idx} bi  NB. remove index of passed from inr
  n=. +/bi                 NB. number left to do 
end.           
acos _1>. 1<. mus          NB. convert mu's to theta's
)
NB. Is 2*r0 <= (1+mu^2)-alpha*(1-mu^2) ?
th_test=: 4 : '(+: rand #y) <: (1+*:y)-x*(1-*:y)'

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=: 4 : 0      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) under counter-clockwise rotation 
NB. of the basis vectors (looking along propagation direction)
NB. Rotate (3,n) set of Stokes parameters by angles th(n)
NB. Usage: th(n) Rot S(3,n) --> S'(3,n)
Rot=: 4 : 0
c=. cos +:x
s=. sin +:x
q=. (c* 1{y)+ s* 2{y
u=. (c* 2{y)- s* 1{y
>(0{y);q;u
)

NB. operation which transforms (I,Q,U) under Rayleigh scattering
NB. but keeps total intensity constant
RayN=: 4 : 0  NB. theta RayN S ---> S'
c2p=. >: c2=. *:cx=. cos x
c2m=. <: c2
Ix=. 0{y
scale=. Ix% (c2p* Ix)+(c2m* 1{y)
Qsp=. (c2m* Ix)+(c2p* 1{y)
Usp=. 2*cx* 2{y
>Ix;(scale*Qsp);(scale*Usp)
)

NB. DCs n  --- gets array of n random direction cosines 
DCs=: 3 : 0
dc=. _1+ 2*rand y  NB. random mu= cos theta = z direction cosine
st=. %: 1- *:dc    NB. sin theta(s)
xc=. st*cos phi=. 2p1* rand y
yc=. st*sin phi
xc,yc,:dc
)

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=: 3 : 0
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=: 3 : 0 
theta=. acos 2{y
theta,: atan2r (1 0){y
)

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=: 4 : 0
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=: 4 : 0
s=. cos -: y
f=. sin -: y
s, (f*"1 x) 
)

NB. The conjuate of the quaternions: Qb q --> q_bar
Qb=: 3 : '(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=: 4 : '}. 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=: 4 : '(*y)* x<. |y'  NB. limits range of y to (-x) <= y <= (+x)
abslim=: ([: * ])* [<. [: |]
abslim1=: * * 1<. |
NB. atan2r=: 3 : '2p1|((0>1{y)*1p1)+ 2p1+ atan %/y=. Re y'
atan2r=: 13 : '2p1|((0>1{y)*1p1)+ 2p1+ atan %/y=. Re y'
cbrt=: 13 : '(*y)*(|y)^%3'  NB. real cube root of (possibly) negative  
map=: 3 : '(1- +:1p1<y)*1p1-|(y- 1p1)'

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

NB. Given the array A, which is a function of 3 dimensions, sort it into  
NB. bins where the bin edges are provided in the arrays fr1,fr2 & fr3, 
NB. and sum the contents of each bin. The arrays a1, a2, and a3 are the  
NB. coordinates of A that we are sorting on. The output Sp has dimensions 
NB. (1+#fr1,1+#fr2,1+#fr3). Usage: Sp=.  (fr1;fr2;fr3) Sort a1;a2;a3;A 
NB. (Thanks to Raul Miller and Roger Hui.)
Sort=: 4 : 0
b=. |:x I.&> }:y
idx=. <"1 ~. b 
w=. b +//. |: >{:y         
Shape=. 0$~(>:#&>x),3        
NB. |:"_1 |: w idx} Shape        
w idx} Shape        
)

trim=: 3 : 0  NB. fix edges & reorder for use with FracPol 
'a0 a1 a2'=. |: y      NB. transpose and assign
 a0=. }. (+/><0 1{a0) 1}a0 NB. add 1st (y<0) bin to 2nd & trim off 1st
 a1=. }. (+/><0 1{a1) 1}a1 
 a2=. }. (+/><0 1{a2) 1}a2 
a0=. |: a0         NB. join to reflection in x-axis, undo transpose
a1=. |: a1        
a2=. |: a2        
>a0;a1;a2              NB. join to make (3, mu, nx, ny) array
)

trim2=: 3 : 0  NB. fix edges & reorder for use with FracPol 
'a0 a1 a2'=. |: y      NB. transpose and assign
 a0=. }. (+/><0 1{a0) 1}a0 NB. add 1st (y<0) bin to 2nd & trim off 1st
 a1=. }. (+/><0 1{a1) 1}a1 
 a2=. }. (+/><0 1{a2) 1}a2 
a0=. |: a0         NB. join to reflection in x-axis, undo transpose
a1=. |: a1        
a2=. |: a2        
>a0;"_1 a1;"_1 a2      NB. join to make (nkeep,3,mu,nx,ny) array
)

NB. Reflect the right (y>0) side through xz plane and prepend for full image.
NB. The Stokes U must be reversed in the reflection for correct angles.
NB. Usage: full SSup --> full xy plane
full=: 3 : 0     
ar=. |."1 (1 1 _1)* y  NB. take negative of U and transpose y-axis
-: ar,"1 y             NB. join to right side; 1/2 is so sum remains unity
)

NB. make polarization magnitude & angle a complex number for use
NB. in viewmat. Usage: (polarization) cvec (angle) --> complex
NB. I.e., viewmat 10{ (1{FracPol Su) cvec (2{FracPol Su)
cvec=: 4 : 'x* (-sin a)j. cos a=. (1p1%180)*y'

NB. a_series=: 4 : '(i. >:y)^~(>:y)#x' NB. generate series 1,x,x^2,x^3,...,x^y
a_series=: ([: i. [:>:]) ^~ [ #~ [:>:] NB. Usage: albedo a_series nkeep

rand=: ?@# 0:           NB. rand n -- provides "n" random numbers on [0,1].
Dt=: [: -[:^. 1-]       NB. "Dt rand n" maps [0,1] into travel on 0<tau<inf.
int01=: i.@>:%]         NB. int01 n --> divides [0 ... 1] into n intervals.
atan2=: 13 : '2p1 | ((0>1{"1 y)*1p1)+2p1+_3 o. %/"1 y'
BIN =: -:@(}. + }:)     NB. midpoints within series

NB. Run "PPaSrj" ncycle times:
PPaSrj_CYC=: 3 : 0
'T z_s hnmu nxy dxy nphoton nkeep nscat ncycle'=. y
z=. (T;z_s;hnmu;nxy;dxy) PPaSrj nphoton;nkeep;nscat 
'mubin xbins ybins SSup SSdn SSu SSd'=. z
remain=. remn
count=. 1
while. count< ncycle do.
  z=. (T;z_s;hnmu;nxy;dxy) PPaSrj nphoton;nkeep;nscat 
  'mubin xbins ybins Sup Sdn Su Sd'=. z
  SSup =. SSup  + Sup   
  SSdn =. SSdn  + Sdn   
  SSu=. SSu + Su
  SSd=. SSd + Sd
  remain=. remain + remn
  count=. >: count
  ((": count),LF) (1!:2)<Cycle_log
end.
SSup=. count%~ SSup    NB. divide by number of cycles
SSdn=. count%~ SSdn
SSu=. count%~ SSu
SSd=. count%~ SSu
remain=. count%~ remain
ZZ=. mubin;xbins;ybins;SSup;SSdn;SSu;SSd;remain
ID=. (":T),'_',(":z_s),'_',(":hnmu),'_',(":nxy),'_',(":dxy)
ID=. ID,'_',(":nphoton),'_',(":nkeep),'_',(":nscat),'_',(":ncycle)
ZZ writeB Path,ID,'.dat'
ZZ
)
Cycle_log=: '/your_path/PPaSrj_CYC.log'
Path=: '/your_path/PPaSrj_'
writeB=: ((1) 3!:1 [) 1!:2 [: < ]  

NB. For example consider a slab of total thickness T=4 with the source
NB. located at z_s=1. Then the source is tau=1 from top and tau=3 from botton
NB. face of the slab. (at the midplane). Lets scatter 200*2e5 = 4e7 photons
NB. 65 times each, looking at a 80x40 grid of the surface with a x-range of
NB. of (-6,6)=12 in tau. Let's consider 20 angles for the emerging radiation.
NB. We will keep the first 12 scatterings, plus the sum of the 65. Thus we run:
NB.    z=. PPaSrj_CYC 4 1 20 40 0.15 2e5 12 65 200
NB. This may take days to run! Check progress by "cat /your_path/PPaS_CYC.log".
NB. When complete, results are written to a file. But to examine them, set
NB.    'mubin xbins ybins SSup SSdn Su Sd remain'=. z
NB. To look at the results, consider mu index 12, which we see corresponds
NB. to an angle with respect to the perpendicular to the slab of
NB.    (180%1p1)*acos 12{mubin
NB. 51.318 degrees. (Actually the middle of the bin 51.3 - 53.1 deg.)
NB. We can then plot the log of the intensity of the emergent radiation:
NB.    'surface' plot 10^. 1e_6+ 12{0{ SSup    NB. from the upper face
NB.    'surface' plot 10^. 1e_6+ 12{0{ SSdn    NB. from the lower face
NB. The flange around the surface is the radiation beyond the border which
NB. we force into the last bin. The unscattered radiation is not included
NB. in "SSup" or "SSdn" -- it is however the first elements 0{Su and 0{Sd.
NB. The radiation remaining after the nth scattering is given in "remain":
NB.    'bar' plot remain
NB. If we look at the intensity as a contour plot:
NB.    'contour' plot xbins;ybins; 10^. 1e_6+ 12{0{ SSdn
NB. The peak is at the same location as the emerging unscattered radiation:
NB.    'contour' plot xbins;ybins; 10^. 1e_6+ 12{0{ (0{Sd)
NB. The peak is at y=0, x=3.75. This due to the radiation that escapes
NB. directly from the source without scattering. The source is (T/2)+z_s =
NB. 2.0+1.0=3.0 units above the bottom surface, so from this angle, we see it
NB. displaced 3.0*tan acos 15{mubin = 3.747 units.
NB. Next, look at the polarization. We use "FracPol" to transform the Stokes
NB. parameters in to a triplet of intensity, fractional polarization, and the
NB. angle of polarization (in degrees):
NB.    Wu=. FracPol full SSup
NB.    Wd=. FracPol full SSdn
NB. (0{Wu) is intensity, (1{Wu) is polarization magnitude, (2{Wu) is angle.
NB. Here, we've used "full" to reflect the output across the y-axis.  
NB. Looking at the polarization, we see it's strongest along the x-axis:
NB.    'surface' plot xbins; xbins;  (12{1{Wu)
NB.    'surface' plot xbins; xbins;  (12{1{Wd)
NB. We can combine the polarization magnitude and angle into a complex number:
NB.    cpol=. (12{1{Wu) cvec (12{2{Wu)
NB. We then use the J utility "veiwmat" to look at the results;
NB.    load'viewmat'
NB.    viewmat cpol
NB. The little arrows show the direction of polarization and the colors the
NB. magnitude. The verb "scale" (below) can put a color scale on the image:
NB.    viewmat (5,0,0.5) scale cpol
NB. We see the plane of polarization is generally perpendicular to the line
NB. to the projected position of the internal source.
NB. We can also look at the radiation escaping after exactly "n" scatterings.
NB. This is n{Su. For example, radiation that has been scattered once is just
NB.    W1u=. FracPol full S1u=. 1{Su
NB.    W1d=. FracPol full S1d=. 1{Sd
NB.    viewmat (12{1{W1u) cvec (12{2{W1u)
NB.    viewmat (12{1{W1d) cvec (12{2{W1d)
NB.
NB. We can average over the whole surface easily. But note that we must average
NB. the Stokes parameters, and average over the full plane not half-plane:
NB.    WAu=. FracPol SAu=. +/"1 +/"1 full SSup
NB.    WAd=. FracPol SAd=. +/"1 +/"1 full SSdn
NB.    'pensize 3' plot mubin;>(1{WAu);(1{WAd);((#mubin)#0) 
NB. We see the curve has kinks -- look at the angles (2{WAu) to see that they
NB. are either 0 (=180) or 90, as they must be by symmetry. The 90's are negative:
NB.    sg1=. _1,19#1         NB. set the polarization negative when perpendicular
NB.    sg2=. (11#_1),(9#1)   NB. to the meridian plane (angle = 90 deg).          
NB.    'pensize 3' plot mubin;>(sg1*(1{WAu));(sg2*(1{WAd));((#mubin)#0) 
NB. Note that this averaged polarization is less that 10%. The polarization at
NB. specific points on the surface may be much higher -- ~50% in some places.
NB. Imagine the mid-plane covered with these point sources -- the *average* will
NB. remain the same. Thus the value of the point source averaged over the plane
NB. must equal the solution for the semi-infinite plane parallel case with a
NB. mid-plane source. We can verify this with the code for a plane-parallel source.

NB. put scale stripes in upper left of image
NB. Usage: (n;low;high) scale image 
NB. low=lowest level, high=highest
NB. number of levels including ends= n+1
scale=: 4 : 0
'n sl su'=. x
sc=. sl+ (su-sl)*int01 n
print 'Scale Levels: ';sc
np=. n+1
tt=. y
k=.0
while. k<np do.
ind=. (k,0);(k,1);(k,2);(k,3);(k,4)
tt=. (k{sc) ind } tt
k=. k+1
end.
tt
)