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