-
Notifications
You must be signed in to change notification settings - Fork 2
/
photon.m
164 lines (163 loc) · 5.11 KB
/
photon.m
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
classdef photon
properties
x, y, z
ux, uy, uz
w
dead
s
scatters
layer
end
methods
function obj = initialize(obj, x, y, ly_ls)
% initialize photon
obj.x = x;
obj.y = y;
obj.z = 0;
obj.ux = 0;
obj.uy = 0;
obj.uz = 1;
obj.dead = false;
obj.s = 0;
obj.scatters = 0;
obj.layer = 2;
% in case that first layer is clear, do reflection
amb_ly = ly_ls(1);
this_ly = ly_ls(2);
bot_ly = ly_ls(3);
n1 = amb_ly.n;
n2 = this_ly.n;
n3 = bot_ly.n;
if this_ly.clear
r1 = ((n1 - n2)/(n1 + n2))^2;
r2 = ((n3 - n2)/(n3 + n2))^2;
r_sp = r1 + (1-r1)^2*r2/(1-r1*r2);
obj.w = 1 - r_sp;
else
if n1 ~= n2
r_sp = ((n1 - n2)/(n1 + n2))^2;
obj.w = 1 - r_sp;
else
obj.w = 1;
end
end
end
function obj = get_s(obj)
% get random s
xi = rand();
obj.s = -log(xi);
end
function [obj, delta_w] = absorb(obj, ly_ls)
% deduce weight due to absorption
this_ly = ly_ls(obj.layer);
mu_ratio = this_ly.mu_a/this_ly.mu_t;
delta_w = obj.w * mu_ratio;
obj.w = obj.w - delta_w;
end
function obj = scatter(obj, ly_ls)
% update direction cosines upon scattering
xi = rand();
this_ly = ly_ls(obj.layer);
g = this_ly.g;
if g ~= 0
cosine = 1 + g^2 - ((1-g^2)/(1-g+2*g*xi))^2;
cosine = cosine / (2*g);
else
cosine = 2*xi - 1;
end
theta = acos(cosine);
xi = rand();
phi = 2 * pi * xi;
mux = obj.ux;
muy = obj.uy;
muz = obj.uz;
if abs(muz) < 0.99999
obj.ux = sin(theta)*(mux*muz*cos(phi)-muy*sin(phi));
obj.ux = obj.ux/sqrt(1-muz^2) + mux*cosine;
obj.uy = sin(theta)*(muy*muz*cos(phi)+mux*sin(phi));
obj.uy = obj.uy/sqrt(1-muz^2) + muy*cosine;
obj.uz = -sin(theta)*cos(phi)*sqrt(1-muz^2) + muz*cosine;
else
obj.ux = sin(theta) * cos(phi);
obj.uy = sin(theta) * sin(phi);
obj.uz = sign(muz) * cosine;
end
obj.scatters = obj.scatters + 1;
end
function obj = move(obj, ly_ls)
% judge if the photon will hit boundary then move it
this_ly = ly_ls(obj.layer);
mu_t = this_ly.mu_t;
z_bound = [this_ly.z0, this_ly.z1];
if obj.uz < 0
db = (z_bound(1)-obj.z) / obj.uz;
elseif obj.uz == 0
db = inf;
else
db = (z_bound(2)-obj.z) / obj.uz;
end
if db*mu_t <= obj.s
% hit boundary
obj.x = obj.x + obj.ux*db;
obj.y = obj.y + obj.uy*db;
obj.z = obj.z + obj.uz*db;
obj.s = obj.s - db*mu_t;
else
% does not hit boundary
obj.x = obj.x + obj.ux*obj.s/mu_t;
obj.y = obj.y + obj.uy*obj.s/mu_t;
obj.z = obj.z + obj.uz*obj.s/mu_t;
obj.s = 0;
end
end
function obj = reflect_transmit(obj, ly_ls)
% do reflection or transmission
this_ly = ly_ls(obj.layer);
a_i = acos(abs(obj.uz));
n_i = this_ly.n;
if obj.uz < 0
next_ly = ly_ls(obj.layer-1);
dir = -1;
else
next_ly = ly_ls(obj.layer+1);
dir = 1;
end
n_t = next_ly.n;
a_t = asin(n_i*sin(a_i)/n_t);
if n_i > n_t && a_i > asin(n_t/n_i)
r = 1;
else
r = (sin(a_i-a_t)/sin(a_i+a_t))^2 ;
r = 0.5 * (r+(tan(a_i-a_t)/tan(a_i+a_t))^2);
end
xi = rand();
if xi <= r
% internally reflected
obj.uz = -obj.uz;
else
% transmitted
if dir == -1
obj.layer = obj.layer - 1;
else
obj.layer = obj.layer + 1;
end
if obj.layer == 1 || obj.layer == length(ly_ls)
obj.dead = true;
else
obj.ux = obj.ux*n_i/n_t;
obj.uy = obj.uy*n_i/n_t;
obj.uz = sign(obj.uz)*cos(a_t);
end
end
end
function obj = terminate(obj, m)
xi = rand();
if xi <= 1/m
obj.w = m*obj.w;
else
obj.w = 0;
obj.dead = true;
end
end
end
end