NB. Monte Carlo polarized scattering of beam by a slab, using two NB. precomputed (e.g., Rayleigh & Mie) scattering matrices. The scattering albedos NB. and fractional contribution of the two scattering phase functions are assumed NB. to vary with depth z. Usage: NB. (half_tau;mu0;hnmu;nang) PB_Mie3 (# photons),(# scatterings),(# scat kept) PB_Mie3=: dyad define 'T0 mu0 hnmu nang'=: x NB. T0= midplane optical depth; mu0= mu of beam 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 z=. z+ mu* Dt (1-extn)*rand n NB. z = coordinates of new scat point tau=: T0- z NB. GLOBALS: vertical optical depth from surface cfr=: 1- fr=: haze tau NB. fraction of scattering by type A vs type B S12b=: (cfr*S12Ab)+ fr*S12Bb NB. each scattering point has different S12b ! S=. S*"1 albedo'' NB. reduce intensity by single scat albedo remn=: remn, (n%~ +/0{S) NB. integrated intensity still in slab (global) 'S C'=. C SCAT S NB. scatter photons into new direction C S;z;C ) NB. --------------- just examples for testing ------------------------------ NB. this is just an example of a depth dependent function for the albedo: NB. the albedo is unity at the surface and decreases to 0.75 at bottom of slab albedo=: monad define fract=. tau% +:T0 1-0.25*fract ) NB. this is just a toy haze layer, where the layer is centered at tau=0.8 and NB. reaches a maximum of 50% of the opacity at that depth haze=: monad define NB. proportion of scattering due to 2nd component maxh* ^-(|tau- tauh)^6 ) tauh=: 0.8 maxh=: 0.5 NB. ------------------------------------------------------------------------ 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) 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=. (n0=. #y)#1 NB. inr are indices still to do phip=. n0#0 NB. placeholder for phi' results NB. p_count=: 0 NB. initialize counter while. n0>0 do. r=. 2p1* rand n0 NB. trial value of phi' b2=. ((inr{y);(inr{S12b)) 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" n0=. +/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 define 'xx s12b'=. x ((1+ |xx*s12b)*(rand #y))<: 1+ s12b*xx*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=. (TAB51A;angrad) lintrpX rd NB. thetas for 51 alphas at each rd TA=. (aa;toobig) lintrpY y NB. & interpolate for the correct alpha toobig=. (TAB51B;angrad) lintrpX rd NB. thetas for 51 alphas at each rd TB=. (aa;toobig) lintrpY y NB. & interpolate for the correct alpha (cfr*TA)+ fr*TB ) 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 cf=. 1- f=. haze tau NB. get relative fractions for each point s11=. cfr* (angrad;S11A;S11DA) spltrp x NB. 1st S11 scat matrix element s11=. s11+ fr* (angrad;S11B;S11DB) spltrp x NB. added to 2nd S11 element s12=. cfr* (angrad;S12A;S12DA) spltrp x NB. ditto for S12 elements ... s12=. s12+ fr* (angrad;S12B;S12DB) spltrp x s33=. cfr* (angrad;S33A;S33DA) spltrp x s33=. s33+ fr* (angrad;S33B;S33DB) spltrp x s34=. cfr* (angrad;S34A;S34DA) spltrp x s34=. s34+ fr* (angrad;S34B;S34DB) 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 PBM3_dir,ffff 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 1st file with tabulated scattering matrix: ' Smat=. enter'' file=. dirct,Smat NB. full path to scattering matrix 'angdeg S11A S12A S33A S34A'=: |:mread file NB. read in scattering matrix table angrad=: (1p1%180)*angdeg NB. convert angles to radians S11b=: -: angrad sinteg (S11A* sin angrad) NB. should be normalized to unity 'S11A S12B S33B S34B'=: S11b%~ S11A,S12A,S33A,:S34A S11DA=: angrad spline S11A NB. get 2nd derivatives for spline interpolation S12DA=: angrad spline S12A S33DA=: angrad spline S33A S34DA=: angrad spline S34A F11=: -:angrad spintg S11A* sin angrad NB. Cumulative S11A element S12Ab=: {: F12=: -:angrad spintg S12A* sin angrad NB. Cumulative integral of S12A aa=: _1++:int01 50 NB. Grid of alpha values (_1< alpha <1) eqn51=: dyad : '(1+y*S12Ab)%~ ((0{x)+y*(1{x))' NB. equation 51 solved for r TAB51A=: (>F11;F12)&eqn51"0 aa NB. Table of r as function of alpha & theta print 'S12Ab = '; S12Ab;' .... OK ' dirct=. '/your_path/' NB. directory with scattering matrix print 'Enter name of 2nd file with tabulated scattering matrix: ' Smat=. enter'' file=. dirct,Smat NB. full path to scattering matrix 'angdeg S11B S12B S33B S34B'=: |: mread file NB. read in scattering matrix table angrad=: (1p1%180)*angdeg NB. convert angles to radians S11b=: -: angrad sinteg (S11B* sin angrad) NB. should be normalized to unity 'S11B S12B S33B S34B'=: S11b%~ S11B,S12B,S33B,:S34B S11DB=: angrad spline S11B NB. get 2nd derivatives for spline interpolation S12DB=: angrad spline S12B S33DB=: angrad spline S33B S34DB=: angrad spline S34B F11=: -:angrad spintg S11B* sin angrad NB. Cumulative S11B element S12Bb=: {: F12=: -:angrad spintg S12B* sin angrad NB. Cumulative integral of S12B aa=: _1++:int01 50 NB. Grid of alpha values (_1< alpha <1) eqn51=: dyad : '(1+y*S12Bb)%~ ((0{x)+y*(1{x))' NB. equation 51 solved for r TAB51B=: (>F11;F12)&eqn51"0 aa NB. Table of r as function of alpha & theta print 'S12Bb = '; S12Bb;' .... 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 ) PBM3_dir=: '/your_path/' Cyc_Log=: PBM3_dir,'PBM3_CYC.log' load'/your_path/splines.ijs' NB. load spline routines NB. make polarization magnitude & angle a complex number for use NB. in viewmat. Usage: (polarization) cvec (angle) --> complex NB. I.e., viewmat (1{FracPol SSu) cvec (2{FracPol SSu) cvec=: dyad : 'x* (-sin a)j. cos a=. (1p1%180)*y' 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=. PBM3_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)