diff options
| -rw-r--r-- | face_beam_interaction.m | 35 |
1 files changed, 21 insertions, 14 deletions
diff --git a/face_beam_interaction.m b/face_beam_interaction.m index f7fbb2b..a94b2f2 100644 --- a/face_beam_interaction.m +++ b/face_beam_interaction.m @@ -4,17 +4,22 @@ function [is_face_hit, hit_position, hit_distance, new_beams] = face_beam_inter % 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 ibeam starts from face then its index is here + % 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 - index of refraction on the left hand side + % face.n_left - indexes of refraction on the left hand side % with respect to 1st -> 2nd vertex direction - % face.n_right - index of refraction on the right hand side + % [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); @@ -60,15 +65,15 @@ function [is_face_hit, hit_position, hit_distance, new_beams] = face_beam_inter % 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; - n2=face.n_right; + 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; - n2=face.n_left; + 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; @@ -93,30 +98,32 @@ function [is_face_hit, hit_position, hit_distance, new_beams] = face_beam_inter % angle of the normal with respect to horizon theta_normal = atan2(nf(2), nf(1)); - % reflected beam direction + %% 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; - reflectivity = .1; + beam_reflected.polarization = polarization; beam_reflected.intensity = beam.intensity * reflectivity; new_beams{1} = beam_reflected; - % refracted beam direction + %% refracted beam direction % refracted angle with respect to normal - sin_theta_refracted_rel2normal = n1/n2*sin(theta_i); - if ( abs(sin_theta_refracted_rel2normal) >=1 ) + if ( transmission == 0 ) % total internal reflection else % beam refracts - theta_refracted_rel2normal = asin( sin_theta_refracted_rel2normal ); if (are_nf_and_k_in_the_same_plane) theta_refracted = theta_normal + theta_refracted_rel2normal; else @@ -126,7 +133,7 @@ function [is_face_hit, hit_position, hit_distance, new_beams] = face_beam_inter beam_refracted.origin = hit_position; beam_refracted.k = [cos(theta_refracted), sin(theta_refracted)]; beam_refracted.face=closest_face_index; - transmission = .1; + beam_refracted.polarization = polarization; beam_refracted.intensity = beam.intensity * transmission; new_beams{2} = beam_refracted; end |