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
|
function [is_face_hit, hit_position, hit_distance, new_beams] = face_beam_interaction(beam, faces)
%% calculates refracted and reflected beam after interaction with a face
% beam - structure defining the light beam
% beam.origin - array [x,y] origin/soutce of the light beam
% beam.k - k vector i.e. direction [kx,ky]
% beam.intensity - intensity of the beam
% beam.face - if beam starts from face then its index is here
% beam.polarization - can be either 1 for 's' or 2 for 'p'
% faces cell array of face structures
% face - structure definiong the beam
% face.vertex1 - [x,y] of the 1st point/vertex of the face
% face.vertex2 - [x,y] of the 2nd point/vertex of the face
% face.n_left - indexes of refraction on the left hand side
% with respect to 1st -> 2nd vertex direction
% [ns,np] - for s and p polarizations
% face.n_right - indexes of refraction on the right hand side
% [ns,np] - for s and p polarizations
k=beam.k;
polarization=beam.polarization;
%% we go over all faces to find the closest which beam hits
Nfaces=size(faces)(2);
hit_distance=Inf;
is_face_hit = false;
hit_position = [NA, NA];
closest_face_index=NA;
for i=1:Nfaces
if ( beam.face == i) continue; end
face=faces{i};
[hit_distance_tmp, hit_position_tmp, is_face_hit_tmp] = beam2face_distance(beam,face);
if ( hit_distance_tmp < hit_distance )
% this is the closer face
is_face_hit=is_face_hit_tmp;
hit_position=hit_position_tmp;
hit_distance=hit_distance_tmp;
closest_face_index=i;
end
end
if (!is_face_hit)
new_beams={};
return;
end
%% closest face
face=faces{closest_face_index};
kf=face.vertex2 - face.vertex1; % not a unit vector
hold on;
% draw face
t=linspace(0,1);
x=face.vertex1(1)+kf(1)*t;
y=face.vertex1(2)+kf(2)*t;
plot(x,y,'k-');
t=linspace(0,hit_distance);
% draw beam
x=beam.origin(1)+k(1)*t;
y=beam.origin(2)+k(2)*t;
plot(x,y,'r-');
% find is beam arriving from left or right. I will use vectors cross product property.
% if z component of the product 'k x kf' is positive then beam arrives from the left
if ( ( k(1)*kf(2)-k(2)*kf(1) ) > 0 )
% beam coming from the left
n1=face.n_left(polarization);
n2=face.n_right(polarization);
% this means that notmal and beam k vector point to the same half plane
% relative to the face
are_nf_and_k_in_the_same_plane=true;
else
% beam coming from the right
n1=face.n_right(polarization);
n2=face.n_left(polarization);
% this means that notmal and beam k vector point to the different half plane
% relative to the face
are_nf_and_k_in_the_same_plane=false;
end
% normal vector to the face, looks to the left of it
nf=[ kf(2), -kf(1) ] / norm(kf);
% incidence angle calculation
cos_theta_i = dot(k, nf) / (norm(k)*norm(nf));
sin_theta_i = - ( k(1)*nf(2)-k(2)*nf(1) ) / (norm(k)*norm(nf));
% positive angle to the right from normal before incidence to the face
theta_i = atan2(sin_theta_i, cos_theta_i);
% we need to make sure that angle of incidence belong to [-pi/2,pi/2] interval
if( theta_i > pi/2)
theta_i=pi-theta_i;
end
if( theta_i < -pi/2)
theta_i=-pi+theta_i;
end
% angle of the normal with respect to horizon
theta_normal = atan2(nf(2), nf(1));
%% reflected beam direction
if (are_nf_and_k_in_the_same_plane)
theta_reflected = theta_normal - theta_i + pi;
else
theta_reflected = theta_normal + theta_i;
end
%% coefficients of reflection and transmission for given polarization
[R,theta_refracted_rel2normal]=fresnel_reflection(n1, n2, theta_i);
reflectivity = R(polarization);
transmission = 1 - reflectivity;
beam_reflected.origin = hit_position;
beam_reflected.k = [cos(theta_reflected), sin(theta_reflected)];
beam_reflected.face=closest_face_index;
beam_reflected.polarization = polarization;
beam_reflected.intensity = beam.intensity * reflectivity;
new_beams{1} = beam_reflected;
%% refracted beam direction
% refracted angle with respect to normal
if ( transmission == 0 )
% total internal reflection
else
% beam refracts
if (are_nf_and_k_in_the_same_plane)
theta_refracted = theta_normal + theta_refracted_rel2normal;
else
theta_refracted = theta_normal - theta_refracted_rel2normal + pi;
end
beam_refracted.origin = hit_position;
beam_refracted.k = [cos(theta_refracted), sin(theta_refracted)];
beam_refracted.face=closest_face_index;
beam_refracted.polarization = polarization;
beam_refracted.intensity = beam.intensity * transmission;
new_beams{2} = beam_refracted;
end
end
|
