diff options
Diffstat (limited to 'beam_tracing/face_beam_interaction.m')
| -rw-r--r-- | beam_tracing/face_beam_interaction.m | 151 |
1 files changed, 151 insertions, 0 deletions
diff --git a/beam_tracing/face_beam_interaction.m b/beam_tracing/face_beam_interaction.m new file mode 100644 index 0000000..3bc3e4a --- /dev/null +++ b/beam_tracing/face_beam_interaction.m @@ -0,0 +1,151 @@ +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' + % beam.status - could be 'reflected', 'refracted', 'incoming' + % 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 + line([face.vertex1(1), face.vertex2(1)], [face.vertex1(2), face.vertex2(2)] , 'linestyle', '-', 'color', [0,0,0] ); + % draw beam + Ncolors=256; + linewidth=2; + if ( polarization == 1) + % s-polarization + line_base_color=[1,0,0]; % RGB - red + else + % p-polarization + line_base_color=[0,0,1]; % RGB - blue + end + %plot(x,y,line_str); + line_color = color_gradient( beam.intensity*Ncolors, line_base_color, Ncolors); + line([beam.origin(1), hit_position(1)], [beam.origin(2), hit_position(2)], 'color', line_color, 'linewidth', 3, 'linestyle', '-'); + + % 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; + beam_reflected.status = 'reflected'; + 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; + beam_refracted.status = 'refracted'; + new_beams{2} = beam_refracted; + end +end + + + |