reflection coefficient according to Fresnel equation

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