module bhmie use math, only: pi implicit none integer, parameter :: & mxnang = 1000, & ! Note: important that MXNANG be consistent with dimension of S1 and S2 in calling routine! s_size = 2 * MXNANG - 1 save contains SUBROUTINE BH_MIE(X, REFREL, NANG, S1, S2, QEXT, QSCA, QBACK, GSCA) implicit none ! Declare parameters: INTEGER NMXX ! PARAMETER(MXNANG=1000,NMXX=15000) PARAMETER(NMXX = 150000) ! Arguments: INTEGER NANG doubleprecision GSCA, QBACK, QEXT, QSCA, X COMPLEX(kind=8) REFREL COMPLEX(kind=8) S1(s_size), S2(s_size) ! Local variables: INTEGER J, JJ, N, NSTOP, NMX, NN DOUBLE PRECISION CHI, CHI0, CHI1, DANG, DX, EN, FN, P, PSI, PSI0, PSI1, & THETA, XSTOP, YMOD DOUBLE PRECISION AMU(MXNANG), pi_(MXNANG), PI0(MXNANG), PI1(MXNANG), & TAU(MXNANG) COMPLEX(kind=8) AN, AN1, BN, BN1, DREFRL, XI, XI1, Y COMPLEX(kind=8) D(NMXX) !*********************************************************************** ! 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 ! 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 !*********************************************************************** !*** Safety checks IF(NANG>MXNANG)STOP '***Error: NANG > MXNANG in bhmie' IF(NANG<2)NANG = 2 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 + 4d0 * X ** (1d0 / 3d0) + 2d0 NMX = int(MAX(XSTOP, YMOD)) + 15 ! 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 = int(XSTOP) ! IF(NMX>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 = 0d0 IF(NANG>1)DANG = 0.5d0 * pi / DBLE(NANG - 1) DO J = 1, NANG THETA = DBLE(J - 1) * DANG AMU(J) = COS(THETA) end do DO J = 1, NANG PI0(J) = 0d0 PI1(J) = 1d0 end do NN = 2 * NANG - 1 DO J = 1, NN S1(J) = (0d0, 0d0) S2(J) = (0d0, 0d0) end do ! !*** Logarithmic derivative D(J) calculated by downward recurrence ! beginning with initial value (0.,0.) at J=NMX ! D(NMX) = (0d0, 0d0) NN = NMX - 1 DO N = 1, NN EN = NMX - N + 1 D(NMX - N) = (EN / Y) - (1. / (D(NMX - N + 1) + EN / Y)) end do ! !*** Riccati-Bessel functions with real argument X ! calculated by upward recurrence ! PSI0 = COS(DX) PSI1 = SIN(DX) CHI0 = -SIN(DX) CHI1 = COS(DX) XI1 = DCMPLX(PSI1, -CHI1) QSCA = 0d0 GSCA = 0d0 P = -1d0 DO N = 1, NSTOP EN = N FN = (2d0 * EN + 1d0) / (EN * (EN + 1d0)) ! 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 = (2d0 * EN - 1d0) * PSI1 / DX - PSI0 CHI = (2d0 * EN - 1d0) * CHI1 / DX - CHI0 XI = DCMPLX(PSI, -CHI) ! !*** Store previous values of AN and BN for use ! in computation of g= IF(N>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 + (2d0 * EN + 1d0) * (ABS(AN)**2 + ABS(BN)**2) GSCA = GSCA + ((2d0 * EN + 1d0) / (EN * (EN + 1d0))) * & (dble(AN) * dble(BN) + AIMAG(AN) * AIMAG(BN)) IF(N>1)THEN GSCA = GSCA + ((EN - 1d0) * (EN + 1d0) / EN) * & (dble(AN1) * dble(AN) + AIMAG(AN1) * AIMAG(AN) + & dble(BN1) * dble(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 + 1d0) * PI0(J) S1(J) = S1(J) + FN * (AN * pi_(J) + BN * TAU(J)) S2(J) = S2(J) + FN * (AN * TAU(J) + BN * pi_(J)) end do ! !*** 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)) end do PSI0 = PSI1 PSI1 = PSI CHI0 = CHI1 CHI1 = CHI XI1 = DCMPLX(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. * EN + 1.) * AMU(J) * pi_(J) - (EN + 1.) * PI0(J)) / EN PI0(J) = pi_(J) end do end do ! !*** Have summed sufficient terms. ! Now compute QSCA,QEXT,QBACK,and GSCA GSCA = 2d0 * GSCA / QSCA QSCA = (2d0 / (DX * DX)) * QSCA QEXT = (4d0 / (DX * DX)) * dble(S1(1)) QBACK = (ABS(S1(2 * NANG - 1)) / DX)**2 / pi RETURN END end module bhmie