arbitrary beam tracing

6 files changed, 251 insertions, 1 deletions

diff --git a/beam2face_distance.m b/beam2face_distance.m
new file mode 100644
index 0000000..ba3d16b
--- /dev/null
+++ b/beam2face_distance.m
@@ -0,0 +1,36 @@
+function [hit_distance, hit_position, is_face_hit] = beam2face_distance(beam,face)
+	%% return  distance to face if beam hits it or infinity otherwise
+	% for beams and faces structure description see face_beam_interaction.m
+
+	k = beam.k;
+	% face direction or kf
+	kf=face.vertex2 - face.vertex1; % not a unit vector
+
+	%% let's find the intersection of lines passing through face and light beam
+	% we introduce parameters tb and tf 
+	% which define distance (in arb. units) along k vectors
+	% we are solving 
+	% xb_o+ k_x*tb = v1_x+ kf_x*tf
+	% yb_o+ k_y*tb = v1_y+ kf_y*tf
+	% A*t = B
+	A = [ k(1) , -kf(1); k(2), -kf(2)];
+	B=  [ face.vertex1(1)-beam.origin(1); face.vertex1(2)-beam.origin(2) ];
+	t = A\B; tb=t(1); tf=t(2);
+
+	% beam intersects face only if 0<=tf<=1 and tb>=0
+	if ( (0 <= tf) && (tf<=1) && tb >=0 )
+		% intersection,  beam hits the face
+		is_face_hit = true;
+	else
+		% beam misses the face
+		is_face_hit = false;
+		hit_position = [NA, NA];
+		hit_distance = Inf;
+		return;
+	end
+
+	% calculating hit position
+	hit_position=face.vertex1+kf*tf;
+	hit_distance = tb; % distance along the light beam
+end
+
diff --git a/beam_trace.m b/beam_trace.m
new file mode 100644
index 0000000..401391a
--- /dev/null
+++ b/beam_trace.m
@@ -0,0 +1,31 @@
+function ret = beam_trace(beams, faces, borders)
+	%% trace beam and all it reflections and refractions on faces
+	% for beams and faces structure description see face_beam_interaction.m
+	% border similar to faces cell array which enclose the beam
+	% i.e. it does not leave the polygon consisting of border faces
+
+	Nbeams=size(beams)(2);
+	if ( Nbeams == 0 )
+		% no more beam to trace
+		return;
+	end
+
+	% trace each beam in beams cell array
+	for i=1:Nbeams
+		beam=beams{i};
+		[is_face_hit, hit_position, hit_distance, new_beams] =  face_beam_interaction(beam, faces);
+		if (is_face_hit)
+			beam_trace(new_beams, faces, borders );
+			continue;
+		end
+
+		% beam does not hit face but it should stop and borders
+		beam.face=NA;
+		[is_face_hit, hit_position, hit_distance, new_beams] =  face_beam_interaction(beam, borders);
+		if (!is_face_hit)
+			error('borders are badly defined, the beam misses them');
+		end
+	end
+end
+
+
diff --git a/coupling_angles.m b/coupling_angles.m
index df7c5f1..ca7dcf0 100644
--- a/coupling_angles.m
+++ b/coupling_angles.m
@@ -144,4 +144,5 @@ text(0,1.2, str);
 axis([-1,1,-0.5,1.5],'equal');
 
 %% output of the plot to the file
-print('prism_disk_coupling.pdf');
+print('prism_disk_coupling.eps','-depsc2');
+print('prism_disk_coupling.png');
diff --git a/example.m b/example.m
new file mode 100644
index 0000000..0b1e855
--- /dev/null
+++ b/example.m
@@ -0,0 +1,30 @@
+beam.k=[1,-.5];
+beam.origin=[0,0];
+beam.face=NA;
+
+face.vertex1=[1.1,-1];
+face.vertex2=[2.7,1];
+face.n_right = 1;
+face.n_left= 2.5;
+
+face2.vertex1=[-.1,-3];
+face2.vertex2=[-.1,3];
+face2.n_right = 1;
+face2.n_left= 1;
+
+face3.vertex1=[2.1,-3];
+face3.vertex2=[2.1,3];
+face3.n_right = 1.5;
+face3.n_left= 1;
+
+faces={face,face2,face3};
+
+[is_face_hit, hit_position, hit_distance, new_beams] =  face_beam_interaction(beam, faces);
+beam_reflected=new_beams{1};
+beam_refracted=new_beams{2};
+
+
+[is_face_hit, hit_position, hit_distance, new_beams] =  face_beam_interaction(beam_reflected, faces);
+
+[is_face_hit, hit_position, hit_distance, new_beams] =  face_beam_interaction(beam_refracted, faces);
+
diff --git a/example_beam_trace.m b/example_beam_trace.m
new file mode 100644
index 0000000..ada35c9
--- /dev/null
+++ b/example_beam_trace.m
@@ -0,0 +1,45 @@
+beam.k=[1,-.5];
+beam.origin=[0,0];
+beam.face=NA;
+
+face.vertex1=[1.1,-1];
+face.vertex2=[2.7,1];
+face.n_right = 1;
+face.n_left= 2.5;
+
+face2.vertex1=[-.1,-3];
+face2.vertex2=[-.1,3];
+face2.n_right = 1;
+face2.n_left= 1;
+
+face3.vertex1=[2.1,-3];
+face3.vertex2=[2.1,3];
+face3.n_right = 1.5;
+face3.n_left= 1;
+
+faces={face,face2,face3};
+
+
+border1.vertex1=[-5,-5];
+border1.vertex2=[5,-5];
+border1.n_right = 1.;
+border1.n_left= 1;
+
+border2.vertex1=[5,-5];
+border2.vertex2=[5,5];
+border2.n_right = 1.;
+border2.n_left= 1;
+
+border3.vertex1=[5,5];
+border3.vertex2=[-5,5];
+border3.n_right = 1.;
+border3.n_left= 1;
+
+border4.vertex1=[-5,5];
+border4.vertex2=[-5,-5];
+border4.n_right = 1.;
+border4.n_left= 1;
+
+borders={border1,border2, border3, border4};
+beams={beam};
+beam_trace(beams, faces, borders )
diff --git a/face_beam_interaction.m b/face_beam_interaction.m
new file mode 100644
index 0000000..4294c11
--- /dev/null
+++ b/face_beam_interaction.m
@@ -0,0 +1,107 @@
+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 ibeam starts from face then its index is here
+	% 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 
+	%                 with respect to 1st -> 2nd vertex direction
+	%  face.n_right - index of refraction on the right hand side 
+
+
+	k=beam.k;
+
+	%% 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
+	t=linspace(0,1);
+	x=face.vertex1(1)+kf(1)*t;
+	y=face.vertex1(2)+kf(2)*t;
+	plot(x,y,'k-');
+	t=linspace(0,hit_distance);
+	% draw beam
+	x=beam.origin(1)+k(1)*t;
+	y=beam.origin(2)+k(2)*t;
+	plot(x,y,'r-');
+
+	% 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;
+	       n2=face.n_right;
+	else	       
+	       % beam coming from the right
+	       n1=face.n_right;
+	       n2=face.n_left;
+	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);
+	
+	% reflected beam direction
+	theta_normal = atan2(nf(2), nf(1));
+	theta_reflected = theta_normal + pi - theta_i;
+	
+	beam_reflected.origin = hit_position;
+	beam_reflected.k = [cos(theta_reflected), sin(theta_reflected)];
+	beam_reflected.face=closest_face_index;
+	new_beams{1} = beam_reflected;
+
+
+	% 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 ) 
+		% total internal reflection
+	else
+		% beam refracts
+		theta_refracted_rel2normal = asin( sin_theta_refracted_rel2normal );
+		theta_refracted = theta_normal + theta_refracted_rel2normal;
+
+		beam_refracted.origin = hit_position;
+		beam_refracted.k = [cos(theta_refracted), sin(theta_refracted)];
+		beam_refracted.face=closest_face_index;
+		new_beams{2} = beam_refracted;
+	end
+end
+
+
+