proper reflection and refraction

1 files changed, 32 insertions, 3 deletions

diff --git a/face_beam_interaction.m b/face_beam_interaction.m
index 4294c11..f7fbb2b 100644
--- a/face_beam_interaction.m
+++ b/face_beam_interaction.m
@@ -62,10 +62,16 @@ function [is_face_hit, hit_position, hit_distance, new_beams] =  face_beam_inter
 	       % beam coming from the left
 	       n1=face.n_left;
 	       n2=face.n_right;
+	       % 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;
+	       % 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 
@@ -75,14 +81,31 @@ function [is_face_hit, hit_position, hit_distance, new_beams] =  face_beam_inter
 	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
 	
-	% reflected beam direction
+	% angle of the normal with respect to horizon
 	theta_normal = atan2(nf(2), nf(1));
-	theta_reflected = theta_normal + pi - theta_i;
+
+	% 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
+
 	
 	beam_reflected.origin = hit_position;
 	beam_reflected.k = [cos(theta_reflected), sin(theta_reflected)];
 	beam_reflected.face=closest_face_index;
+	reflectivity = .1;
+	beam_reflected.intensity = beam.intensity * reflectivity;
 	new_beams{1} = beam_reflected;
 
 
@@ -94,11 +117,17 @@ function [is_face_hit, hit_position, hit_distance, new_beams] =  face_beam_inter
 	else
 		% beam refracts
 		theta_refracted_rel2normal = asin( sin_theta_refracted_rel2normal );
-		theta_refracted = theta_normal + theta_refracted_rel2normal;
+		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;
+		transmission = .1;
+		beam_refracted.intensity = beam.intensity * transmission;
 		new_beams{2} = beam_refracted;
 	end
 end