Line data Source code
1 : !! See Bohren and Huffman, "Absorption and Scattering of Light by
2 : !! Small Particles", 1983, pg 480 (in Appendix A).
3 : !!
4 : !! Subroutine bhmie calculates amplitude scattering matrix
5 : !! elements and efficiencies for extinction, total scattering
6 : !! and backscattering for a given size parameter and
7 : !! relative refractive index.
8 : !!
9 : !! From the main program:
10 : !! refrel = cmplx(refre,refim) / refmed
11 : !!
12 : !! @author Chuck Bardeen
13 : !! @version 2011
14 0 : subroutine bhmie(carma, x, refrel, nang, s1, s2, Qext, Qsca, Qback, gfac, rc)
15 :
16 : ! types
17 : use carma_precision_mod
18 : use carma_enums_mod, only : RC_ERROR
19 : use carma_types_mod, only : carma_type
20 : use carma_mod
21 :
22 : implicit none
23 :
24 : type(carma_type), intent(in) :: carma !! the carma object
25 : real(kind=f), intent(in) :: x !! radius / wavelength
26 : complex(kind=f), intent(in) :: refrel !! refractive index particle / reference refractive index
27 : integer, intent(in) :: nang !! number of angles in s1 and s2
28 : complex(kind=f), intent(out) :: s1(2*nang-1) !! CORE RADIUS
29 : complex(kind=f), intent(out) :: s2(2*nang-1) !! REAL PART OF THE CORE INDEX OF REFRACTION
30 : real(kind=f), intent(out) :: Qext !! EFFICIENCY FACTOR FOR EXTINCTION
31 : real(kind=f), intent(out) :: Qsca !! EFFICIENCY FACTOR FOR SCATTERING
32 : real(kind=f), intent(out) :: Qback !! BACK SCATTER CROSS SECTION.
33 : real(kind=f), intent(out) :: gfac !! asymmetry factor
34 : integer, intent(inout) :: rc !! return code, negative indicates failure
35 :
36 :
37 : real(kind=f) :: amu(100), theta(100), pi(100), tau(100), pi0(100), pi1(100)
38 : complex(kind=f) :: y, xi, xi0, xi1, an, bn
39 0 : complex(kind=f), allocatable :: d(:)
40 : complex(kind=f) :: ccan, ccbn, anmi1, bnmi1
41 : real(kind=f) :: psi0, psi1, psi, dn, dx, chi0, chi1, apsi0, apsi1, g1, g2
42 : real(kind=f) :: dang, fn, ffn, apsi, chi, p, t, xstop, ymod
43 : integer :: j, jj, n, nn, rn, nmx, nstop
44 :
45 :
46 : ! Mie x and y values.
47 0 : dx = x
48 0 : y = x * refrel
49 :
50 : ! Series terminated after nstop terms
51 0 : xstop = x + 4._f * x**0.3333_f + 2.0_f
52 0 : nstop = xstop
53 :
54 : ! Will loop over nang angles.
55 0 : ymod = int(abs(y))
56 0 : nmx = max(xstop, ymod) + 15
57 0 : dang = 1.570796327_f / real(nang - 1, kind=f)
58 0 : allocate(d(nmx))
59 :
60 0 : do j = 1, nang
61 0 : theta(j) = (real(j, kind=f) - 1._f) * dang
62 0 : amu(j) = cos(theta(j))
63 : end do
64 :
65 : ! Logarithmic derivative d(j) calculated by downword
66 : ! recurrence beginning with initial value 0.0 + i*0.0
67 : ! at j = nmx
68 0 : d(nmx) = cmplx(0.0_f, 0.0_f, kind=f)
69 0 : nn = nmx-1
70 : ! write(*,*) 'nmx=',nmx,' d(nmx)=',d(nmx), ' nn=',nn
71 :
72 0 : do n = 1, nn
73 0 : rn = nmx - n + 1
74 0 : d(nmx-n) = (rn/y) - (1._f / (d(nmx - n + 1) + rn / y))
75 :
76 : ! write(*,*) 'n=',n,' rn=',rn,' y=', y,' d(nmx-n)=',d(nmx-n)
77 : ! write(*,*) 'rn/y=',rn/y, 'd(nmx-n+1)=',d(nmx-n+1),'(d(nmx-n+1)+rn/y)', &
78 : ! (d(nmx-n+1)+rn/y),'1./(d(nmx-n+1)+rn/y)',1./(d(nmx-n+1)+rn/y)
79 : end do
80 :
81 0 : pi0(1:nang) = 0.0_f
82 0 : pi1(1:nang) = 1.0_f
83 :
84 0 : nn = 2 * nang-1
85 0 : s1(1:nn) = cmplx(0.0_f, 0.0_f, kind=f)
86 0 : s2(1:nn) = cmplx(0.0_f, 0.0_f, kind=f)
87 :
88 : ! Riccati-Bessel functions with real argument x
89 : ! calculated by upward recurrence
90 0 : psi0 = cos(dx)
91 0 : psi1 = sin(dx)
92 0 : chi0 = -sin(x)
93 0 : chi1 = cos(x)
94 0 : apsi0 = psi0
95 0 : apsi1 = psi1
96 0 : xi0 = cmplx(apsi0,-chi0, kind=f)
97 0 : xi1 = cmplx(apsi1,-chi1, kind=f)
98 0 : Qsca = 0.0_f
99 0 : g1 = 0.0_f
100 0 : g2 = 0.0_f
101 0 : n = 1
102 :
103 : ! Loop over the terms n in the Mie series
104 : do while (.true.)
105 0 : dn = n
106 0 : rn = n
107 0 : fn = (2._f * rn + 1._f) / (rn * (rn + 1._f))
108 0 : ffn = (rn - 1._f) * (rn + 1._f) / rn
109 0 : psi = (2._f * dn - 1._f) * psi1 / dx - psi0
110 0 : apsi = psi
111 0 : chi = (2._f * rn - 1._f) * chi1 / x - chi0
112 0 : xi = cmplx(apsi, -chi, kind=f)
113 : ! write(*,*) 'n=', n
114 : ! write(*,*) 'd(n)=',d(n),' refrel=',refrel,' rn=',rn, ' x=',x,'apsi=',apsi,' apsi1=',apsi1
115 :
116 0 : an = (d(n) / refrel + rn / x) * apsi - apsi1
117 : ! write(*,*) 'an=',an,' xi=',xi,' xi1=',xi1
118 :
119 0 : an = an / ((d(n) / refrel + rn / x) * xi - xi1)
120 0 : bn = (refrel * d(n) + rn / x) * apsi - apsi1
121 0 : bn = bn / ((refrel * d(n) + rn / x) * xi - xi1)
122 0 : ccan = conjg(an)
123 0 : ccbn = conjg(bn)
124 0 : g2 = g2 + fn * real(an * ccbn)
125 :
126 0 : if (n-1 > 0) then
127 0 : g1 = g1 + ffn * real(anmi1 * ccan + bnmi1 * ccbn)
128 : end if
129 0 : Qsca = Qsca + (2._f * rn + 1._f) * (abs(an) * abs(an) + abs(bn) * abs(bn))
130 :
131 0 : do j = 1, nang
132 0 : jj = 2 * nang-j
133 0 : pi(j) = pi1(j)
134 0 : tau(j) = rn * amu(j) * pi(j) - (rn + 1._f) * pi0(j)
135 0 : p = (-1._f)**(n-1)
136 : ! write(*,*) 'fn=',fn,' an=',an,' bn=',bn,' pi(j)=',pi(j),' tau(j)=',tau(j)
137 :
138 0 : s1(j) = s1(j) + fn * (an * pi(j) + bn * tau(j))
139 0 : t = (-1._f)**n
140 0 : s2(j) = s2(j) + fn * (an * tau(j) + bn * pi(j))
141 :
142 0 : if (j.ne.jj) then
143 0 : s1(jj)=s1(jj) + fn*(an*pi(j)*p+bn*tau(j)*t)
144 0 : s2(jj)=s2(jj) + fn*(an*tau(j)*t+bn*pi(j)*p)
145 : ! write(*,*) 'j=',j,' s1(j)=',s1(j),' s2(j)=',s2(j)
146 : end if
147 : end do
148 :
149 0 : psi0 = psi1
150 0 : psi1 = psi
151 0 : apsi1 = psi1
152 0 : chi0 = chi1
153 0 : chi1 = chi
154 0 : xi1 = cmplx(apsi1, -chi1, kind=f)
155 0 : n = n+1
156 0 : rn = n
157 :
158 0 : do j = 1, nang
159 0 : pi1(j) = ((2._f * rn - 1._f) / (rn - 1._f)) * amu(j) * pi(j)
160 0 : pi1(j) = pi1(j) -rn * pi0(j) / (rn - 1._f)
161 0 : pi0(j) = pi(j)
162 : end do
163 :
164 0 : anmi1 = an
165 0 : bnmi1 = bn
166 :
167 0 : if (n - 1 - nstop >= 0) exit
168 :
169 : end do
170 :
171 0 : Qsca = (2._f / (x * x)) * Qsca
172 0 : gfac = (4._f / (x * x * Qsca)) * (g1+g2)
173 0 : Qext = (4._f / (x * x)) * real(s1(1))
174 0 : Qback = (4._f / (x * x)) * abs(s1(2 * nang - 1)) * abs(s1(2 * nang - 1))
175 :
176 : ! write(*,*) 'x',x,' s1(1)=',s1(1),' real(s1(1))=',real(s1(1))
177 : ! write(*,*) 'Qsca=',Qsca,' gfac=',gfac,' Qext=',Qext,'Qback=',Qback
178 :
179 0 : deallocate(d)
180 :
181 0 : return
182 0 : end subroutine bhmie
|
