PRO BHMIE, X,REFREL,ANGLES,S1,S2,QEXT,QSCA,QBACK,GSCA ;*********************************************************************** ; 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) ; ANGLES : The angles at which you would like the S1, S2 matrices evaluated. They may be unsorted. ; ; Returns: ; S1 The S1 amplitude scattering matrix evaluted at ANGLES ; = -i*f_22 (incid. E perp. to scatt. plane, ; scatt. E perp. to scatt. plane) ; S2 The S2 amplitude scattering matrix evaluted at ANGLES ; = -i*f_11 (incid. E parr. to scatt. pl ane, ; scatt. E parr. to scatt. p lane) ; 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] ; [BTJ: perhaps 1/(4*pi) smaller than equation 4.82, page 121] ; 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 Fortran 90 by Michael A. Walters, SSEC, UW-Madison ; January 18, 2002. Also changed modified error messages for use ; with web program ; ; 2004ish : Converted to IDL by Chris O'Dell ; ; Modified to accept unequally spaced angles. ; ; All calculations done in double precision ; ; 02/01/18 (MAW): Converted to Fortran 90 ; 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 ;^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^* ; implicit none ; Declare parameters: ; Note: important that MXNANG be consistent with dimension of S1 and S2 ; in calling routine; ; integer, parameter :: ; Arguments: ; real, intent (in) :: X ; complex, intent (in) :: REFREL ; integer, intent (in) :: NANGIN ; complex, intent (out) :: S1(2*MXNANG-1), S2(2*MXNANG-1) ; real, intent (out) :: QEXT, QSCA, QBACK, GSCA SMALL = 1e-6 ZERO = 0.0d ONE = 1.0d TWO = 2.0d FOUR = 4.0d HALF = 0.5d PII = !dpi angles_ = angles nai = n_elements(angles) if abs(min(angles)) GT SMALL then angles_ = [angles_, ZERO] if abs(max(angles) - 180.0d) GT SMALL then angles_ = [angles_, 180.0d] na = n_elements(angles) sortang = sort(angles_) angles_ = angles_(sortang) fw = where(angles_ LE 90.0d, nf) bw = where(angles_ GT 90.0d, nb) nfi = nf ; initial number of front-side angles front = angles_[fw] flip = 180.0d - angles_[bw] if nf eq 1 then v = lonarr(nb) else v = value_locate(front, flip) whave = where( abs(front[v] - flip) LT SMALL, nhave, ncomp=nnot, comp=wnot ) ; these guys I already have wb = v*0 if nhave GT 0 then wb[whave] = v[whave] if nnot GT 0 then begin front = [front, flip[wnot]] ; new front angles wb[wnot] = lindgen(nnot) + nf nf = nf + nnot ; final number of front-side elements endif S1f = dcomplexarr(nf) S2f = S1f S1b = dcomplexarr(nb) S2b = S1b ; Local variables: ; integer :: NANG ; integer :: J, JJ, N, NSTOP, NMX, NN ; double precision :: CHI,CHI0,CHI1,DANG,DX,EN,FN,P,PII,PSI,PSI0,PSI1, & ; THETA,XSTOP,YMOD ; double precision :: AMU(MXNANG),PI(MXNANG),PI0(MXNANG),PI1(MXNANG), & ; TAU(MXNANG) ; complex (kind = kind(1.0d0)) :: AN, AN1, BN, BN1, DREFRL, XI, XI1, Y ; complex (kind = kind(1.0d0)) :: D(NMXX) ; character (len = 64) :: message DX = X DREFRL = REFREL Y = X * DREFRL YMOD = abs(Y) ; Series expansion terminated after NSTOP terms ; Logarithmic derivatives calculated from NMX on down XSTOP = X + FOUR * X^0.3333333d0 + 2.0d0 NMX = fix(max(XSTOP,YMOD) + 15) ; BTD experiment 91/1/15: add one more term to series and compare result ; 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 = fix(XSTOP) if (NMX GT 150000) then begin print, 'Error: NMX > NMXX=', 150000,' for |m|x=', YMOD stop endif ; Require NANG.GE.1 in order to calculate scattering intensities AMU = cos(front * PII/180.0d) PI0 = ZERO + dblarr(nf) PI1 = PI0 + ONE ; Logarithmic derivative D(J) calculated by downward recurrence ; beginning with initial value (0.,0.) at J=NMX D = dcomplexarr(nmx+1) D[NMX] = dcomplex(ZERO, ZERO) NN = NMX - 1 for N = 1, NMX-1 do begin last = NMX - N + 1 D[last-1] = last/Y - (ONE/(D[last] + last/Y)) endfor ; note, D[0] is unused in this algorithm. ; Riccati-Bessel functions with real argument X ; calculated by upward recurrence PSI0 = cos(DX) PSI1 = sin(DX) CHI0 = -PSI1 CHI1 = PSI0 ;;$ XI1 = cmplx(PSI1, -CHI1, kind = kind(1.0d0)) XI1 = dcomplex(PSI1, -CHI1) QSCA = ZERO GSCA = ZERO P = -ONE for N = 1, NSTOP do begin EN = N FN = (TWO*EN + ONE) / (EN*(EN + ONE)) ; 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 = (TWO*EN - ONE) * PSI1/DX - PSI0 CHI = (TWO*EN - ONE) * CHI1/DX - CHI0 ;;$ XI = cmplx(PSI, -CHI, kind = kind(1.0d0)) 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= QSCA = QSCA + (TWO*EN + ONE) * (abs(AN)^2 + abs(BN)^2) GSCA = GSCA + ((TWO*EN + ONE)/(EN*(EN + ONE))) * $ (real(AN)*real(BN) + imaginary(AN)*imaginary(BN)) if (N GT 1) then begin GSCA = GSCA + ((EN - ONE) * (EN + ONE) / EN) * $ (real(AN1)*real(AN) + imaginary(AN1)*imaginary(AN) + $ real(BN1)*real(BN) + imaginary(BN1)*imaginary(BN)) endif ; Now calculate scattering intensity pattern ; First do angles from 0 to 90 ; for J = 1, NANG do begin ; JJ = 2*NANG - J ; PI(J) = PI1(J) ; TAU(J) = EN*AMU(J)*PI(J) - (EN + 1.0)*PI0(J) ; S1(J) = S1(J) + FN*(AN*PI(J) + BN*TAU(J)) ; S2(J) = S2(J) + FN*(AN*TAU(J) + BN*PI(J)) ; endfor PI = PI1 TAU = EN * AMU * PI - (EN+ONE)*PI0 S1f = S1f + FN*( AN*PI + BN*TAU ) S2f = S2f + FN*( AN*TAU + BN*PI ) ; Now do angles greater than 90 using PI and TAU from ; angles less than 90. ; P=1 for N=1,3,...; P=-1 for N=2,4,... P = -P PIR = PI[wb] ; corresponding frontside elements TAUR = TAU[wb] ; corresponding frontside elements S1b = S1b + FN*P*(AN*PIR - BN*TAUR) S2b = S2b + FN*P*(BN*PIR - AN*TAUR) PSI0 = PSI1 PSI1 = PSI CHI0 = CHI1 CHI1 = CHI ;;$ XI1 = cmplx(PSI1, -CHI1, kind = kind(1.0d0)) 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 = ((TWO*EN + ONE)*AMU*PI - (EN + ONE)*PI0)/EN PI0 = PI endfor ; Have summed sufficient terms. ; Now compute QSCA,QEXT,QBACK,and GSCA ;print, 'gsca =', gsca ;print, 'qsca = ', qsca if x LT 0.05d0 then gsca = ZERO else $ GSCA = TWO * GSCA / QSCA QSCA = (TWO/(DX*DX)) * QSCA QEXT = (FOUR/(DX*DX)) * real(S1f[0]) S1 = [S1f[0:(nfi-1)], S1b] S2 = [S2f[0:(nfi-1)], S2b] QBACK = (abs(S1[na-1]) / DX)^2 / PII ; unsort everything in angle space unsort = (sort(sortang))[0:(nai-1)] S1 = S1[unsort] S2 = S2[unsort] end ; subroutine BHMIE