NB. Mie scattering from a point source at x=y=0, z=z_s of slab. NB. Usage: (T; z_s; hnmu; nxy; dxy) PPaSmix 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 scatterings. We need to set nkeep < nscat. NB. The albedo is unity - other albedos can recovered from seperate scatterings. PPaSmix=: dyad define '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 and y-axis trimmed) Sorted=: 0$~(>(1+#)&.>bfr),4 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,4)$1 0 0 0 NB. Initial Stokes parameters (unpolarized) P=. >0;0;(k#z_s) 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 C P n0'=. Step^:nkeep S;C;P;0 NB. do "nkeep" escape/scatter cycles SS=: k%~ }. SS NB. trim off empty 1st item & normalize Step^:(nscat-nkeep) S;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=: |:"3 |: upind{ Sorted NB. seperate upward and downward streams Sdn=: |:"3 |: dnind{ Sorted NB. trim empty borders & reflect in x-axis Su=. |:"3 |:"_1 upind{"4 SS NB. Stokes parameters of upward radiation Sd=. |:"3 |:"_1 dnind{"4 SS NB. Stokes of downward, U fliped mubin=. upind{(%nmu)+ }: fr1 NB. the value of mu at the bin centers xbin=. ({.fr2),(BIN fr2),{:fr2 NB. x-coordinates of bin centers + ends ybin=. (BIN fr3),{:fr3 NB. y-coordinates of bin centers + far end mubin;xbin;ybin;Sup;Sdn;Su;Sd ) NB. One cycle of escape and re-scatter. Step=: monad define 'S C P key'=. y z=. {: P NB. z coordinate of scattering points 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 phi_ray=. atan2 |: |. }: C NB. angular coordinate of escape trajectory QQ=.(QY-acos mu)QM (QZ-phi_ray) NB. quaternion to rotate about z, then y Pnew=. QQ Pr Pout NB. rotate emergent points about z, then y 'xout yout'=. }: Pnew NB. xy-coordinates 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 P=. P+ Dtau *"1 C NB. new position of next scattering 'S C'=. C SCAT S S;C;P;key ) 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 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. 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 scattering NB. 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. 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. For rotations about y (& z) axis: QY angles --> rotation quaternions QY=: [:>([:cos-:);0;0;~[:sin-: NB. quaternion for rotation about y-axis QZ=: [:>([:cos-:);0;0;[:sin-: NB. quaternion for rotation about z-axis 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=: dyad : '}. 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=: monad : '(1- +:1p1 }:y NB. idx=. <"1 ~. b NB. w=. b +//. |: >{:y NB. Shape=. 0$~(>:#&>x),3 NB. |:"_1 |: w idx} Shape NB. w idx} Shape NB. ) NB. Given a disordered array A of results depending on 2 variables, sort NB. it into bins where the bin edges are provided in the arrays fr1, fr2 & fr3 NB. while summing the contents of each bin. The vectors a1, a2 & a3 are the NB. variables attached to the values of A by which we sort them. NB. The output has dimensions (1+#fr1,1+#fr2,1+#fr3). NB. Usage: Sp=. (fr1;fr2;fr3) Sort a1;a2;a3;A Sort=: dyad define b=. |: > x I. &.> }:y NB. get cell indicies of each a1,a2,a3 triplet 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 ) 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=: monad define ar=. |."1 (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=: dyad : 'x* (-sin a)j. cos a=. (1p1%180)*y' NB. a_series=: dyad : '(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:%] 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 "PPaS" ncycle times: PPaSmix_CYC=: monad define Setup'' 'T z_s hnmu nxy dxy nphoton nkeep nscat ncycle'=. y zres=. (T;z_s;hnmu;nxy;dxy) PPaSmix nphoton;nkeep;nscat 'mubin xbins ybins SSup SSdn SSu SSd'=. zres file=. 'PPaSmix_',(":T),'_',(":z_s),'_',(":hnmu),'_',(":nxy),'_',(":dxy) file=. file,'_',(":nphoton),'_',(":nkeep),'_',(":nscat),'_',(":ncycle),'.dat' ((":1),LF) (1!:2) < Cyc_Log (zres=. 1;zres;remn) writeB PPaS_dir,file remain=. remn count=. 2 while. count<: ncycle do. zres=. (T;z_s;hnmu;nxy;dxy) PPaSmix nphoton;nkeep;nscat 'mubin xbins ybins Sup Sdn Su Sd'=. zres SSup =. SSup + Sup SSdn =. SSdn + Sdn SSu=. SSu + Su SSd=. SSd + Sd remain=. remain + remn zres=. count;mubin;xbins;ybins;count%~ &.> SSup;SSdn;SSu;SSd;remain zres writeB PPaS_dir,file ((": 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 load'/your_path/splines.ijs' NB. load spline routines 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 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. (}.}: t) I. x0 NB. linear interpolation: (x;y) lintrpX 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 ) PPaS_dir=: '/your_path/' Cyc_Log=: PPaS_dir,'PPaSmix_CYC.log' 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=. PPaSmix_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. 'count mubin xbin ybin 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. load'plot' 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=: dyad define 'n sl su'=. x sc=. sl+ (su-sl)*int01 n print 'Scale Levels: ';sc np=. n+1 tt=. y k=.0 while. k