moved octave using code to separate folder

1 files changed, 151 insertions, 0 deletions

diff --git a/beam_tracing/octave/face_beam_interaction.m b/beam_tracing/octave/face_beam_interaction.m
new file mode 100644
index 0000000..3bc3e4a
--- /dev/null
+++ b/beam_tracing/octave/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
+
+
+