-
Notifications
You must be signed in to change notification settings - Fork 0
/
BoundaryInteractions.py
200 lines (145 loc) · 6.05 KB
/
BoundaryInteractions.py
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
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
"""
Class simulating boundary interactions. Assumes particles are
already at boundary.
Author: Jyotirmai (Joe) Singh 10/7/18
"""
from UtilityMethods import *
import AnharmonicDecay
def specular_scatter(particle, box):
"""
Simulates specular scattering.
:param particle: The particle undergoing the interaction
:param box: The box in which the particle is in
"""
print("SPECULAR SCATTER")
x, y, z = particle.get_x(), particle.get_y(), particle.get_z()
vx, vy, vz = particle.get_vx(), particle.get_vy(), particle.get_vz()
if x <= 0 or x >= box.width:
vx = -vx
elif y <= 0 or y >= box.height:
vy = -vy
elif z <= 0 or z >= box.depth:
vz = -vz
particle.set_velocity(vx, vy, vz)
def get_lambertian_probability(particle):
# Convert frequency to GHz
f = particle.get_f()/1e9
return -2.98e-11 * (f**4) + 1.71e-8 * (f**3) - 2.47e-6 * (f**2) + 7.83e-4 * f + 5.88e-2
def get_sad_probability(particle):
# Convert frequency to GHz
f = particle.get_f()/1e9
return 1.51e-14 * (f**5)
def get_specular_probability(particle):
# Convert frequency to GHz
f = particle.get_f()/1e9
return 2.9e-13 * (f**4) + 3.1e-9 * (f**3) - 3.21e-6 * (f**2) - 2.03e-4 * f + 0.928
def lambertian_scatter(particle, box):
"""
Simulates scattering by a Lambertian process with
a polar angle dependent on a cos(theta) distribution.
:param particle: The particle undergoing the interaction
:param box: The box in which the particle is in
"""
diffusive_angle = get_cos_angle()
azimuthal_angle = np.random.uniform(0, 2 * PI)
print("DIFFUSIVE SCATTER")
adjust_boundary_velocity(particle, box, diffusive_angle, azimuthal_angle)
def cylindrical_vol_and_sa(r, h):
"""
Calculate volume and surface area of a cylinder.
:param r: The radius in m
:param h: The height in m
:return: The volume and surface area
"""
return PI * h * r ** 2, 2 * PI * (r * h + r ** 2)
def convert_w_to_T(w):
"""
Converts from angular frequency to temperature as in Eq (4)
in Klistner and Pohl.
:param w: Angular frequency, 2 * PI * v_dom where v_dom is the dominant frequency
Klistner and Pohl mention, Eq. 4 with the factor of 4.25
:return: The temperature.
"""
return (h * w) / (2 * PI * k_b * 4.25)
def diffusive_scatter_rate(T, particle, box):
"""
The total diffusive scatter rate (Lambertian + Surface Anharmonic) as
defined by Fig. 1 in Klistner and Pohl. The equation here is obtained by
extracting the points from their plot and then fitting via matplotlib.
Multiply by 100 to correct for cm vs m.
:param T: Temperature
:param particle: The particle
:param box: The box
:return: The decay rate.
"""
material = box.get_material()
phonon_velocity = material.get_particle_velocity(particle.get_type())
return 100 * (0.036452 * (T ** 2) + 0.216544) * phonon_velocity
def boundary_interaction(particle, box, points, colours):
"""
Method to simulate boundary interactions. Assumes the particle
is already present at the boundary and then simulates Lambertian
or SAD scattering as necessary.
:param particle: The particle undergoing the interactions
:param box: The box in which the particle is in
:param points: The array of points associated with the phonon locations
:param colours: The colour configuration of the phonons
"""
# First check if we are going to potentially be absorbed by the phonon detectors.
if np.random.rand() < box.get_coverage():
# Particle incident on aluminium
if np.random.rand() < box.get_material().get_sensor_absorb_probability():
# Particle absorbed by aluminium
remove_particle(particle, box, points)
return
# If no absorption, continue the process.
material = box.get_material()
lambertian_prob = get_lambertian_probability(particle)
sad_prob = get_sad_probability(particle)
specular_prob = get_specular_probability(particle)
# Rescale probabilities to account for minor errors that mean
# probabilities do not sum exactly to 1
total_prob = lambertian_prob + sad_prob + specular_prob
lambertian_prob /= total_prob
sad_prob /= total_prob
specular_prob /= total_prob
if particle.get_type() != 3:
process = (np.random.choice(2, 1, p=[1-specular_prob, specular_prob]) + 1)[0]
else:
process = (np.random.choice(3, 1, p=[lambertian_prob, specular_prob, sad_prob]) + 1)[0]
if process == 1:
lambertian_scatter(particle, box)
elif process == 2:
specular_scatter(particle, box)
else:
assert process == 3
LTT_prob = material.get_LTT_ratio()
if np.random.rand() < LTT_prob:
# Do LTT boundary decay
print("SURFACE ANHARMONIC LTT")
AnharmonicDecay.anharmonic_decay_LTT(particle, box, 0, points, colours, 1)
else:
# Do LLT boundary decay
print("SURFACE ANHARMONIC LLT")
AnharmonicDecay.anharmonic_decay_LLT(particle, box, 0, points, colours, 1)
def propagate(particle, box, t, points, colours):
"""
Simulate particle moving forward for time t, taking hitting
wall into account.
:param particle: Particle moving.
:param box: Box in which particle is in.
:param t: Time for which box is simulated to move.
"""
t_boundary, x_boundary, y_boundary, z_boundary = time_to_boundary(box, particle)
# If particle propagates for time t and goes beyond the boundary, stop it at the
# boundary and change its velocity to simulate bouncing off the wall
if t_boundary <= t:
particle.set_x(x_boundary)
particle.set_y(y_boundary)
particle.set_z(z_boundary)
# Simulates specular/diffusive scattering.
boundary_interaction(particle, box, points, colours)
# Otherwise can safely propagate for time t without hitting the wall.
else:
print("In propagation non boundary case.")
particle.advance(t)