NB. This code evaluates the Rayleigh scattering from a semi-infinite atmosphere NB. illuminated by a beam of arbitrary polarization. Impliments equations given NB. by Abhyankar & Fymat, ApJ Supp, 23, 35 (1971). Usage: NB. 'I Q U V'=. (albedo;mu_in;FI,FQ,FU,FV) ARayf (mu_out's; phi_out's) NB. where the Stokes parameters of the incoming beam are (FI,FQ,FU,FV). load '/your_path/splines.ijs' ARayf=: dyad define 'om mu0 FF'=. x 'FI FQ FU FV'=. FF 'FL FR'=. -: FI (+,-) FQ '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'=. |: (N mm NT) mm (FL,FR) NB. H^(1)(mu), H^(2)(mu), H_V^(0)(mu) & H_V^(1)(mu): 'H1 H2 HV0 HV1'=. om H12f mu NB. H^(1)(mu_0), H^(2)(mu_0), H_V^(0)(mu_0) & H_V^(1)(mu_0): 'h1 h2 hv0 hv1'=. om H12f mu0 mu2=. *:mu mu02=. *:mu0 muf=. (%:1-mu02)*(%:1-mu2) E=. -:3* D=. om*mu0%4*mu0+mu z2=. 1.5*h1*mu*muf*H1 'sph cph'=. scfuns phi z2=. z2*/ -(FU*sph)+ +:FL*mu0*cph z4=. 0.75*h2*mu2*H2 'sph2 cph2'=. scfuns +:phi z4=. z4*/ ((-FR- mu02*FL)*cph2)+ mu0*FU*sph2 w2=. _0.75*h2*H2 w2=. w2*/ ((-FR- mu02*FL)*cph2)+ mu0*FU*sph2 u2=. h1*muf*H1 u2=. u2*/ (_2*FL*mu0*sph)+ FU*cph u4=. -h2*mu*H2 u4=. u4*/ ((FR- mu02*FL)*sph2)+ mu0*FU*cph2 I=. D* (Al+Ar)+ z2+ z4+ w2 Q=. D* (Al-Ar)+ z2+ z4- w2 U=. E* u2+ u4 V=. E*FV*(-mu0*hv0*mu*HV0)+ hv1*muf*HV1 */cph NB. Stokes parameters of scattered radiation: >I;Q;U;V ) scfuns=: [: > sin;cos 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 integral eqns for H^(1), H^(2), H_V^(0) and H_V^(1): H12f=: dyad define om=: x Cz=: %:1- 0.7*om Cq=: %:1- -: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 HV0=. mu& HV0_step ^:_ (#mu)#2 HV1=. mu& HV1_step ^:_ (#mu)#1 NB. Interpolate H(mu) grid for value(s) of mu requested: H1y=. (mu;H1) splint y H2y=. (mu;H2) splint y HV0y=. (mu;HV0) splint y >H1y;H2y;HV0y; (mu;HV1) 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. Do the HV^(0) integrals over mu' : HV0_step=: dyad define Phi=. (3*om%4)* *:x HH=. |: y*Phi*x%x +/x %Cq+ x sintegA HH ) NB. Do the HV^(1) integrals over mu' : HV1_step=: dyad define Phi=. (3*om%8)*(1- *:x) HH=. |: y*Phi*x%x +/x %Cq+ 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 and polarization of incoming NB. beam given by (1,_0.2,0.3,0.05) 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;(1,_0.2,0.3,0.05)) ARayf mubin ; |. abin NB. 'surface' plot mubin;abin; 1{W0