NB. Monte Carlo polarized scattering of beam by a slab, using NB. precomputed (e.g., Mie) scattering matrix. Usage: NB. (half_tau;mu0;hnmu;nang) PB_Mie (# photons),(# scatterings),(# scat kept) PB_Mie=: dyad define '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 trimmed frets Sorted=: 0$~(>(1+#)&.>bfr),4 NB. place holder for sorted results SS=: (1,$Sorted)$ ,Sorted NB. SortSum framework, new axis to acumulate SS 'n nscat nkeep'=: 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 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. 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 Theta & phi for random scattering MonCar=: monad define 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) 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' 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 *linear* polarization from (I,Q,U,V) vector or NB. (4,n) array. added 1<. so result can't exceed unity. POL=: 1 <. {. %"1~ [: L 1 2 { ] 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 NB. scattering, 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 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 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=: 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 01{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 SSu;SSd;SSNu;SSNd;remain zres writeB PB_dir,file NB. overwrite old file with latest accumulated results ((": count),LF) (1!:2) < Cyc_Log count=. >: count end. zres ) Setup=: monad define dirct=. '/your_path/' NB. directory with scattering matrix print 'Enter name of file with tabulated scattering matrix: ' Smat=. enter'' file=. dirct,Smat NB. full path to scattering matrix 'angdeg S11 S12 S33 S34'=: |: mread file NB. read in scattering matrix table 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 NB. (but if not, normalize.) 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 print 'S12b = '; S12b;' .... OK ' ) 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) lintrpY 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 ) PBM_dir=: '/your_path/' Cyc_Log=: PBM_dir,'PBM_CYC.log' load'/your_path/splines.ijs' NB. load spline routines 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 1st 10; repeat 100 times: NB. Z=. PBM_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)