-
Notifications
You must be signed in to change notification settings - Fork 0
/
bhcoat.jl
149 lines (136 loc) · 4.88 KB
/
bhcoat.jl
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
@doc raw"""
```
bhcoat([T=Float64,], xᵢₙ, xₒᵤₜ, mᵢₙ, mₒᵤₜ; nₘₐₓ, tolerance = 1e-8)
```
Inputs:
- `T`: Type used for calculation. All real numbers will be stored as `T`, while complex numbers will be stored as `C = complex(T)`.
- `xᵢₙ`: Size parameter of the inner sphere. Defined as ``\frac{2\pi r}{\lambda}``
- `xₒᵤₜ`: Size parameter of the coated sphere. `xₒᵤₜ >= xᵢₙ` should hold.
- `mᵢₙ`: Refractive index of the inner sphere, relative to the host medium.
- `mₒᵤₜ`: Refractive index of the mantle, relative to the host medium.
Keyword arguments:
- `nₘₐₓ`: Maximum order of the Mie coefficients. Default to ``\max(x_m + 4\sqrt{3}{x_m} + 2, \max(|m_c|, |m_m|)x_m)``.
- `tolerance`: Error tolerance. Default to `1e-8`.
Outputs:
- `a, b`: Mie coefficients. Both are of type `Vector{C}`.
References:
- Bohren, C.F., Huffman, D.R., 1983. Absorption and scattering of light by small particles. John Wiley & Sons.
"""
function bhcoat(T, xᵢₙ,
xₒᵤₜ, mᵢₙ, mₒᵤₜ;
nₘₐₓ = ceil(Int,
max(xₒᵤₜ + 4 * ∛xₒᵤₜ + 2,
xₒᵤₜ * max(abs(mᵢₙ), abs(mₒᵤₜ)))),
tolerance = 1e-8)
@assert xₒᵤₜ>=xᵢₙ "xₒᵤₜ must be greater than or equal to xᵢₙ"
C = complex(T)
x = T(xᵢₙ)
y = T(xₒᵤₜ)
m₁ = C(mᵢₙ)
m₂ = C(mₒᵤₜ)
x₁ = m₁ * x
x₂ = m₂ * x
y₂ = m₂ * y
m = m₂ / m₁
a = zeros(C, nₘₐₓ)
b = zeros(C, nₘₐₓ)
d0x1 = cot(x₁)
d0x2 = cot(x₂)
d0y2 = cot(y₂)
d1x1 = zero(T)
d1x2 = zero(T)
brack = zero(T)
crack = zero(T)
psi0y = cos(y)
psi1y = sin(y)
chi0y = -sin(y)
chi1y = cos(y)
xi0y = psi0y - 1.0im * chi0y
xi1y = psi1y - 1.0im * chi1y
chi0y2 = -sin(y₂)
chi1y2 = cos(y₂)
chipy2 = zero(T)
chiy2 = zero(T)
chi0x2 = -sin(x₂)
chi1x2 = cos(x₂)
chipx2 = zero(T)
chix2 = zero(T)
amess1 = zero(T)
amess2 = zero(T)
amess3 = zero(T)
amess4 = zero(T)
iflag = false
for n in 1:nₘₐₓ
rn = T(n)
psiy = (2rn - 1) * psi1y / y - psi0y
chiy = (2rn - 1) * chi1y / y - chi0y
xiy = psiy - 1.0im * chiy
d1y2 = 1 / (rn / y₂ - d0y2) - rn / y₂
if !iflag
d1x1 = 1 / (rn / x₁ - d0x1) - rn / x₁
d1x2 = 1 / (rn / x₂ - d0x2) - rn / x₂
chix2 = (2rn - 1) * chi1x2 / x₂ - chi0x2
chiy2 = (2rn - 1) * chi1y2 / y₂ - chi0y2
chipx2 = chi1x2 - rn * chix2 / x₂
chipy2 = chi1y2 - rn * chiy2 / y₂
ancap = m * d1x1 - d1x2
ancap = ancap / (m * d1x1 * chix2 - chipx2)
ancap = ancap / (chix2 * d1x2 - chipx2)
brack = ancap * (chiy2 * d1y2 - chipy2)
bncap = m * d1x2 - d1x1
bncap = bncap / (m * chipx2 - d1x1 * chix2)
bncap = bncap / (chix2 * d1x2 - chipx2)
crack = bncap * (chiy2 * d1y2 - chipy2)
amess1 = brack * chipy2
amess2 = brack * chiy2
amess3 = crack * chipy2
amess4 = crack * chiy2
end
if abs(amess1) < tolerance * abs(d1y2) && abs(amess2) < tolerance &&
abs(amess3) < tolerance * abs(d1y2) && abs(amess4) < tolerance
brack = zero(T)
crack = zero(T)
iflag = true
else
iflag = false
end
dnbar = d1y2 - brack * chipy2
dnbar = dnbar / (1.0 - brack * chiy2)
gnbar = d1y2 - crack * chipy2
gnbar = gnbar / (1.0 - crack * chiy2)
a[n] = ((dnbar / m₂ + rn / y) * psiy - psi1y) / ((dnbar / m₂ + rn / y) * xiy - xi1y)
b[n] = ((m₂ * gnbar + rn / y) * psiy - psi1y) / ((m₂ * gnbar + rn / y) * xiy - xi1y)
psi0y = psi1y
psi1y = psiy
chi0y = chi1y
chi1y = chiy
xi1y = psi1y - 1.0im * chi1y
chi0x2 = chi1x2
chi1x2 = chix2
chi0y2 = chi1y2
chi1y2 = chiy2
d0x1 = d1x1
d0x2 = d1x2
d0y2 = d1y2
end
return a, b
end
function bhcoat(xᵢₙ, xₒᵤₜ, mᵢₙ, mₒᵤₜ;
nₘₐₓ = ceil(Int,
max(xₒᵤₜ + 4 * ∛xₒᵤₜ + 2,
xₒᵤₜ * max(abs(mᵢₙ), abs(mₒᵤₜ)))),
tolerance = 1e-8)
bhcoat(Float64, xᵢₙ, xₒᵤₜ, mᵢₙ, mₒᵤₜ; nₘₐₓ = nₘₐₓ, tolerance = tolerance)
end
@testitem "bhcoat" begin
using TransitionMatrices: bhmie, bhcoat
@testset "converges to bhmie when x = $x, m = $m" for (x, m) in [
(1.0, 1.311),
(2.0, 1.5 + 0.01im),
]
am, bm = bhmie(x, m)
ac, bc = bhcoat(x, x, m, m)
@test all(isapprox.(am, ac; atol = 1e-12))
@test all(isapprox.(bm, bc; atol = 1e-12))
end
end