NB. Pure Rayleigh scattering of an unpolarized beam. NB. This gives the polarization and intensity of the NB. scattered radiation. The *first scattering only*! NB. Usage: mu_in Pure_RS theta_out, phi_out NB. input angles in radians; output is (I,Q,U) Pure_RS=: dyad define th_i=. acos mu_i=. x sin_i=. %:1- *:mu_i 'th_o ph_o'=. y mu_o=. cos th_o NB. calculate cosine of scattering angle, gamma: cos_G=. (mu_i*mu_o)+sin_i*(sin th_o)*cos ph_o sin_G=. sin gam=: acos cos_G cos2=. *: cos_G NB. cosine^2 of scattering angle I0=. 0.75* cos2+ 1 NB. intensity of scattered radiation Q0=. 0.75* cos2- 1 NB. Q wrt the plane of scattering cos_xi=. (mu_i- cos_G*mu_o)%(sin_G*sin th_o) sign=. 1+ _2*(ph_o >1p1) NB. negative if phi_out > 1p1 xi=. sign* acos cos_xi NB. angle between scat plane & meridian NB. rotate reference plane of I,Q,U to meridian: S=. I0, Q0* (cos +:xi),(-sin +:xi) S*(mu_o%(mu_i+mu_o)) NB. internal extinction of scattered ray ) NB. Compute Rayleigh scattering for a 2D array of emergent angles. NB. Output format is like that of "Pol_beam.ijs". Run_RS=: dyad define mu_in=. x 'theta phi'=. y nt=. #theta np=. #phi sc=. 2*nt*np NB. scale intensities by # bins SS=. (np,nt,3)$0 n1=. _1 while. (n1=. >:n1)< nt do. th=. n1{theta n2=. _1 while. (n2=. >:n2)< np do. ph=. n2{phi S=. mu_in Pure_RS th,ph SS=. (S%sc) (1{y)*1p1)+ 2p1+ atan %/y=. Re y' NB. I, Fractional polarization & angle from (I,Q,U) vector (or array) FracPol=: monad define pol=. POL y chi=. CHI y (0{y), pol,: (180%1p1)*chi ) NB. Example: NB. load 'plot' NB. th=. }. }: 0.5p1*int01 100 NB. define theta array NB. ph=. }. }: 1p1*int01 200 NB. phi array NB. W=. FracPol Z=. 0.7 Run_RS th;ph NB. 'surface' plot th;ph;(0{W) NB. plot of "intensity" from unit area NB. The specific intensity should include the increase in projected area: NB. 'surface' plot th;ph;((cos th)%~ 0{W) NB. Intensity plot NB. 'surface' plot th;ph;(1{W) NB. Polarization plot NB. 'surface' plot th;ph;(2{W)+((2{W)>90)*_180 NB. angle plot