NB. This code evaluates the Rayleigh scattering from an illuminated, NB. semi-infinite atmosphere using the analytic formulae as found in NB. ---> Madhusudhan & Burrows, ApJ 747, 25 (2012). NB. Usage: 'I Q U'=. (albedo;mu_in) ARay (mu_out's; phi_out's) load '/your_path/splines.ijs' ARay=: dyad define 'om mu0'=. x 'mu phi'=. y N=. om MNf mu NB. Get N(mu) matrices NT=. |: om MNf mu0 NB. Transpose of N(mu_0) matrix 'Al Ar'=. |: +/"1 N mm NT NB. sum of matrix product rows = Al & Ar 'H1 H2'=. om H12f mu NB. H^(1)(mu) & H^(2)(mu) 'h1 h2'=. om H12f mu0 NB. H^(1)(mu_0) & H^(2)(mu_0) mu2=. *:mu BU=. 3*mu0*(%:1-*:mu0)*h1*(%:1-mu2)*H1 Cr=. 0.75*(1-*:mu0)*h2*H2 CU=. 2*mu*Cr D=. om*mu0%8*mu0+mu NB. note Abhyankar & Fymat erattum Blc=. -mu*BU*/cos phi CC=. Cr*/cos +:phi I=. D*(Al+Ar)+ Blc+ (1-mu2)*CC Q=. D*(Al-Ar)+ Blc- (1+mu2)*CC U=. -D*(BU*/sin phi)+ CU*/sin +:phi >I;Q;U ) NB. Solve integral equations for N(mu) matrix(s) MNf=: dyad define om=: x NB. Define array of 44 mu values, emphasis near mu=0: mu=. (1e_6,0.0005*10^(10%~ i. 26)),0.2+0.05*i. 17 MMM=: MM mu N=. mu& MN_step ^:_ MMM NB. iterate equation to convergence. 'n1 n2 n3 n4'=. |: ,"2 N n1y=. (mu;n1) splint y n2y=. (mu;n2) splint y n3y=. (mu;n3) splint y n4y=. (mu;n4) splint y shape >n1y;n2y;n3y;n4y ) MN_step=: dyad define NT=. |:"2 y MP=. NT mm"2 MMM mx=. |: (2,2,(#x))$x Q=. MP%"2 (mx +"2/mx) 'a11 a12 a21 a22'=. |: ,"2 Q Q11=: x sintegA a11 Q12=: x sintegA a12 Q21=: x sintegA a21 Q22=: x sintegA a22 IM=. shape >Q11;Q12;Q21;Q22 MMM+ -:om*x* y mm"2 IM ) mcn=: -: %:3 NB. (1/2) sqrt(3) shape=: [:(2 2$,)"1 |: NB. form 4 items into 2x2 matrix MM=: monad define NB. construct M(omega) matrix m1m=. (%:2)* 1- m2=. *:y n=. # m2=. mcn*m2 shape >m2;(mcn*m1m);(n#mcn);n#0 ) NB. Solve the integral equations for H^(1) & H^(2): H12f=: dyad define om=: x Cz=: %:1- 0.7*om NB. Define array of 44 mu values, emphasis near mu=0: mu=. (0,0.0005*10^(10%~ i. 26)),0.2+0.05*i. 17 H1=. mu& H1_step ^:_ (#mu)#2 NB. iterate to convergence H2=. mu& H2_step ^:_ (#mu)#2 NB. Interpolate H(mu) grid for value(s) of mu requested: H1y=. (mu;H1) splint y >H1y; (mu;H2) splint y ) NB. Do the H^(1) integrals over mu' : H1_step=: dyad define Phi=. (3*om%8)*(1-x*x)*(1+ +:x*x) HH=. |: y*Phi*x%x +/x %Cz+ x sintegA HH ) NB. Do the H^(2) integrals over mu' : H2_step=: dyad define Phi=. (3*om%16)* *:(1+x*x) HH=. |: y*Phi*x%x +/x %Cz+ x sintegA HH ) NB. ---------------------------------------------------------------- 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 ) L=: +/ &. (*:"_) NB. Gives the lengths of vector arrays 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 }. atan2r=: 13 : '2p1|((0>1{y)*1p1)+ 2p1+ atan %/y=. Re y' NB. Angle of plane of polarization wrt reference plane of NB. (I,Q,U) vector or (3,n) array CHI=: [: -: [: atan2r 2 1 { ] NB. Example: case of omega=0.95, mu_in=0.5 NB. for grid of mu's and phi's used for Monte Carlo runs: NB. mubin=. }. (%_60)+ int01 30 NB. abin=. (1p1%180)* }. _2.5+ 180*int01 36 NB. W0=. FracPol Z0=. (0.95 0.5) ARay mubin ; |. abin NB. Plot the polarization: NB. 'surface' plot mubin;abin; 1{W0