| [3240] | 1 | module bhmie |
|---|
| 2 | use math, only: pi |
|---|
| 3 | |
|---|
| 4 | implicit none |
|---|
| 5 | |
|---|
| 6 | integer, parameter :: & |
|---|
| 7 | mxnang = 1000, & ! Note: important that MXNANG be consistent with dimension of S1 and S2 in calling routine! |
|---|
| 8 | s_size = 2 * MXNANG - 1 |
|---|
| 9 | |
|---|
| 10 | save |
|---|
| 11 | |
|---|
| 12 | contains |
|---|
| 13 | SUBROUTINE BH_MIE(X, REFREL, NANG, S1, S2, QEXT, QSCA, QBACK, GSCA) |
|---|
| 14 | implicit none |
|---|
| 15 | ! Declare parameters: |
|---|
| 16 | |
|---|
| 17 | INTEGER NMXX |
|---|
| 18 | ! PARAMETER(MXNANG=1000,NMXX=15000) |
|---|
| 19 | PARAMETER(NMXX = 150000) |
|---|
| 20 | ! Arguments: |
|---|
| 21 | INTEGER NANG |
|---|
| 22 | doubleprecision GSCA, QBACK, QEXT, QSCA, X |
|---|
| 23 | COMPLEX(kind=8) REFREL |
|---|
| 24 | COMPLEX(kind=8) S1(s_size), S2(s_size) |
|---|
| 25 | ! Local variables: |
|---|
| 26 | INTEGER J, JJ, N, NSTOP, NMX, NN |
|---|
| 27 | DOUBLE PRECISION CHI, CHI0, CHI1, DANG, DX, EN, FN, P, PSI, PSI0, PSI1, & |
|---|
| 28 | THETA, XSTOP, YMOD |
|---|
| 29 | DOUBLE PRECISION AMU(MXNANG), pi_(MXNANG), PI0(MXNANG), PI1(MXNANG), & |
|---|
| 30 | TAU(MXNANG) |
|---|
| 31 | COMPLEX(kind=8) AN, AN1, BN, BN1, DREFRL, XI, XI1, Y |
|---|
| 32 | COMPLEX(kind=8) D(NMXX) |
|---|
| 33 | !*********************************************************************** |
|---|
| 34 | ! Subroutine BHMIE is the Bohren-Huffman Mie scattering subroutine |
|---|
| 35 | ! to calculate scattering and absorption by a homogenous isotropic |
|---|
| 36 | ! sphere. |
|---|
| 37 | ! Given: |
|---|
| 38 | ! X = 2*pi*a/lambda |
|---|
| 39 | ! REFREL = (complex refr. index of sphere)/(real index of medium) |
|---|
| 40 | ! NANG = number of angles between 0 and 90 degrees |
|---|
| 41 | ! (will calculate 2*NANG-1 directions from 0 to 180 deg.) |
|---|
| 42 | ! if called with NANG<2, will set NANG=2 and will compute |
|---|
| 43 | ! scattering for theta=0,90,180. |
|---|
| 44 | ! Returns: |
|---|
| 45 | ! S1(1 - 2*NANG-1) = -i*f_22 (incid. E perp. to scatt. plane, |
|---|
| 46 | ! scatt. E perp. to scatt. plane) |
|---|
| 47 | ! S2(1 - 2*NANG-1) = -i*f_11 (incid. E parr. to scatt. plane, |
|---|
| 48 | ! scatt. E parr. to scatt. plane) |
|---|
| 49 | ! QEXT = C_ext/pi*a**2 = efficiency factor for extinction |
|---|
| 50 | ! QSCA = C_sca/pi*a**2 = efficiency factor for scattering |
|---|
| 51 | ! QBACK = (dC_sca/domega)/pi*a**2 |
|---|
| 52 | ! = backscattering efficiency [NB: this is (1/4*pi) smaller |
|---|
| 53 | ! than the "radar backscattering efficiency"; see Bohren & |
|---|
| 54 | ! Huffman 1983 pp. 120-123] |
|---|
| 55 | ! GSCA = <cos(theta)> for scattering |
|---|
| 56 | ! |
|---|
| 57 | ! Original program taken from Bohren and Huffman (1983), Appendix A |
|---|
| 58 | ! Modified by B.T.Draine, Princeton Univ. Obs., 90/10/26 |
|---|
| 59 | ! in order to compute <cos(theta)> |
|---|
| 60 | ! 91/05/07 (BTD): Modified to allow NANG=1 |
|---|
| 61 | ! 91/08/15 (BTD): Corrected error (failure to initialize P) |
|---|
| 62 | ! 91/08/15 (BTD): Modified to enhance vectorizability. |
|---|
| 63 | ! 91/08/15 (BTD): Modified to make NANG=2 if called with NANG=1 |
|---|
| 64 | ! 91/08/15 (BTD): Changed definition of QBACK. |
|---|
| 65 | ! 92/01/08 (BTD): Converted to full double precision and double complex |
|---|
| 66 | ! eliminated 2 unneed lines of code |
|---|
| 67 | ! eliminated redundant variables (e.g. APSI,APSI0) |
|---|
| 68 | ! renamed RN -> EN = double precision N |
|---|
| 69 | ! Note that DOUBLE COMPLEX and DCMPLX are not part |
|---|
| 70 | ! of f77 standard, so this version may not be fully |
|---|
| 71 | ! portable. In event that portable version is |
|---|
| 72 | ! needed, use src/bhmie_f77.f |
|---|
| 73 | ! 93/06/01 (BTD): Changed AMAX1 to generic function MAX |
|---|
| 74 | !*********************************************************************** |
|---|
| 75 | !*** Safety checks |
|---|
| 76 | IF(NANG>MXNANG)STOP '***Error: NANG > MXNANG in bhmie' |
|---|
| 77 | IF(NANG<2)NANG = 2 |
|---|
| 78 | |
|---|
| 79 | DX = X |
|---|
| 80 | DREFRL = REFREL |
|---|
| 81 | Y = X * DREFRL |
|---|
| 82 | YMOD = ABS(Y) |
|---|
| 83 | ! |
|---|
| 84 | !*** Series expansion terminated after NSTOP terms |
|---|
| 85 | ! Logarithmic derivatives calculated from NMX on down |
|---|
| 86 | XSTOP = X + 4d0 * X ** (1d0 / 3d0) + 2d0 |
|---|
| 87 | NMX = int(MAX(XSTOP, YMOD)) + 15 |
|---|
| 88 | ! BTD experiment 91/1/15: add one more term to series and compare results |
|---|
| 89 | ! NMX=AMAX1(XSTOP,YMOD)+16 |
|---|
| 90 | ! test: compute 7001 wavelengths between .0001 and 1000 micron |
|---|
| 91 | ! for a=1.0micron SiC grain. When NMX increased by 1, only a single |
|---|
| 92 | ! computed number changed (out of 4*7001) and it only changed by 1/8387 |
|---|
| 93 | ! conclusion: we are indeed retaining enough terms in series! |
|---|
| 94 | NSTOP = int(XSTOP) |
|---|
| 95 | ! |
|---|
| 96 | IF(NMX>NMXX)THEN |
|---|
| 97 | WRITE(0, *)'Error: NMX > NMXX=', NMXX, ' for |m|x=', YMOD |
|---|
| 98 | STOP |
|---|
| 99 | ENDIF |
|---|
| 100 | !*** Require NANG.GE.1 in order to calculate scattering intensities |
|---|
| 101 | DANG = 0d0 |
|---|
| 102 | IF(NANG>1)DANG = 0.5d0 * pi / DBLE(NANG - 1) |
|---|
| 103 | DO J = 1, NANG |
|---|
| 104 | THETA = DBLE(J - 1) * DANG |
|---|
| 105 | AMU(J) = COS(THETA) |
|---|
| 106 | end do |
|---|
| 107 | DO J = 1, NANG |
|---|
| 108 | PI0(J) = 0d0 |
|---|
| 109 | PI1(J) = 1d0 |
|---|
| 110 | end do |
|---|
| 111 | NN = 2 * NANG - 1 |
|---|
| 112 | DO J = 1, NN |
|---|
| 113 | S1(J) = (0d0, 0d0) |
|---|
| 114 | S2(J) = (0d0, 0d0) |
|---|
| 115 | end do |
|---|
| 116 | ! |
|---|
| 117 | !*** Logarithmic derivative D(J) calculated by downward recurrence |
|---|
| 118 | ! beginning with initial value (0.,0.) at J=NMX |
|---|
| 119 | ! |
|---|
| 120 | D(NMX) = (0d0, 0d0) |
|---|
| 121 | NN = NMX - 1 |
|---|
| 122 | DO N = 1, NN |
|---|
| 123 | EN = NMX - N + 1 |
|---|
| 124 | D(NMX - N) = (EN / Y) - (1. / (D(NMX - N + 1) + EN / Y)) |
|---|
| 125 | end do |
|---|
| 126 | ! |
|---|
| 127 | !*** Riccati-Bessel functions with real argument X |
|---|
| 128 | ! calculated by upward recurrence |
|---|
| 129 | ! |
|---|
| 130 | PSI0 = COS(DX) |
|---|
| 131 | PSI1 = SIN(DX) |
|---|
| 132 | CHI0 = -SIN(DX) |
|---|
| 133 | CHI1 = COS(DX) |
|---|
| 134 | XI1 = DCMPLX(PSI1, -CHI1) |
|---|
| 135 | QSCA = 0d0 |
|---|
| 136 | GSCA = 0d0 |
|---|
| 137 | P = -1d0 |
|---|
| 138 | DO N = 1, NSTOP |
|---|
| 139 | EN = N |
|---|
| 140 | FN = (2d0 * EN + 1d0) / (EN * (EN + 1d0)) |
|---|
| 141 | ! for given N, PSI = psi_n CHI = chi_n |
|---|
| 142 | ! PSI1 = psi_{n-1} CHI1 = chi_{n-1} |
|---|
| 143 | ! PSI0 = psi_{n-2} CHI0 = chi_{n-2} |
|---|
| 144 | ! Calculate psi_n and chi_n |
|---|
| 145 | PSI = (2d0 * EN - 1d0) * PSI1 / DX - PSI0 |
|---|
| 146 | CHI = (2d0 * EN - 1d0) * CHI1 / DX - CHI0 |
|---|
| 147 | XI = DCMPLX(PSI, -CHI) |
|---|
| 148 | ! |
|---|
| 149 | !*** Store previous values of AN and BN for use |
|---|
| 150 | ! in computation of g=<cos(theta)> |
|---|
| 151 | IF(N>1)THEN |
|---|
| 152 | AN1 = AN |
|---|
| 153 | BN1 = BN |
|---|
| 154 | ENDIF |
|---|
| 155 | ! |
|---|
| 156 | !*** Compute AN and BN: |
|---|
| 157 | AN = (D(N) / DREFRL + EN / DX) * PSI - PSI1 |
|---|
| 158 | AN = AN / ((D(N) / DREFRL + EN / DX) * XI - XI1) |
|---|
| 159 | BN = (DREFRL * D(N) + EN / DX) * PSI - PSI1 |
|---|
| 160 | BN = BN / ((DREFRL * D(N) + EN / DX) * XI - XI1) |
|---|
| 161 | ! |
|---|
| 162 | !*** Augment sums for Qsca and g=<cos(theta)> |
|---|
| 163 | QSCA = QSCA + (2d0 * EN + 1d0) * (ABS(AN)**2 + ABS(BN)**2) |
|---|
| 164 | GSCA = GSCA + ((2d0 * EN + 1d0) / (EN * (EN + 1d0))) * & |
|---|
| 165 | (dble(AN) * dble(BN) + AIMAG(AN) * AIMAG(BN)) |
|---|
| 166 | IF(N>1)THEN |
|---|
| 167 | GSCA = GSCA + ((EN - 1d0) * (EN + 1d0) / EN) * & |
|---|
| 168 | (dble(AN1) * dble(AN) + AIMAG(AN1) * AIMAG(AN) + & |
|---|
| 169 | dble(BN1) * dble(BN) + AIMAG(BN1) * AIMAG(BN)) |
|---|
| 170 | ENDIF |
|---|
| 171 | ! |
|---|
| 172 | !*** Now calculate scattering intensity pattern |
|---|
| 173 | ! First do angles from 0 to 90 |
|---|
| 174 | DO J = 1, NANG |
|---|
| 175 | JJ = 2 * NANG - J |
|---|
| 176 | pi_(J) = PI1(J) |
|---|
| 177 | TAU(J) = EN * AMU(J) * pi_(J) - (EN + 1d0) * PI0(J) |
|---|
| 178 | S1(J) = S1(J) + FN * (AN * pi_(J) + BN * TAU(J)) |
|---|
| 179 | S2(J) = S2(J) + FN * (AN * TAU(J) + BN * pi_(J)) |
|---|
| 180 | end do |
|---|
| 181 | ! |
|---|
| 182 | !*** Now do angles greater than 90 using pi_ and TAU from |
|---|
| 183 | ! angles less than 90. |
|---|
| 184 | ! P=1 for N=1,3,...; P=-1 for N=2,4,... |
|---|
| 185 | P = -P |
|---|
| 186 | DO J = 1, NANG - 1 |
|---|
| 187 | JJ = 2 * NANG - J |
|---|
| 188 | S1(JJ) = S1(JJ) + FN * P * (AN * pi_(J) - BN * TAU(J)) |
|---|
| 189 | S2(JJ) = S2(JJ) + FN * P * (BN * pi_(J) - AN * TAU(J)) |
|---|
| 190 | end do |
|---|
| 191 | PSI0 = PSI1 |
|---|
| 192 | PSI1 = PSI |
|---|
| 193 | CHI0 = CHI1 |
|---|
| 194 | CHI1 = CHI |
|---|
| 195 | XI1 = DCMPLX(PSI1, -CHI1) |
|---|
| 196 | ! |
|---|
| 197 | !*** Compute pi_n for next value of n |
|---|
| 198 | ! For each angle J, compute pi_n+1 |
|---|
| 199 | ! from pi_ = pi_n , PI0 = pi_n-1 |
|---|
| 200 | DO J = 1, NANG |
|---|
| 201 | PI1(J) = ((2. * EN + 1.) * AMU(J) * pi_(J) - (EN + 1.) * PI0(J)) / EN |
|---|
| 202 | PI0(J) = pi_(J) |
|---|
| 203 | end do |
|---|
| 204 | end do |
|---|
| 205 | ! |
|---|
| 206 | !*** Have summed sufficient terms. |
|---|
| 207 | ! Now compute QSCA,QEXT,QBACK,and GSCA |
|---|
| 208 | GSCA = 2d0 * GSCA / QSCA |
|---|
| 209 | QSCA = (2d0 / (DX * DX)) * QSCA |
|---|
| 210 | QEXT = (4d0 / (DX * DX)) * dble(S1(1)) |
|---|
| 211 | QBACK = (ABS(S1(2 * NANG - 1)) / DX)**2 / pi |
|---|
| 212 | RETURN |
|---|
| 213 | END |
|---|
| 214 | end module bhmie |
|---|
