NB. Rayleigh scattering from a point source at x=y=0, z=z_s of slab. NB. Usage: (T; z_s; hnmu; nxy; dxy) PPaS 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. PPaS=: 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 frets trimmed Sorted=: 0$~(>: #&>bfr),3 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,3)$1 0 0 NB. Initial Stokes parameters (unpolarized) z=. k#z_s NB. initially emitted photons at z=z_s. P=. >0;0;z 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 z C P n0'=. Step^:nkeep S;z;C;P;0 NB. do "nkeep" escape/scatter cycles SS=: k%~ }. SS NB. trim off empty 1st item & normalize Step^:(nscat-nkeep) S;z;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=: trim upind{ Sorted NB. seperate upward and downward streams Sdn=: trim dnind{ Sorted NB. trim empty borders & 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)* Sdn Su=. trim2 upind{"4 SS NB. Stokes parameters of upward radiation Sd=. (1 1 _1)*"1 dnind{"4 SS NB. Stokes of downward, U fliped Sd=. trim2 Sd NB. trim & reorder axes mubin=. upind{(%nmu)+ }: fr1 NB. the value of mu at the bin centers xbins=. ({.fr2),(BIN fr2),{:fr2 NB. x-coordinates of bin centers + ends ybins=. (BIN fr3),{:fr3 NB. y-coordinates of bin centers + far end mubin;xbins;ybins;Sup;Sdn;Su;Sd ) NB. One cycle of escape and re-scatter. Step=: monad define 'S z C P key'=. y 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 Rout=. L }: Pout NB. radial distance from z-axis of escape point phi_out=. atan2 |: }: Pout NB. angular coordinate of escape point phi_ray=. atan2 |: }: C NB. angular coordinate of escape trajectory rphi=. 2p1| phi_out - phi_ray NB. angle between position & trajectory xout=. Rout* cos rphi NB. x-coordinate of escape point yout=. Rout* sin rphi NB. y-coordinate 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 z=. 2{ P=. P+ Dtau *"1 C NB. new position of next scattering 'S C'=. C SCAT S S;z;C;P;key ) 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=: dyad : '(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=: 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 idx=. <"1 ~. b w=. b +//. |: >{:y Shape=. 0$~(>:#&>x),3 NB. |:"_1 |: w idx} Shape w idx} Shape ) trim=: monad define NB. fix edges & reorder for use with FracPol 'a0 a1 a2'=. |: y NB. transpose and assign a0=. }. (+/><0 1{a0) 1}a0 NB. add 1st (y<0) bin to 2nd & trim off 1st a1=. }. (+/><0 1{a1) 1}a1 a2=. }. (+/><0 1{a2) 1}a2 a0=. |: a0 NB. join to reflection in x-axis, undo transpose a1=. |: a1 a2=. |: a2 >a0;a1;a2 NB. join to make (3, mu, nx, ny) array ) trim2=: monad define NB. fix edges & reorder for use with FracPol 'a0 a1 a2'=. |: y NB. transpose and assign a0=. }. (+/><0 1{a0) 1}a0 NB. add 1st (y<0) bin to 2nd & trim off 1st a1=. }. (+/><0 1{a1) 1}a1 a2=. }. (+/><0 1{a2) 1}a2 a0=. |: a0 NB. join to reflection in x-axis, undo transpose a1=. |: a1 a2=. |: a2 >a0;"_1 a1;"_1 a2 NB. join to make (nkeep,3,mu,nx,ny) array ) 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)* 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: PPaS_CYC=: monad define 'T z_s hnmu nxy dxy nphoton nkeep nscat ncycle'=. y zres=. (T;z_s;hnmu;nxy;dxy) PPaS nphoton;nkeep;nscat 'mubin xbins ybins SSup SSdn SSu SSd'=. zres file=. 'PPaS_',(":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) PPaS 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 ) PPaS_dir=: '/your_path/' Cyc_Log=: PPaS_dir,'PPaS_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=. PPaS_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 xbins ybins 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