NB. Monte Carlo polarized scattering in a slab, using the NB. Rayleigh scattering polarization matrix. No circular NB. polarization, so we do not include 4th Stokes parameter. NB. the Stokes arrays, etc., are (3 n) not (n 3) NB. The source is assumed to be uniformly distributed in the slab. NB. Usage: (half_tau; hnmu) slabU (number of photons), (number of times scattered) slabU=: dyad define 'n nscat'=: y 'T0 hnmu'=: x NB. T0 = optical depth to midplane nmu=. +: hnmu NB. total number of mu bins (+ and -) frets=: _1+ 2*int01 nmu NB. frets in mu over [-1,1] S=. |: (n,3)$1 0 0 NB. Initial Stokes parameters (unpolarized) z=. (-T0)+2*T0*rand n NB. photon source is uniform along z-axis C=. DCs n NB. get n random initial direction cosigns SS=: 1 3$SortSum=: |:(nmu,3)$0 NB. SortSum framework, new axis for SS to acumulate Round^:(nscat) S;z;C NB. do "nscat" escape/scatter cycles SS=: n%~ }. SS NB. trim off empty 1st item & normalize mubin=. (hnmu+ i. hnmu){ (%nmu)+ }: frets NB. mu values at bin centers Su=. n%~ (hnmu+ i. hnmu){"1 SortSum Sd=. n%~ |."1 (i. hnmu){"1 SortSum SNu=. (hnmu+ i. hnmu){"1 SS SNd=. |."1 (i. hnmu){"1 SS mubin;Su;Sd;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 extn=. ^-dst NB. extn= extinction to edge of slab Sout=. S*"1 extn NB. photons (S) exiting at mu Sortn=. frets Bins mu;Sout NB. sort into bins by mu and sum SortSum=: SortSum + Sortn NB. add nth scattering SS=: SS, Sortn NB. appends nth scattering to leading axis S=. S *"1 (1-extn) NB. fractional intensities if scattered again remain=: n%~ +/ 0{ S NB. fraction of original intensity still trapped 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) 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) >(TH alpha); phi ) NB. PHI S --> random phi angles PHI=: monad define p=. POL y NB. fractional polarization chi=. CHI y NB. angle of polarization 2p1| chi+ -: p Kepler 4p1*rand #p NB. solve for random phi ) 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. Solve Kepler's equation M = E - (p/2)sin(E) by N-R NB. iteration. Usage: p Kepler M --> E Kepler=: 4 : '(x;y)&Newton^:_ y' NB. iterate to convergence Newton=: dyad define hp=. -: > {. x NB. 1/2 polarization p M=. > {: x NB. "mean anomaly" y + ((M-y)+ hp*sin y)%(1- hp*cos y) NB. N-R iteration ) NB. For Rayleigh scattering: explicit solution of NB. the cubic equation (alpha=-1 is a special case) NB. TH alpha --> random Theta TH=: monad define am=. 1- y Q=: I. 0= ap=. >: y NB. Q= indices of alpha=-1 cases z2=. *: z=.(2-y)* r=. 1- +: rand #y sqr=. %: z2+ ap%~am^3 sqr=. 1 Q} sqr NB. block potential NaN in next line AB=. (cbrt ap)%~(cbrt z+sqr)+cbrt z-sqr AB=. (Q{r) Q} AB NB. solution is r if (alpha+1)=0 acos abslim1 AB ) 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) 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=: 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 ) NB. operation which transforms (I,Q,U) under Rayleigh scattering NB. but keeps total intensity constant RayN=: dyad define 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=: monad define 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 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 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=: 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 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: count end. SSu=. count%~ SSu SSd=. count%~ SSd SSNu=. count%~ SSNu SSNd=. count%~ SSNd zres=. mubin;SSu;SSd;SSNu;SSNd zres writeB Results zres ) Cyc_Log=: '/your_path/UCYC.log' Results=: '/your_path/UResults.dat' NB. Example: Source in slab of half-thickness 1.5. NB. scatter 3e5 photons 35 times; repeat 20 times: NB. 'mubin SSu SSd SSNu SSNd'=. UCYC 1.5 20 3e5 35 20 NB. get polarization fraction & angle from Stokes parameters: NB. Wu=. FracPol SSu NB. Wd=. FracPol SSd NB. plot fraction polarization vs. mu: NB. 'pensize 2' plot mubin;>(1{Wu);(1{Wd) NB. plot polarization angle (should be 0(=180) or 90 deg): NB. 'pensize 2' plot mubin;>(2{Wu);(2{Wd) NB. plot the Stokes parameters for escaping radiation after each scatter: NB. 'pensize 2' plot 0{"2 SSNd+SSNu NB. I(mu) (sum upwards+downwards) NB. 'pensize 2' plot 1{"2 SSNd+SSNu NB. Q(mu) NB. 'pensize 2' plot 2{"2 SSNd+SSNu NB. U(mu) NB. get polarization fraction and angle for escape after each scattering NB. WN=. FracPol"2 SSNu+SSNd NB. plot polarization fraction vs. mu: NB. 'pensize 2' plot 1{"2 WN NB. just for first 3 scatters (0th is direct escape = zero polarization) NB. 'pensize 2' plot 1{"2 (i. 4){ WN