-
Notifications
You must be signed in to change notification settings - Fork 0
/
cylinder.jl
105 lines (88 loc) · 2.78 KB
/
cylinder.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
@doc raw"""
A cylindrical scatterer.
Attributes:
- `r`: the radius of the cylinder base.
- `h`: the height of the cylinder.
- `m`: the relative complex refractive index.
"""
struct Cylinder{T, CT} <: AbstractAxisymmetricShape{T, CT}
r::T
h::T
m::CT
end
volume(c::Cylinder) = π * c.r^2 * c.h
volume_equivalent_radius(c::Cylinder) = ∛(3 * c.r^2 * c.h / 4)
has_symmetric_plane(::Cylinder) = true
@testitem "Utility functions are correct" begin
using TransitionMatrices: Cylinder, volume, volume_equivalent_radius,
has_symmetric_plane
@testset "r = $r, h = $h" for (r, h, V, rᵥ) in [
(2.0, 0.5, 6.283185307179586, 1.1447142425533319),
(1.0, 1.0, 3.141592653589793, 0.9085602964160698),
(0.5, 2.0, 1.5707963267948966, 0.7211247851537042),
]
c = Cylinder(r, h, 1.311)
@test volume(c) ≈ V
@test volume_equivalent_radius(c) ≈ rᵥ
@test has_symmetric_plane(c)
end
end
@doc raw"""
```
gaussquad(c::Cylinder{T}, ngauss) where {T}
```
Evaluate the quadrature points, weights and the corresponding radius and radius derivative (to ``\vartheta``) for a cylinder.
Returns: (`x`, `w`, `r`, `r′`)
- `x`: the quadrature points.
- `w`: the quadrature weights.
- `r`: the radius at each quadrature point.
- `r′`: the radius derivative at each quadrature point.
"""
function gaussquad(c::Cylinder{T}, ngauss) where {T}
ng = ngauss ÷ 2
ng1 = ng ÷ 2
ng2 = ng - ng1
x = zeros(T, ngauss)
w = zeros(T, ngauss)
x1, w1 = gausslegendre(T, ng1)
x2, w2 = gausslegendre(T, ng2)
xx = -cos(atan(2c.r / c.h))
@. x[1:ng1] = 0.5(xx + 1.0) * x1 + 0.5(xx - 1.0)
@. w[1:ng1] = 0.5(xx + 1.0) * w1
@. x[(ng1 + 1):ng] = -0.5xx * x2 + 0.5xx
@. w[(ng1 + 1):ng] = -0.5xx * w2
@. x[(ng + 1):ngauss] = (-1.0) * x[ng:-1:1]
@. w[(ng + 1):ngauss] = w[ng:-1:1]
h = c.h / 2
d = c.r
r = similar(x)
r′ = similar(x)
# We could have used @turbo for primitive types like
# `Float64`, but this is not the bottleneck, so we only use
# `@simd` here.
@simd for i in 1:(ngauss ÷ 2)
cosϑ = abs(x[i])
sinϑ = √(1 - cosϑ^2)
if h / cosϑ < d / sinϑ
r[i] = h / cosϑ
r′[i] = h * sinϑ / cosϑ^2
else
r[i] = d / sinϑ
r′[i] = -d * cosϑ / sinϑ^2
end
r[ngauss + 1 - i] = r[i]
r′[ngauss + 1 - i] = r′[i]
r′[i] = -r′[i]
end
return x, w, r, r′
end
function rmax(c::Cylinder)
return hypot(c.r, c.h / 2)
end
function rmin(c::Cylinder)
return min(c.r, c.h / 2)
end
function Base.:∈(x, c::Cylinder)
return abs2(x[1]) + abs2(x[2]) <= c.r^2 && abs(x[3]) <= c.h / 2
end
refractive_index(c::Cylinder, x) = x ∈ c ? c.m : one(c.m)