pro smbmie, x,refrel,theta,nang,s1,s2,qext,qsca,qback,gsca ; WARNING your nang should be *smaller* by 1 than mxnang ; because of IDL array indexinf (from 0). ;*********************************************************************** ; Subroutine BHMIE is the Bohren-Huffman Mie scattering subroutine ; to calculate scattering and absorption by a homogenous isotropic ; sphere. ; Given: ; X = 2*pi*a/lambda ; REFREL = (complex refr. index of sphere)/(real index of medium) ; NANG = number of angles between 0 and 90 degrees ; (will calculate 2*NANG-1 directions from 0 to 180 deg.) ; if called with NANG<2, will set NANG=2 and will compute ; scattering for theta=0,90,180. ; Returns: ; S1(1 - 2*NANG-1) = -i*f_22 (incid. E perp. to scatt. plane, ; scatt. E perp. to scatt. plane) ; S2(1 - 2*NANG-1) = -i*f_11 (incid. E parr. to scatt. plane, ; scatt. E parr. to scatt. plane) ; QEXT = C_ext/pi*a**2 = efficiency factor for extinction ; QSCA = C_sca/pi*a**2 = efficiency factor for scattering ; QBACK = (dC_sca/domega)/pi*a**2 ; = backscattering efficiency [NB: this is (1/4*pi) smaller ; than the "radar backscattering efficiency"; see Bohren & ; Huffman 1983 pp. 120-123] ; GSCA = for scattering ; ; Original program taken from Bohren and Huffman (1983), Appendix A ; Modified by B.T.Draine, Princeton Univ. Obs., 90/10/26 ; in order to compute ; Converted to IDL by P. J. Flatau, Scripps Inst. Oceanography. UCSD 96/11/14 ; Modified by Henry Throop, U. Colorado, 99/04/10 ; ; 91/05/07 (BTD): Modified to allow NANG=1 ; 91/08/15 (BTD): Corrected error (failure to initialize P) ; 91/08/15 (BTD): Modified to enhance vectorizability. ; 91/08/15 (BTD): Modified to make NANG=2 if called with NANG=1 ; 91/08/15 (BTD): Changed definition of QBACK. ; 92/01/08 (BTD): Converted to full double precision and double complex ; eliminated 2 unneed lines of code ; eliminated redundant variables (e.g. APSI,APSI0) ; renamed RN -> EN = double precision N ; Note that DOUBLE COMPLEX and DCMPLX are not part ; of f77 standard, so this version may not be fully ; portable. In event that portable version is ; needed, use src/bhmie_f77.f ; 93/06/01 (BTD): Changed AMAX1 to generic function MAX ; 96/11/14 (PJF): Converted to double precision (again) ; added declaration standarization (strong typing) ; some code "polishing" using Toolpack ; Converted to IDL ; 99/04/10 (HBT): Vectorized several (not all) for-next loops. ; For X=1000, speed-up is ~ factor of 10. ;*********************************************************************** x = x*1d ; added by ht, error checking -- ints will confuse refrel = refrel*1d ; .. Parameters .. ; INTEGER mxnang, nmxx ; PARAMETER (mxnang=1000,nmxx=15000) ; mxnang=10000 nmxx=15000 ; .. Scalar Arguments .. ; DOUBLE COMPLEX refrel ; DOUBLE PRECISION gsca, qback, qext, qsca, x ; INTEGER nang ; .. Array Arguments .. ; DOUBLE COMPLEX s1(2*mxnang-1), s2(2*mxnang-1) ; s1 and s2 are defined here not in the subroutine calling "bhmie" s1=dcomplexarr(nang+1) s2=dcomplexarr(nang+1) ; .. Local Scalars .. ; DOUBLE COMPLEX an, an1, bn, bn1, drefrl, xi, xi1, y ; DOUBLE PRECISION chi, chi0, chi1, dang, dx, en, fn, p, pii, psi,$ ; psi0, psi1, theta, xstop, ymod ; INTEGER j, jj, n, nmx, nn, nstop ; .. Local Arrays .. ; DOUBLE COMPLEX d(nmxx) d=dcomplexarr(nmxx) ; DOUBLE PRECISION amu(mxnang), pi(mxnang), pi0(mxnang),$ ; pi1(mxnang), tau(mxnang) amu=dblarr(nang+1) pi=dblarr(nang+1) pi0=dblarr(nang+1) pi1=dblarr(nang+1) tau=dblarr(nang+1) ; .. Intrinsic Functions .. ; INTRINSIC abs, atan, cos, dble, dcmplx, max, sin ; IF (nang GT mxnang) then begin ; print, ' error: nang > mxnang in bhmie' ; stop ; endif IF (nang LT 2) then nang = 2 dx = x ;flatau I have converted everything to double so drefrl is really not needed drefrl = refrel y = x*drefrl ymod = abs(y) ; ;*** Series expansion terminated after NSTOP terms ; Logarithmic derivatives calculated from NMX on down xstop = x + 4.D0*x^0.3333D0 + 2.D0 nmx = max([xstop,ymod]) + 15 ;flatau IDL fix (convert to integer) nmx=fix(nmx) ; print, ' nmx ', nmx ; commented out by HT ; BTD experiment 91/1/15: add one more term to series and compare results ; NMX=AMAX1(XSTOP,YMOD)+16 ; test: compute 7001 wavelengths between .0001 and 1000 micron ; for a=1.0micron SiC grain. When NMX increased by 1, only a single ; computed number changed (out of 4*7001) and it only changed by 1/8387 ; conclusion: we are indeed retaining enough terms in series! nstop = xstop ; IF (nmx GT nmxx) THEN begin print, 'error: nmx > nmxx=', nmxx, ' for |m|x=', ymod STOP ENDIF amu[1:nang] = cos(theta*!pi/180.) pi0[1:nang] = 0.d0 ; ht pi1[1:nang] = 1.d0 ; ht s1[1:nang] = dcomplex(0.d0, 0.d0) s2[1:nang] = dcomplex(0.d0, 0.d0) ; ;*** Logarithmic derivative D(J) calculated by downward recurrence ; beginning with initial value (0.,0.) at J=NMX ; d(nmx) = dcomplex(0.D0,0.D0) nn = nmx - 1 for n=1, nn do begin ;DO 40 n = 1, nn en = nmx - n + 1 d(nmx-n) = (en/y) - (1.D0/ (d(nmx-n+1)+en/y)) endfor ;40 CONTINUE ;*** Riccati-Bessel functions with real argument X ; calculated by upward recurrence ; psi0 = cos(dx) psi1 = sin(dx) chi0 = -sin(dx) chi1 = cos(dx) xi1 = dcomplex(psi1,-chi1) qsca = 0.D0 gsca = 0.D0 ; ht -- the code could be sped up significantly ; by rewriting following loop over n into vector arithmetic. for n=1, nstop do begin; DO 80 n = 1, nstop en = n fn = (2.D0*en+1.D0)/ (en* (en+1.D0)) ; for given N, PSI = psi_n CHI = chi_n ; PSI1 = psi_{n-1} CHI1 = chi_{n-1} ; PSI0 = psi_{n-2} CHI0 = chi_{n-2} ; Calculate psi_n and chi_n psi = (2.D0*en-1.D0)*psi1/dx - psi0 chi = (2.D0*en-1.D0)*chi1/dx - chi0 xi = dcomplex(psi,-chi) ; ;*** Store previous values of AN and BN for use ; in computation of g= IF (n GT 1) THEN begin an1 = an bn1 = bn ENDIF ; ;*** Compute AN and BN: an = (d(n)/drefrl+en/dx)*psi - psi1 an = an/ ((d(n)/drefrl+en/dx)*xi-xi1) bn = (drefrl*d(n)+en/dx)*psi - psi1 bn = bn/ ((drefrl*d(n)+en/dx)*xi-xi1) ; ;*** Augment sums for Qsca and g= animag=double(an*complex(0.d0,-1.d0)) bnimag=double(bn*complex(0.d0,-1.d0)) qsca = qsca + (2.D0*en+1.D0)* (abs(an)^2+abs(bn)^2) gsca = gsca + ((2.D0*en+1.D0)/ (en* (en+1.D0)))* $ (double(an)*double(bn)+animag*bnimag) IF (n GT 1) THEN begin an1imag=double(an1*complex(0d0,-1.d0)) bn1imag=double(bn1*complex(0d0,-1.d0)) animag=double(an*complex(0d0,-1.d0)) bnimag=double(bn*complex(0d0,-1.d0)) gsca = gsca + ((en-1.D0)* (en+1.D0)/en)*$ (double(an1)*double(an)+an1imag*animag+$ double(bn1)*double(bn)+bn1imag*bnimag) ENDIF ; ;*** Now calculate scattering intensity pattern pi[1:nang] = pi1[1:nang] ; ht tau[1:nang] = en*amu[1:nang]*pi[1:nang] - (en+1.d)*pi0[1:nang] ; ht s1[1:nang] = s1[1:nang] + fn*(an*pi[1:nang] + bn*tau[1:nang]) ; smb s2[1:nang] = s2[1:nang] + fn*(an*tau[1:nang] + bn*pi[1:nang]) ; smb psi0 = psi1 psi1 = psi chi0 = chi1 chi1 = chi xi1 = dcomplex(psi1,-chi1) ; ;*** Compute pi_n for next value of n ; For each angle J, compute pi_n+1 ; from PI = pi_n , PI0 = pi_n-1 pi1[1:nang] = ((2.D0*en+1.D0)*amu[1:nang]*pi[1:nang] - (en+1.D0)*pi0[1:nang])/en ; smb pi0[1:nang] = pi[1:nang] ; smb endfor ; 80 CONTINUE ; ;*** Have summed sufficient terms. ; Now compute QSCA,QEXT,QBACK,and GSCA gsca = 2.D0*gsca/qsca qsca = (2.D0/ (dx*dx))*qsca qext = (4.D0/ (dx*dx))*double(s1(1)) qback = (abs(s1(nang))/dx)^2/!pi RETURN END