NB. This code computes values of the Mie scattering matrix NB. elements S_11, S_12, S_33 and S_34, averaged over one NB. of two power-law distributions of grain radius "a". NB. Usage: n Power_law 'file' -- here the integer n NB. picks the type of distribution: n=1 sets a sharp NB. upper radius a_max, while n=2 provides an exponential NB. roll off. 'file' is the name of a ASCII file with NB. the real and imaginary parts of the grain's dielectric NB. constant as a function of energy in Ev. See appendix of NB. Harrington, Lame, White, & Borkowski (1997) AJ,113,2147 NB. for the equations defining the grain distributions. load '/your_path/splines.ijs' Power_law=: dyad define file=. '/your_path/',y 'evs eps1 eps2'=. |: mread file NB. Use lots of grain sizes to smooth over Mie variations with radius! nbin=. 5000 if. x=1 do. print 'Enter min grain radius (a_min) and' print 'maximum radius (a_max) in microns:' 'A1 A2'=. enter'' print 'Enter gamma for M(a)=a**(-gamma):' Gam=. enter'' FM=. Gam % (A1^-Gam)-(A2^-Gam) abar=. (FM%1-Gam)*((A2^1-Gam)-(A1^1-Gam)) print ' Mean grain radius = ',(":abar),' microns' print ' m(a) = const* a^(-',(":Gam),'), for ',(":A1),' < a < ',(":A2) delt=. (^.A2)- ^.A1 rayn=. aray=. ^(^.A1)+delt*int01 nbin-1 elseif. x=2 do. print 'Enter min grain radius (a_min) and' print 'rolloff radius a_b in microns:' 'A1 Ab'=. enter'' print 'Enter gamma for M(a)=a**(-gamma)*exp(-a/a_b):' Gam=. enter'' X0=. A1%Ab FM=. % GamInc (-Gam),X0 abar=. Ab*FM* GamInc (1-Gam),X0 print ' Mean grain radius = ',(":abar),' microns' print ' m(a) = exp[-a/',(":Ab),']*a^(-',(":Gam),'), for a > ',(":A1) dell=. (nbin-1)%~ (^. 4.6*Ab)- ^.A1 aray=. ^ (^.A1)+(1+i. (nbin-2))*dell rayn=. Ab%~ aray=. A1,aray,4.6*Ab NB. Re following 2 distributions: Seager "Exoplanet Atmospheres" p. 72. elseif. x=3 do. NB. distribution for H2O clouds, where a0= ~4 microns. print 'Enter droplet scale radius a0:' a0=. enter'' FM=. 8748%35 abar=. 1.5*a0 print ' Peak grain radius = ',(":abar),' microns' print ' N(a) = const*exp[-6a/',(":a0),']*(a/',(":a0),')^6 , for all a' aray=. a0* rayn=. }. 4.5*int01 nbin elseif. x=4 do. NB. stratospheric aerosols print 'Enter aerosol scale radius (a_0):' a0=. enter'' abar=. 18*a0 FM=. 8%315 print ' Peak grain radius = ',(":abar),' microns' print ' N(a) = const*(a/',(":a0),')*exp[-2(a/',(":a0),')^(1/2)], for all a' aray=. a0* rayn=. }. 120*int01 nbin end. print 'Enter wavelength (A): ' E=. 12398% waveA=. enter'' NB. energy in Ev corresponding to wavelength wave_mic=. 1.2398%E NB. wavelength in microns XX=. 2p1*aray%wave_mic NB. dimensionless parameters x = (2 pi lambda)/a ep1=. (evs;eps1) splint E NB. spline interpolation is table 'file' ep2=. |(evs;eps2) splint E NB. Qscat expects positive imaginary part cdie=. ep1 + (ep2*0j1) NB. cdie is a complex number print 'lambda = ',(":wave_mic),' microns, E = ',(":E),' eV, eps = ',(":cdie) SS=. cdie&Qscat"0 XX NB. get the Mie solutions for each x 'Qe Qs Qa'=. |: >>(0{"1 SS) 'S11 S12 S33 S34'=. |: }."1 SS NB. now integrate over all (nbin=1000) grain radii. The must be NB. weighted by the geometric area of dust as a function of "a": if. x=1 do. G=. % aray^(Gam+1) NB. dust area goes as exp[-(gamma+1)] elseif. x=2 do. G=. %(^rayn)*rayn^Gam+1 NB. exp[-a/a_b] * (a/a_b)^[-(gamma+1)] elseif. x=3 do. G=. (^_6*rayn)*rayn^8 NB. exp[-6a/a0] * (a/a0)^8 elseif. x=4 do. G=. (^_2* %:rayn)*rayn^3 NB. (a/a0)^3 * exp[-2(a/a0)^(1/2)] end. QBe=: FM* (rayn sinteg Qe*G) NB. sinteg = spline integration QBs=: FM* (rayn sinteg Qs*G) QBa=: FM* (rayn sinteg Qa*G) albd=: QBs%QBe p1=. 'Q_ex = ',(":QBe),', Q_sc = ',(":QBs),', Q_ab = ',(":QBa) print p1,', albedo = ',(":albd) NB. The average scattering matrix elements must also be weighted by NB. the dimensionless scattering cross-section Qs: G=. G*Qs QBS11=. FM*(rayn sinteg S11*G) QBS12=. FM*(rayn sinteg S12*G) QBS33=. FM*(rayn sinteg S33*G) QBS34=. FM*(rayn sinteg S34*G) NB. th=. 180*int01 180 NB. the 181 scattering angles used by "Qscat" th=. 180*int01 900 NB. the 901 scattering angles used by "Qscat" > th;QBS11;QBS12;QBS33;QBS34 ) NB. Compute the matrix elements for polarized scattering as a NB. function of 181 scattering angles "theta". NB. Usage: (diel_const) ang_scat (x = 2p1*a%lambda) NB. For this routine, diel_const should have *positive* imag. part Qscat=: dyad define ci=. 0j1 cm=. %: x cnrm=. cm sphere yy=. y nci=. >: nc n=. qext=. 0 NB. evaluate total extinction: while. (n=. >:n) <: nc do. nm=. <: n cim=. % ci^(>: n) rf=. (1+2*n)%(n*n+1) z=. cim*(nm{f + ci* nm{g) qext=. qext + rf* 1{ +. z end. qext=. qext% *: yy n=. t=. 0 NB. evaluate total scattering: while. (n=. >:n) <: nc do. nm=. <: n z1=. *: | nm{f z2=. *: | nm{g t=. t + (nm{cnrm)*z1+z2 end. qsca=. t % *: yy qabs=. qext - qsca NB. th=. 180*int01 180 NB. 181 angles: each degree 0,1,2, ... 179,180 th=. 180*int01 900 NB. 901 angles: each degree 0,0.2,0.4, ... 179.8,180 theta=. (1p1%180)*th costh=. cos theta pnmllg=. nci genlgp theta n=. s1=. s2=. 0 while. (n=. >:n) <: nc do. nm=. <: n cim=. % ci^(>: n) p1=. (n*costh*n{pnmllg)-((n+1)*(n-1){pnmllg) p2=. n{pnmllg s1=. s1 + cim*((p2* nm{g)+(ci*p1* nm{(-f)))* nm{cnrm s2=. s2 + cim*((p1* nm{g)-(ci*p2* nm{f))* nm{cnrm end. norm=. 4%t NB. t=(x^2 Qscat); norm so 1/2 integral S11(mu) dmu = 1 S11=. -:((*: |s2)+(*: |s1))*norm NB. ref. Bohren & Huffman 4.77 S12=. -:((*: |s2)-(*: |s1))*norm NB. note that (+s1) is the complex conjuate of s1, etc. S33=. -:Re(((+s2)*s1)+(s2*(+s1)))*norm S34=. -:Im((s1*(+s2))-(s2*(+s1)))*norm QQ=. qext;qsca;qabs QQ;S11;S12;S33;S34 ) sphere=: dyad define ci=. 0j1 nc=: <. xc =. 2+y+4.05*y^%3 z=. x*y nmx=. 15 + <.(xc >. |z) amat=. 0 n=. nmx while. n > nc do. amat=. (n%z)-(%amat+n%z) n=. <: n end. n=. nc while. n > 1 do. amat=. amat,((n%z)-(%(_1{amat)+n%z)) n=. <: n end. amat=. |. amat hkl=. (>: nc) besh y bj=. 0{ +. 0{ hkl fl=. gl=. cnrml=. n=. 1 while. n <: nc do. rf=. 2*n*n+1 bjm=. bj bj=. 0{ +. n{hkl b=. (n%y)+ x*(n-1){amat fl=. fl,((-ci^n)*rf*((-bjm)+b*bj)%((-(n-1){hkl)+b*n{hkl)) b=. (n%y)+(((n-1){amat)%x) gl=. gl,((ci^n+1)*rf*((-bjm)+b*bj)%((-(n-1){hkl)+b*n{hkl)) cnrml=. cnrml,(1+2*n)%(rf*n*n+1) n=. >: n end. f=: }. fl g=: }. gl cnrm=: }. cnrml ) NB. Generate the Hankel functions besh=: dyad define a=. (1 o. y)%y by=. -(2 o. y)%y by=. by,((by%y)-a) n=. 1 while. n < (x-1) do. z=. (1+2*n)*(_1{by)%y by=. by,(z - _2{by) n=. >: n end. nst=. x+ <. %:y +101 t3=. 0 t2=. 1e_35 i=. nst - 1 while. i > (x- 2) do. t1=. ((1+2*i)*t2%y)-t3 t3=. t2 t2=. t1 i=. <: i end. bj=. t3,t2 i=. (x- 2) while. i > 0 do. z=. (1+2*i)*(_1{bj)%y bj=. bj,(z - _2{bj) i=. <: i end. bj=. |. bj ((a% 0{bj)*bj)+0j1*by ) NB. Generate associated Legendre functions. NB. Usage: (number in series) genlgp (angle in radians) genlgp=: dyad define m=. #y costh=. 2 o. y pnmllg =. > (m#0) ; (m#1) n=. 2 while. n < x do. a=. (_1+2*n)*costh*_1{pnmllg pnmllg =. pnmllg,((%n-1)*a - n*_2{pnmllg) n=. >: n end. pnmllg ) NB. Incomplete Gamma function. NB. Usage GamInc a,x --> Gamma[a,x] GamInc=: monad define 'a t'=. y b0=: g=. 0 a0=: b1=: fac=: 1 a1=: t n_gam=: 1 g=. {:Gam_Iter^:_ a,t,g g* ^(-t)+ a*^.t ) Gam_Iter=: monad define 'a t g'=. y a0=: fac*a1+ a0*an=. n_gam - a b0=: fac*b1+ b0*an anf=. n_gam*fac a1=: (t*a0)+ anf*a1 b1=: (t*b0)+ anf*b1 g=. b1*fac=: %a1 n_gam=: >:n_gam a,t,g ) NB. Example: NB. load'Power_law.ijs' NB. C=. 2 Power_law 'BE1' NB. Enter min grain radius (a_min) and NB. rolloff radius a_b in microns: NB. 0.001 0.1 NB. Enter gamma for M(a)=a**(-gamma)*exp(-a/a_b): NB. 0.2 NB. Mean grain radius = 0.0167361 microns NB. m(a) = exp(-a/0.1)*a^(-0.2), for a > 0.001 NB. Enter wavelength (A): NB. 5000 NB. lambda = 0.5 microns, E = 2.4796 eV, eps = 3.43663j3.21454 NB. Q_ex = 0.303789, Q_sc = 0.0777388, Q_ab = 0.22605, albedo = 0.255897 NB. 'theta S11 S12 S33 S34'=. C NB. 'pensize 3'plot theta;>(S11+S12);(S11-S12);(-S12%S11) NB. (This last plots scattering of radiation polarized perpendicular and parallel NB. to the plane of scattering, and the fractional polarization, as function of NB. the scattering angle "theta".)