subroutine BHMIE(X,REFREL,NANG,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) ! 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] ! [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 ! ! 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 :: MXNANG = 1000 !!$ integer, parameter :: NMXX = 15000 integer, parameter :: NMXX = 150000 ! Arguments: real, intent (in) :: X complex, intent (in) :: REFREL integer, intent (in) :: NANG complex, intent (out) :: S1(2*MXNANG-1), S2(2*MXNANG-1) real, intent (out) :: QEXT, QSCA, QBACK, GSCA ! Local variables: 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 PII = 4.0D0 * atan(1.0d0) ! Obtain pi ! Safety checks if (NANG .gt. MXNANG) stop 'Error: NANG > MXNANG in bhmie' if (NANG .lt. 2) stop 'Error: NANG must be greater than 1.' 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 + 4.0 * X**0.3333 + 2.0 NMX = 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 = XSTOP if (NMX .gt. NMXX) then write(0,*) 'Error: NMX > NMXX=', NMXX,' for |m|x=', YMOD stop endif ! Require NANG.GE.1 in order to calculate scattering intensities DANG = 0.0 if (NANG .gt. 1) DANG = 0.5 * PII / dble(NANG - 1) do J = 1, NANG THETA = dble(J - 1) * DANG AMU(J) = cos(THETA) enddo do J = 1, NANG PI0(J) = 0.0 PI1(J) = 1.0 enddo NN = 2*NANG - 1 do J = 1, NN S1(J) = (0.0, 0.0) S2(J) = (0.0, 0.0) enddo ! Logarithmic derivative D(J) calculated by downward recurrence ! beginning with initial value (0.,0.) at J=NMX D(NMX) = (0.0, 0.0) NN = NMX - 1 do N = 1, NN EN = NMX - N + 1 D(NMX-N) = (EN/Y) - (1.0/(D(NMX-N+1) + EN/Y)) enddo ! Riccati-Bessel functions with real argument X ! calculated by upward recurrence PSI0 = cos(DX) PSI1 = sin(DX) CHI0 = -sin(DX) CHI1 = cos(DX) !!$ XI1 = cmplx(PSI1, -CHI1, kind = kind(1.0d0)) XI1 = cmplx(PSI1, -CHI1) QSCA = 0.0E0 GSCA = 0.0E0 P = -1.0 do N = 1, NSTOP EN = N FN = (2.0E0*EN + 1.0) / (EN*(EN + 1.0)) ! 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.0E0*EN - 1.0) * PSI1/DX - PSI0 CHI = (2.0E0*EN - 1.0) * CHI1/DX - CHI0 !!$ XI = cmplx(PSI, -CHI, kind = kind(1.0d0)) XI = cmplx(PSI, -CHI) ! Store previous values of AN and BN for use ! in computation of g= if (N .gt. 1) then 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 + (2.0*EN + 1.0) * (abs(AN)**2 + abs(BN)**2) GSCA = GSCA + ((2.0*EN + 1.0)/(EN*(EN + 1.0))) * & (real(AN)*real(BN) + aimag(AN)*aimag(BN)) if (N .gt. 1) then GSCA = GSCA + ((EN - 1.0) * (EN + 1.0) / EN) * & (real(AN1)*real(AN) + aimag(AN1)*aimag(AN) + & real(BN1)*real(BN) + aimag(BN1)*aimag(BN)) endif ! Now calculate scattering intensity pattern ! First do angles from 0 to 90 do J = 1, NANG 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)) enddo ! 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 do J = 1, NANG-1 JJ = 2*NANG - J S1(JJ) = S1(JJ) + FN*P*(AN*PI(J) - BN*TAU(J)) S2(JJ) = S2(JJ) + FN*P*(BN*PI(J) - AN*TAU(J)) enddo PSI0 = PSI1 PSI1 = PSI CHI0 = CHI1 CHI1 = CHI !!$ XI1 = cmplx(PSI1, -CHI1, kind = kind(1.0d0)) XI1 = cmplx(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 do J = 1, NANG PI1(J) = ((2.0*EN + 1.0)*AMU(J)*PI(J) - (EN + 1.0)*PI0(J))/EN PI0(J) = PI(J) enddo enddo ! Have summed sufficient terms. ! Now compute QSCA,QEXT,QBACK,and GSCA GSCA = 2.0 * GSCA / QSCA QSCA = (2.0/(DX*DX)) * QSCA QEXT = (4.0/(DX*DX)) * real(S1(1)) QBACK = (abs(S1(2*NANG - 1)) / DX)**2 / PII return end subroutine BHMIE