NB. Monte Carlo polarized scattering in a slab, using the NB. Rayleigh scattering polarization matrix, from a source NB. which is located "z_s" units from the midplane z=0. NB. No circular polarization, so we do not include the 4th NB. Stokes parameter; the Stokes arrays are (3 n) not (n 3) NB. Usage: (T; z_s; hnmu) slabP (number of photons), (times scattered) NB. Here, T=slab thickness, z_s= location of source along z-axis slabP=: 4 : 0 'n nscat'=: y 'T z_s hnmu'=: x NB. T = total optical depth of slab T0=: -: T NB. 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=. n#z_s NB. photons launched from z = z_s 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 =: 3 : 0 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=: 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 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) >(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+ -: 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=: 4 : 0 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=: 3 : 0 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=: 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 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=: 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 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;remain ID=. (":T),'_',(":z_s),'_',(":hnmu),'_' ID=. ID,(":nphoton),'_',(":nscat),'_',(":ncycle) zres writeB Path,ID,'.dat' zres ) Cyc_Log=: '/your_path/CYC_count.log' Path=: '/your_path/slabP_' writeB=: ((1) 3!:1 [) 1!:2 [: < ] NB. see local_defs.ijs NB. Example: Source in slab of half-thickness 1.5. NB. Source at midplane (z=0) NB. scatter 3e5 photons 35 times; repeat 20 times: NB. 'mubin SSu SSd SSNu SSNd rem'=. CYC_P 3 0 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